\jourtitle

Mathematics in Computer Science

^†^†thanks: This is a preprint of a paper whose final and definite form is with Math. Comput. Sci., ISSN 1661-8270 (print), ISSN 1661-8289 (electronic), available at [https://www.springer.com/journal/11786]. Submitted 1-Jun-2019; Revised 5-Feb-2020; Accepted 6-Feb-2020.

Global stability of a Caputo fractional SIRS model
with general incidence rate

Moulay Rchid Sidi Ammi¹ Mostafa Tahiri² Delfim F. M. Torres³
ORCID ¹Department of Mathematics, AMNEA Group, Faculty of Sciences and Techniques,
Moulay Ismail University, B. P. 509, Errachidia, Morocco
[email protected], http://orcid.org/0000-0002-4488-9070
²Department of Mathematics, AMNEA Group, Faculty of Sciences and Techniques,
Moulay Ismail University, B. P. 509, Errachidia, Morocco
[email protected]
ORCID ³Center for Research and Development in Mathematics and Applications (CIDMA),
Department of Mathematics, University of Aveiro, 3810-193 Aveiro, Portugal
[email protected], http://orcid.org/0000-0001-8641-2505

(Received: 1 June 2019 / Revised: 5 February 2020 / Accepted: 6 February 2020)

Abstract

We introduce a fractional order SIRS model with non-linear incidence rate. Existence of a unique positive solution to the model is proved. Stability analysis of the disease free equilibrium and positive fixed points are investigated. Finally, a numerical example is presented.

keywords:

Epidemiology and mathematical modeling and fractional calculus and Caputo derivatives and equilibrium and stability

1991 Mathematics Subject Classification:

26A33 and 34A08 and 92D30

1. Introduction

Fractional differential equations (FDEs) are generalizations of classical differential equations, where the integer-order derivative is replaced by a non-integer one. There has been a significant development in FDEs in recent years due its applicability in different fields of science and engineering [1, 7]. In particular, fractional derivatives are used to describe viscoelastic properties of many polymeric materials [25], in diffusion equations [49], in mechanics [50], and decision-making problems [54].

It is worthwhile to mention that fractional derivatives are non-local operators and thus may be more suitable for modelling systems dependent on past history (memory). More precisely, the fractional derivative of a given function does not depend only on its current state, but also on previous historical states [42, 48].

In epidemiology, most mathematical models descend from the classical SIR model of Kermack and McKendrick, established in 1927 [27, 43]. Recently, fractional derivatives have been used to describe epidemiological models and, in some cases, they have proven to be more accurate when compared to the classical ones [10, 12, 45]. Different models described by fractional derivatives are available in the literature, like the SIR model [17, 18, 39], the SIR model with vaccination [46], the SIRC model [19], and the SEIR model [40].

Since the fractional order can be any positive real $\alpha$ , one can choose the one that better fits available data [3]. Therefore, we can adjust the model to real data and, by doing so, better predict the future evolution of the disease taking into account its past and present [45, 55]. Moreover, virus propagation is typically discontinuous, something the classical differential models cannot describe in a proper way. In contrast, fractional systems deal naturally with such discontinuous properties [15, 47].

The virus propagation is similar to heat transmission or moistness penetrability in a porous medium, which can be exactly modelled by fractional calculus [31, 58]. The authors in [22, 23] give a geometrical description of fractional calculus, concluding that the fractional order can be related with the fractal dimension. The relationship between fractal dimension and fractional calculus has been obtained by several different authors: see [44, 56] and references therein. The fractional complex transform [24, 29] is an approximate transform of a fractal space (time) to a continuous one, and it is now widely used in fractional calculus [8, 53, 57].

There are several definitions of fractional derivatives [7, 33]. In this paper, we choose to work with the celebrated Caputo fractional derivatives. One of the main advantages of such derivatives is allowing us to consider classical initial conditions in the formulation of the problem. Also, the Caputo fractional derivatives of a constant are zero. Such properties of the Caputo fractional derivatives are not true for other fractional operators, for example for the Riemann–Liouville derivatives [42, 48].

Most non-linear fractional differential equations do not have analytic solutions [6, 28]. Therefore, approximations and numerical techniques must be used [16, 37, 47]. The decomposition method [36] and the variational iteration method [21, 35] are relatively new approaches to provide an analytical approximate solution to linear and non-linear problems. For a simple algorithm, based on fractional Euler’s method, to numerically solve non-linear fractional differential equations, in a direct way, without using linearisation, perturbations, or restrictive assumptions, see [37].

Here we propose a fractional SIRS model with the spread of the disease being described by a system of non-linear fractional order differential equations as follows:

\left\{\begin{array}[]{lll}D^{\alpha}S(t)=\Lambda-\mu S(t)-\dfrac{\beta S(t)I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}+\lambda R(t),\\[8.61108pt] D^{\alpha}I(t)=\dfrac{\beta S(t)I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-(\mu+r)I(t),\\[8.61108pt] D^{\alpha}R(t)=rI(t)-(\mu+\lambda)R(t),\end{array}\right.

(1)

where $D^{\alpha}$ denotes the (left) Caputo fractional derivative of order $\alpha$ , $0<\alpha\leq 1$ . The model considers a population that is divided into three subgroups: susceptible $S(t)$ , infective $I(t)$ , and recovered $R(t)$ individuals at time $t$ . The positive constants $\Lambda$ , $\beta$ , $\mu$ , and $r$ , are the recruitment rate of the population, the infection rate, the natural death rate, and the recovery rate of the infective individuals, respectively. The rate that recovered individuals lose immunity and return to the susceptible class is $\lambda$ . While contacting with infected individuals, the susceptible become infected at the incidence rate $\beta SI/(1+k_{1}S+k_{2}I+k_{3}SI)$ , with $k_{1}$ , $k_{2}$ , and $k_{3}$ non-negative constants [20]. This incidence function generalizes several types of incidence rates, for example, the traditional bilinear incidence rate, the saturated incidence rate, the Beddington–DeAngelis functional response proposed in [9, 14], and the Crowley–Martin functional response introduced in [13]. For the advantages of using a general incidence rate, see [32, 34]. For other ways to fractionalize a classical system of differential equations, see the discussion in [4, 11].

The paper is organized as follows. In Section 2, we recall necessary definitions and properties from fractional calculus. Our results begin with Section 3, where we show the existence and uniqueness of positive solution. In Section 4, we study the existence of equilibria and their local stability. The global stability is investigated in Section 5. In order to illustrate our theoretical results, numerical simulations of the model are given in Section 6. We end with Section 7 of conclusions and future perspectives.

2. Basic results of fractional calculus

There are many good books on fractional calculus. For a gentle introduction, we refer the reader to [42]. For an encyclopedic treaty, see [48].

Definition 1 (See [42]).

The Riemann–Liouville fractional integral of order $\alpha>0$ of a function $f:\mathbb{R}^{+}\rightarrow\mathbb{R}$ is given by

I^{\alpha}f(x)=\dfrac{1}{\Gamma(\alpha)}\int_{0}^{x}(x-t)^{\alpha-1}f(t)dt,

where $\Gamma(\alpha)=\int_{0}^{\infty}t^{\alpha-1}e^{-t}dt$ is the Euler Gamma function.

Definition 2 (See [42]).

Let $\alpha>0$ , $n=[\alpha]+1$ , $n-1<\alpha\leq n$ , where $[\alpha]$ denotes the integer part of $\alpha$ . The Caputo fractional derivative of order $\alpha$ for a function $f\in C^{n}([0,+\infty),\mathbb{R})$ is defined by

D^{\alpha}f(u)=I^{n-\alpha}D^{n}f(u)=\dfrac{1}{\Gamma(n-\alpha)}\int_{0}^{u}\dfrac{f^{(n)}(s)}{(u-s)^{\alpha+1-n}}ds,\quad u>0,

where $D$ is the usual differential operator, that is, $D=\dfrac{d}{du}$ . In particular, when $0<\alpha\leq 1$ , one has

D^{\alpha}f(u)=\dfrac{1}{\Gamma(1-\alpha)}\int_{0}^{u}\dfrac{f^{\prime}(s)}{(u-s)^{\alpha}}ds.

Next we recall the definition of the Mittag–Leffler function of parameter $\alpha$ , which is a generalization of the exponential function.

Definition 3 (See [42]).

Let $\alpha>0$ . The function $E_{\alpha}$ defined by $E_{\alpha}(z)=\displaystyle\sum_{j=0}^{\infty}\dfrac{z^{j}}{\Gamma(\alpha j+1)}$ is called the Mittag–Leffler function of parameter $\alpha$ .

Let $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ with $n\geq 1$ . Consider the fractional order system

\left\{\begin{array}[]{lll}D^{\alpha}x(t)=f(x),\\ x(0)=x_{0},\end{array}\right.

(2)

where $0<\alpha\leq 1$ and $x_{0}\in\mathbb{R}^{n}$ . The following lemma, which is a direct corollary from the main result of [30], gives global existence of solution to system (2).

Lemma 4 (See [30]).

Assume that $f$ satisfies the following conditions:

(1)

$f(x)$ and $\dfrac{\partial f}{\partial x}$ are continuous for all $x\in\mathbb{R}^{n}$ ;
(2)

$\|f(x)\|\leq\omega+\lambda\|x\|$ $\forall x\in\mathbb{R}^{n}$ , where $\omega$ and $\lambda$ are two positive constants.

Then, system (2) has a unique solution on $[0,+\infty)$ .

3. Existence and uniqueness of positive solution

Denote $\mathbb{R}_{+}^{3}=\{X\in\mathbb{R}^{3}:X\geq 0\}$ and let $X(t)=(S(t),I(t),R(t))^{T}$ . Then system (1) can be reformulated as follows: $D^{\alpha}X(t)=F(X(t))$ , where

F(X)=\left(\begin{array}[]{ccc}\Lambda-\mu S-\dfrac{\beta SI}{1+k_{1}S+k_{2}I+k_{3}SI}+\lambda R\\[8.61108pt] \dfrac{\beta SI}{1+k_{1}S+k_{2}I+k_{3}SI}-(\mu+r)I\\[8.61108pt] rI-(\mu+\lambda)R\end{array}\right).

(3)

For biological reasons, we consider system (1) with the following initial conditions:

S(0)\geq 0,\quad I(0)\geq 0,\quad R(0)\geq 0.

(4)

To prove the main theorem of this section, i.e., Theorem 7, we need the following generalized mean value theorem and its corollary.

Lemma 5 (Generalized Mean Value Theorem [38]).

Suppose that $f\in C[0,b]$ and $D^{\alpha}f\in C(0,b]$ , $0<\alpha\leq 1$ . Then, one has

f(x)=f(0)+\dfrac{1}{\Gamma(\alpha)}(D^{\alpha}f)(\xi)x^{\alpha}

with $0\leq\xi\leq x$ , $\forall x\in(0,b]$ .

Corollary 6 (See [38]).

Suppose that $f\in C[0,b]$ and $D^{\alpha}f\in C(0,b]$ for $0<\alpha\leq 1$ . If $D^{\alpha}f(x)\geq 0$ $\forall x\in(0,b)$ , then $f(x)$ is non-decreasing for each $x\in[0,b]$ . If $D^{\alpha}f(x)\leq 0$ $\forall\ x\in(0,b)$ , then $f(x)$ is non-increasing for each $x\in[0,b]$ .

Theorem 7.

There is a unique solution for (1) satisfying (4) for $t\geq 0$ and the solution will remain in $\mathbb{R}_{+}^{3}$ for all $t\geq 0$ . Moreover, $N(t)\leq N(0)+\dfrac{\Lambda}{\mu}$ , where $N(t)=S(t)+I(t)+R(t)$ .

Proof.

Since the vector function $F$ (3) satisfies the first condition of Lemma 4, we only need to prove the second one. Denote

	$\displaystyle\varepsilon=\left(\begin{array}[]{ccc}\Lambda\\[8.61108pt] 0\\[8.61108pt] 0\end{array}\right),\quad A_{1}=\left(\begin{array}[]{ccc}-\mu&0&\lambda\\[8.61108pt] 0&-(\mu+r)&0\\[8.61108pt] 0&r&-(\mu+\lambda)\end{array}\right),\quad A_{2}=\left(\begin{array}[]{ccc}0&-\beta/k_{1}&0\\[8.61108pt] 0&\beta/k_{1}&0\\[8.61108pt] 0&0&0\end{array}\right),$
	$\displaystyle A_{3}=\left(\begin{array}[]{ccc}-\beta/k_{2}&0&0\\[8.61108pt] \beta/k_{2}&0&0\\[8.61108pt] 0&0&0\end{array}\right),\quad A_{4}=\left(\begin{array}[]{ccc}-\beta/k_{3}\\[8.61108pt] \beta/k_{3}\\[8.61108pt] 0\end{array}\right),\quad A_{5}=\left(\begin{array}[]{ccc}-\beta&0&0\\[8.61108pt] \beta&0&0\\[8.61108pt] 0&0&0\end{array}\right).$

We discuss four cases, as follows.
Case 1. If $k_{1}\neq 0$ , then we have

F(X)=\varepsilon+A_{1}X+\dfrac{k_{1}S}{1+k_{1}S+k_{2}I+k_{3}SI}A_{2}X.

Therefore,

\|F(X)\|\leq\|\varepsilon\|+\|A_{1}X\|+\|A_{2}X\|=\|\varepsilon\|+(\|A_{1}\|+\|A_{2}\|)\|X\|.

Case 2. If $k_{2}\neq 0$ , then

F(X)=\varepsilon+A_{1}X+\dfrac{k_{2}I}{1+k_{1}S+k_{2}I+k_{3}SI}A_{3}X.

Thus,

\|F(X)\|\leq\|\varepsilon\|+\|A_{1}X\|+\|A_{3}X\|=\|\varepsilon\|+(\|A_{1}\|+\|A_{3}\|)\|X\|.

Case 3. If $k_{3}\neq 0$ , then one has

F(X)=\varepsilon+A_{1}X+\dfrac{k_{3}SI}{1+k_{1}S+k_{2}I+k_{3}SI}A_{4}.

We conclude that

\|F(X)\|\leq\|\varepsilon\|+\|A_{4}\|+\|A_{1}\|\|X\|.

Case 4. If $k_{1}=k_{2}=k_{3}=0$ , then we obtain

F(X)=\varepsilon+A_{1}X+IA_{5}X.

It follows that

\|F(X)\|\leq\|\varepsilon\|+(\|A_{1}\|+\|I\|\|A_{5}\|)\|X\|.

By Lemma 4, it follows that (1) subject to (4) has a unique solution. Now we prove the non-negativity of the solution. Observe first that

D^{\alpha}S/_{S=0}=\Lambda+\lambda R,

D^{\alpha}I/_{I=0}=0,

and

D^{\alpha}R/_{R=0}=rI.

We can prove that the solution of (1) remains non-negative for all $t\geq 0$ by proceeding in a similar way as in [5, Theorem 2], that is, by considering an auxiliary system of fractional differential equations and by reductio ad absurdum. For that we use Corollary 6, which is a consequence of Lemma 5, to get a contradiction and then arriving to the intended conclusion by [5, Lemma 1]. Finally, we establish the boundedness of solution. By summing all the equations of system (1), we obtain that

D^{\alpha}N=\Lambda-\mu N.

Solving this equality, we get

N(t)\leq N(0)E_{\alpha}(-\mu t^{\alpha})+\dfrac{\Lambda}{\mu}(1-E_{\alpha}(-\mu t^{\alpha})).

Because $0\leq E_{\alpha}(-\mu t^{\alpha})\leq 1$ , we have $N(t)\leq N(0)+\dfrac{\Lambda}{\mu}$ . This completes the proof. ∎

4. Local stability

In this section, we firstly discuss the existence of equilibria for model (1). Let

R_{0}=\dfrac{\beta\Lambda}{(\mu+\Lambda k_{1})(\mu+r)}.

We prove that model (1) has two possible equilibria.

Theorem 8.

(i)

There is always a disease-free equilibrium $E_{0}=(S_{0},0,0)$ , where $S_{0}=\dfrac{\Lambda}{\mu}$ .

(ii)

If $R_{0}>1$ , then there exists a unique endemic equilibrium

E^{*}=\left(S^{*},\dfrac{(\mu+\lambda)(\Lambda-\mu S^{*})}{c},\dfrac{r(\Lambda-\mu S^{*})}{c}\right),

where

	$\displaystyle c=a(\mu+\lambda)-\lambda r,\quad a=\mu+r,$
	$\displaystyle S^{*}=\dfrac{k_{3}a\Lambda(\mu+\lambda)+k_{1}ac-\beta c-k_{2}a\mu(\mu+\lambda)+\sqrt{\Delta}}{2k_{3}a\mu(\mu+\lambda)},$
	$\displaystyle\Delta=\left(\beta c-k_{1}ac-k_{3}a\Lambda(\mu+\lambda)+k_{2}a\mu(\mu+\lambda)\right)^{2}+4k_{3}a\mu\left(ac+k_{2}a\Lambda(\mu+\lambda)\right).$

Proof.

(i) By direct calculation, we have that $E_{0}$ is the unique steady state of system (1).
(ii) To find the other equilibrium, we solve the system

F(X)=0.

(5)

Let

f(S,I)=\dfrac{\beta S}{1+k_{1}S+k_{2}I+k_{3}SI}.

From system (5), we obtain $I=\dfrac{(\mu+\lambda)(\Lambda-\mu S)}{c}$ , $R=\dfrac{r(\Lambda-\mu S)}{c}$ , and $f\left(S,\dfrac{(\mu+\lambda)(\Lambda-\mu S^{*})}{c}\right)=a$ . Since $c=\mu^{2}+\mu\lambda+r\mu>0$ , we get $I\geq 0$ if $S\leq\dfrac{\Lambda}{\mu}$ . Now we consider function

g(S)=f\left(S,\dfrac{(\mu+\lambda)(\Lambda-\mu S^{*})}{c}-a\right)

defined on the interval $\left[0,\dfrac{\Lambda}{\mu}\right]$ . One has

\dfrac{\partial f}{\partial S}=\dfrac{\beta(1+k_{2}I)}{(1+k_{1}S+k_{2}I+k_{3}SI)^{2}}>0

and

\dfrac{\partial f}{\partial I}=\dfrac{-\beta S(k_{2}+k_{3}S)}{(1+k_{1}S+k_{2}I+k_{3}SI)^{2}}<0.

Then $g^{\prime}(S)>0$ , which implies that $g$ is strictly increasing on $\left[0,\dfrac{\Lambda}{\mu}\right]$ . Hence, if $R_{0}>1$ , the system admits a unique endemic equilibrium $E^{*}=(S^{*},I^{*},R^{*})$ with $S^{*}\in\left(0,\dfrac{\Lambda}{\mu}\right)$ , $I^{*}>0$ , and $R^{*}>0$ . This completes the proof. ∎

Next, we study the local stability of the disease-free equilibrium $E_{0}$ and the endemic equilibrium $E^{*}$ . The Jacobian matrix of system (1) at any equilibrium $\bar{E}=(\bar{S},\bar{I},\bar{R})$ is given by

J_{\bar{E}}=\left(\begin{array}[]{ccc}-\mu-\dfrac{\beta\bar{I}(1+k_{2}\bar{I})}{(1+k_{1}\bar{S}+k_{2}\bar{I}+k_{3}\bar{S}\bar{I})^{2}}&-\dfrac{\beta\bar{S}(1+k_{1}\bar{S})}{(1+k_{1}\bar{S}+k_{2}\bar{I}+k_{3}\bar{S}\bar{I})^{2}}&\lambda\\[8.61108pt] \dfrac{\beta\bar{I}(1+k_{2}\bar{I})}{(1+k_{1}\bar{S}+k_{2}\bar{I}+k_{3}\bar{S}\bar{I})^{2}}&\dfrac{\beta\bar{S}(1+k_{1}\bar{S})}{(1+k_{1}\bar{S}+k_{2}\bar{I}+k_{3}\bar{S}\bar{I})^{2}}-a&0\\[8.61108pt] 0&r&-(\mu+\lambda)\end{array}\right).

(6)

We recall that a sufficient condition for the local stability of $\bar{E}$ is

|arg(\xi_{i})|>\dfrac{\alpha\pi}{2},\quad i=1,2,

(7)

where $\xi_{i}$ are the eigenvalues of $J_{\bar{E}}$ (see [41]). First, we establish the local stability of $E_{0}$ .

Theorem 9.

If $R_{0}<1$ , then the disease-free equilibrium $E_{0}$ is locally asymptotically stable.

Proof.

At $E_{0}$ , (6) becomes

J_{E_{0}}=\left(\begin{array}[]{ccc}-\mu&-\dfrac{\beta\Lambda}{\mu+k_{1}\Lambda}&\lambda\\[8.61108pt] 0&\dfrac{\beta\Lambda}{\mu+k_{1}\Lambda}-a&0\\[8.61108pt] 0&r&-(\mu+\lambda)\end{array}\right).

Hence, the eigenvalues of $J_{E_{0}}$ are $\xi_{1}=-\mu$ , $\xi_{2}=a(R_{0}-1)$ , and $\xi_{3}=-(\mu+\lambda)$ . Clearly, $\xi_{2}$ satisfies condition (7) if $R_{0}<1$ , since $\xi_{1}$ and $\xi_{3}$ are negative, proving the desired result. ∎

We now establish the local stability of $E^{*}$ .

Theorem 10.

If $R_{0}>1$ , then the endemic equilibrium $E^{*}$ is locally asymptotically stable.

Proof.

At equilibrium $E^{*}$ , the characteristic equation for the corresponding linearised system of model (1) is $\xi^{3}+a_{1}\xi^{2}+a_{2}\xi+a_{3}=0$ , where

a_{1}=a+2\mu+\lambda+\dfrac{\beta I^{*}(1+k_{2}I^{*})-\beta S^{*}(1+k_{1}S^{*})}{(1+k_{1}S^{*}+k_{2}I^{*}+k_{3}S^{*}I^{*})^{2}},

a_{2}=\mu a+(\mu+\lambda)\mu+(\mu+\lambda)a+\dfrac{(a+\mu+\lambda)\beta I^{*}(1+k_{2}I^{*})-(2\mu+\lambda)\beta S^{*}(1+k_{1}S^{*})}{(1+k_{1}S^{*}+k_{2}I^{*}+k_{3}S^{*}I^{*})^{2}},

and

a_{3}=(\mu+\lambda)\mu a+\dfrac{c\beta I^{*}(1+k_{2}I^{*})-(\mu+\lambda)\mu\beta S^{*}(1+k_{1}S^{*})}{(1+k_{1}S^{*}+k_{2}I^{*}+k_{3}S^{*}I^{*})^{2}}.

Let $D(f)$ denote the discriminant of polynomial $f(\xi)=\xi^{3}+a_{1}\xi^{2}+a_{2}\xi+a_{3}$ . Then,

D(f)=18a_{1}a_{2}a_{3}+(a_{1}a_{2})^{2}-4a_{3}a_{1}^{3}-4a_{2}^{3}-27a_{3}^{2}.

Suppose that $E^{*}$ exists in $\mathbb{R}_{+}^{3}$ . Based on [52], we have the following conclusions by using Routh–Hurwitz conditions:

(i)

if $a_{1}>0$ , $a_{3}>0$ , and $a_{1}a_{2}>a_{3}$ , then $E^{*}$ is locally asymptotically stable for all $\alpha\in(0,1]$ ;
(ii)

if $D(f)<0$ , $a_{1}\geq 0$ , $a_{2}\geq 0$ , $a_{3}>0$ , $a_{1}a_{2}<a_{3}$ , and $\alpha<\dfrac{2}{3}$ , then $E^{*}$ is locally asymptotically stable;
(iii)

if $D(f)<0$ , $a_{1}<0$ , $a_{2}<0$ , and $\alpha>\dfrac{2}{3}$ , then $E^{*}$ is unstable;
(iv)

if $D(f)<0$ , $a_{1}>0$ , $a_{2}>0$ , $a_{1}a_{2}=a_{3}$ , and $\alpha\in(0,1]$ , then $E^{*}$ is locally asymptotically stable.
(v)

if $D(f)<0$ , $a_{1}>0$ , $a_{3}=0$ , and $\alpha\in(0,1]$ , then $E^{*}$ is locally stable.

The proof is complete. ∎

5. Global stability

In this section, we investigate the global stability of both equilibria $E_{0}$ and ${E^{*}}$ .

Theorem 11.

The disease-free equilibrium $E_{0}$ is globally asymptotically stable whenever $R_{0}\leq 1$ .

Proof.

Consider the following Lyapunov function:

L_{0}(t)=\dfrac{S_{0}}{1+k_{1}S_{0}}\Psi\left(\dfrac{S(t)}{S_{0}}\right)+I(t)+\dfrac{1}{S_{0}(1+k_{1}S_{0})}\left(\dfrac{\lambda}{r}\dfrac{R^{2}(t)}{2}+\dfrac{\lambda}{4\mu}\dfrac{(N(t)-S_{0})^{2}}{2}\right),

where $\Psi(x)=x-1-\ln(x)$ , $x>0$ . Calculating the derivative of $L_{0}$ along the solution of system (1), and by using Lemma 1 in [2] and Lemma 3.1 in [51], we obtain that

	$\displaystyle D^{\alpha}L_{0}(t)$	$\displaystyle\leq\dfrac{1}{1+k_{1}S_{0}}\left(1-\dfrac{S_{0}}{S(t)}\right)D^{\alpha}S(t)+D^{\alpha}I(t)+\dfrac{\lambda}{rS_{0}(1+k_{1}S_{0})}R(t)D^{\alpha}R(t)$
		$\displaystyle\quad+\dfrac{\lambda}{4\mu S_{0}(1+k_{1}S_{0})}(N(t)-S_{0})D^{\alpha}N(t)$
		$\displaystyle\leq\dfrac{1}{1+k_{1}S_{0}}\left(1-\dfrac{S_{0}}{S(t)}\right)\mu(S_{0}-S(t))-\dfrac{1}{1+k_{1}S_{0}}\left(1-\dfrac{S_{0}}{S(t)}\right)\dfrac{\beta S(t)I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}$
		$\displaystyle\quad+\dfrac{\lambda}{1+k_{1}S_{0}}\left(1-\dfrac{S_{0}}{S(t)}\right)R(t)+\dfrac{\beta S(t)I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-aI(t)$
		$\displaystyle\quad+\dfrac{\lambda}{rS_{0}(1+k_{1}S_{0})}R(t)\left(rI(t)-(\mu+\lambda)R(t)\right)+\dfrac{\lambda}{4\mu S_{0}(1+k_{1}S_{0})}(N(t)-S_{0})(\Lambda-\mu N(t))$
		$\displaystyle\leq\dfrac{-\mu}{(1+k_{1}S_{0})S(t)}(S(t)-S_{0})^{2}+a(R_{0}-1)I(t)+\dfrac{\lambda}{S_{0}(1+k_{1}S_{0})}R(t)(N(t)-S_{0})$
		$\displaystyle\quad+\dfrac{\lambda}{1+k_{1}S_{0}}R(t)\left(\dfrac{S_{0}-S(t)}{S_{0}}+\dfrac{S(t)-S_{0}}{S(t)}\right)-\dfrac{\lambda}{S_{0}(1+k_{1}S_{0})}R^{2}(t)$
		$\displaystyle\quad-\dfrac{\lambda}{S_{0}(1+k_{1}S_{0})r}(\mu+\lambda)R^{2}(t)-\dfrac{\lambda}{4S_{0}(1+k_{1}S_{0})}(N(t)-S_{0})^{2}$
		$\displaystyle\leq\dfrac{-\mu}{(1+k_{1}S_{0})S(t)}(S(t)-S_{0})^{2}+a(R_{0}-1)I(t)-\dfrac{\lambda}{(1+k_{1}S_{0})}R(t)\dfrac{(S_{0}-S(t))^{2}}{S_{0}S(t)}$
		$\displaystyle\quad-\dfrac{\lambda}{S_{0}(1+k_{1}S_{0})}\left(R(t)-\dfrac{N(t)-S_{0}}{2}\right)^{2}-\dfrac{\lambda}{S_{0}(1+k_{1}S_{0})r}(\mu+\lambda)R^{2}(t).$

Therefore, $D^{\alpha}L_{0}(t)\leq 0$ if $R_{0}\leq 1$ . Furthermore, it is not hard to verify that the largest compact invariant set of $\left\{(S,I,R)\in\mathbb{R}_{+}^{3}:D^{\alpha}L_{0}(t)=0\right\}$ is the singleton $\{E_{0}\}$ . Therefore, from LaSalle invariance principle [26], we deduce that ${E_{0}}$ is globally asymptotically stable. ∎

Finally, we investigate the global stability of the endemic equilibrium ${E^{*}}$ .

Theorem 12.

The endemic equilibrium ${E^{*}}$ is globally asymptotically stable if $R_{0}>1$ and

R^{*}\leq\dfrac{\mu}{\lambda}S^{*}.

(8)

Proof.

To study the global stability of ${E^{*}}$ for (1), we propose the following Lyapunov function:

\begin{split}L^{*}(t)&=\dfrac{1+k_{2}S^{*}}{1+k_{1}S^{*}+k_{2}I^{*}+k_{3}S^{*}I^{*}}S^{*^{2}}\Psi\left(\dfrac{S(t)}{S^{*}}\right)+S^{*^{2}}I^{*}\Psi\left(\dfrac{I(t)}{I^{*}}\right)\\ &\quad+\dfrac{\lambda(1+k_{2}S^{*})}{4\mu(1+k_{1}S^{*}+k_{2}I^{*}+k_{3}S^{*}I^{*})}\Big{[}(S(t)-S^{*})+(I(t)-I^{*})+(R(t)-R^{*})\Big{]}^{2}\\ &\quad+\dfrac{\lambda(1+k_{2}S^{*})}{2r(1+k_{1}S^{*}+k_{2}I^{*}+k_{3}S^{*}I^{*})}(R(t)-R^{*})^{2},\end{split}

where $\Psi(x)=x-1-\ln(x)$ , $x>0$ . It follows that

\begin{split}D^{\alpha}L^{*}(t)&\leq\dfrac{1+k_{2}S^{*}}{1+k_{1}S^{*}+k_{2}I^{*}+k_{3}S^{*}I^{*}}S^{*}\left(1-\dfrac{S^{*}}{S(t)}\right)D^{\alpha}S(t)+S^{*}\left(1-\dfrac{I^{*}}{I(t)}\right)D^{\alpha}I(t)\\ &\quad+\dfrac{\lambda(1+k_{2}S^{*})}{2\mu(1+k_{1}S^{*}+k_{2}I^{*}+k_{3}S^{*}I^{*})}\Big{[}(S(t)-S^{*})+(I(t)-I^{*})+(R(t)-R^{*})\Big{]}D^{\alpha}N(t)\\ &\quad+\dfrac{\lambda(1+k_{2}S^{*})}{r(1+k_{1}S^{*}+k_{2}I^{*}+k_{3}S^{*}I^{*})}(R(t)-R^{*})D^{\alpha}R(t).\end{split}

Note that $\Lambda=\mu S^{*}+aI^{*}-\lambda R^{*}=\mu(S^{*}+I^{*}+R^{*})$ ,

\dfrac{\beta S^{*}}{1+k_{1}S^{*}+k_{2}I^{*}+k_{3}S^{*}I^{*}}=a,

and $(\lambda+\mu)R^{*}-rI^{*}=0$ . Then,

	$\displaystyle D^{\alpha}L^{*}(t)$	$\displaystyle\leq\dfrac{1+k_{2}S^{}}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}S^{}\left(1-\dfrac{S^{*}}{S(t)}\right)\Bigg{(}\Lambda-\mu S(t)-\dfrac{\beta S(t)I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}$
		$\displaystyle\quad+\lambda R(t)\Bigg{)}+S^{}\left(1-\dfrac{I^{}}{I(t)}\right)\left(\dfrac{\beta S(t)I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-aI(t)\right)$
		$\displaystyle\quad+\dfrac{\lambda(1+k_{2}S^{})}{2\mu(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})}\Big{[}(S(t)-S^{})+(I(t)-I^{})+(R(t)-R^{})\Big{]}\Big{(}\Lambda-\mu N(t)\Big{)}$
		$\displaystyle\quad+\dfrac{\lambda(1+k_{2}S^{})}{r(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})}\left(R(t)-R^{}\right)\left(rI(t)-(\mu+\lambda)R(t)\right)$
		$\displaystyle\leq\dfrac{\mu S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{}}{S(t)}(S^{}-S(t))+aI^{}\dfrac{S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{}}{S(t)}$
		$\displaystyle\quad-\dfrac{S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{*}}{S(t)}\dfrac{\beta S(t)I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}$
		$\displaystyle\quad+\dfrac{\lambda S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{}}{S(t)}(R(t)-R^{})$
		$\displaystyle\quad+S^{}\left(1-\dfrac{I^{}}{I(t)}\right)\dfrac{\beta S(t)I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-aS^{}(I(t)-I^{})$
		$\displaystyle\quad+\dfrac{\lambda(1+k_{2}S^{})}{2\mu(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})}\Big{[}(S(t)-S^{})+(I(t)-I^{*})$
		$\displaystyle\quad+(R(t)-R^{})\Big{]}\Big{[}\mu(S^{}-S(t))+\mu(I^{}-I(t))+\mu(R^{}-R(t))\Big{]}$
		$\displaystyle\quad+\dfrac{\lambda(1+k_{2}S^{})}{r(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})}\left(R(t)-R^{}\right)\left(r(I(t)-I^{})-(\mu+\lambda)(R(t)-R^{})\right)$
		$\displaystyle\leq-\dfrac{\mu S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}+aI^{}\dfrac{S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{*}}{S(t)}$
		$\displaystyle\quad-\dfrac{a(1+k_{2}S^{})(S(t)-S^{})I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}+\left(1-\dfrac{I^{}}{I(t)}\right)\dfrac{a(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{*})S(t)I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}$
		$\displaystyle\quad-aS^{}I(t)+aS^{}I^{}+\dfrac{\lambda S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{}}{S(t)}(R(t)-R^{*})$
		$\displaystyle\quad-\dfrac{\lambda(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}(S(t)-S^{})(R(t)-R^{*})$
		$\displaystyle\leq-\dfrac{\mu S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}+aI^{}\dfrac{S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{*}}{S(t)}$
		$\displaystyle\quad+\dfrac{aS^{}(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{})I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-\dfrac{aI^{}(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})S(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-aS^{*}I(t)$
		$\displaystyle\quad+aS^{}I^{}+\dfrac{\lambda(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}(S(t)-S^{})(R(t)-R^{})\dfrac{S^{}-S(t)}{S(t)}$
		$\displaystyle\leq-\dfrac{\mu S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}+aI^{}\dfrac{S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{*}}{S(t)}$
		$\displaystyle\quad+\dfrac{aS^{}(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{})I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-\dfrac{aI^{}(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})S(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-aS^{*}I(t)$
		$\displaystyle\quad+aS^{}I^{}-\dfrac{\lambda(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}(R(t)-R^{*}))$
		$\displaystyle\leq-\dfrac{\mu S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{*})^{2}}{S(t)}$
		$\displaystyle\quad+aS^{}I^{}\left(\dfrac{(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{}}{S(t)}+\dfrac{(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{})I(t)}{(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t))I^{*}}\right)$
		$\displaystyle\quad+aS^{}I^{}\left(1-\frac{I(t)}{I^{}}-\dfrac{(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})S(t)}{(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t))S^{}}\right)$
		$\displaystyle\quad+R^{}\dfrac{\lambda(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{*})^{2}}{S(t)}$
		$\displaystyle\leq R^{}\dfrac{\lambda(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}-\dfrac{\mu S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}$
		$\displaystyle\quad+aS^{}I^{}\Bigg{(}\dfrac{(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{}}{S(t)}-1+1-\dfrac{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}{1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}}$
		$\displaystyle\quad-\dfrac{(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})S(t)}{(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t))S^{}}\Bigg{)}+aS^{}I^{}\Bigg{(}2-1-\frac{I(t)}{I^{}}$
		$\displaystyle\quad+\dfrac{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}{1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}}+\dfrac{(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{})I(t)}{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)I^{*}}\Bigg{)}$
		$\displaystyle\leq(\lambda R^{}-\mu S^{})\dfrac{(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}$
		$\displaystyle\quad+aS^{}I^{}\Bigg{(}3-\dfrac{(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{})S^{}}{(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})S(t)}-\dfrac{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}{1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{*}}$
		$\displaystyle\quad-\dfrac{(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})S(t)}{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)S^{}}\Bigg{)}+aS^{}I^{}\Bigg{(}-1-\frac{I(t)}{I^{}}$
		$\displaystyle\quad+\dfrac{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}{1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}}+\dfrac{(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{})I(t)}{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)I^{*}}\Bigg{)}$
		$\displaystyle\leq\left(\lambda R^{}-\mu S^{}\right)\dfrac{\left(1+k_{2}S^{}\right)}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}$
		$\displaystyle\quad+aS^{}I^{}\Bigg{[}\Bigg{(}3-\dfrac{\left(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}\right)S^{}}{\left(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}\right)S(t)}-\dfrac{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}{1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{*}}$
		$\displaystyle\quad-\dfrac{(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})S(t)}{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)S^{}}\Bigg{)}+\ln\dfrac{\left(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}\right)S^{}}{\left(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}\right)S(t)}$
		$\displaystyle\quad+\ln\dfrac{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)S^{}}{\left(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}\right)S^{}}+\ln\dfrac{\left(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}\right)S(t)}{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)S^{*}}\Bigg{]}$
		$\displaystyle\quad+aS^{}I^{}\left[\dfrac{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}{1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}}-1\right]$
		$\displaystyle\quad+aS^{}I(t)\left[\dfrac{1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{*}}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-1\right]$
		$\displaystyle\leq(\lambda R^{}-\mu S^{})\dfrac{(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}$
		$\displaystyle\quad-aS^{}I^{}\Bigg{[}\Psi\left(\dfrac{\left(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}\right)S^{}}{\left(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}\right)S(t)}\right)+\Psi\left(\dfrac{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)S^{}}{\left(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}\right)S^{*}}\right)$
		$\displaystyle\quad+\Psi\left(\dfrac{\left(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}\right)S(t)}{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)S^{*}}\right)\Bigg{]}$
		$\displaystyle\quad-\dfrac{aS^{}\left(k_{2}+k_{3}S(t)\right)\left(1+k_{1}S(t)\right)\left(I(t)-I^{}\right)^{2}}{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)\left(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}\right)}.$

Clearly, $\Psi(x)\geq 0$ . Consequently, $D^{\alpha}L^{*}(t)\leq 0$ if $R_{0}>0$ and $R^{*}\leq\frac{\mu}{\lambda}S^{*}$ . Further, the largest invariant set of $\{(S,I,R)\in\mathbb{R}_{+}^{3}:D^{\alpha}L^{*}(t)=0\}$ is the singleton $\{E^{*}\}$ . Therefore, from LaSalle’s invariance principle, ${E^{*}}$ is globally asymptotically stable. ∎

It is easy to see that condition (8) in Theorem 12 is equivalent to

\displaystyle\begin{aligned} R_{0}\leq 1&+\dfrac{k_{3}a\beta\Lambda^{2}(\lambda^{2}r+a(\mu+r))(\mu+r)(\mu+\lambda)+k_{1}a^{2}\beta\lambda^{2}r(\mu+\lambda)^{2}}{\lambda^{2}k_{3}a(\mu+r)^{2}(\mu+\lambda)(\mu+k_{1}\Lambda)}\\ &+\dfrac{k_{1}a\beta(a^{2}(\mu+\lambda)^{2}+\lambda^{2}r^{2})(\mu+\lambda)+\beta(\mu+r)\sqrt{\Delta}}{\lambda^{2}k_{3}a(\mu+r)^{2}(\mu+\lambda)(\mu+k_{1}\Lambda)}\end{aligned}

(9)

and

\displaystyle\begin{aligned} &\lim\limits_{\lambda\rightarrow 0}\Bigg{(}\dfrac{k_{3}a\beta\Lambda^{2}(\lambda^{2}r+a(\mu+r))(\mu+r)(\mu+\lambda)+k_{1}a^{2}\beta\lambda^{2}r(\mu+\lambda)^{2}}{\lambda^{2}k_{3}a(\mu+r)^{2}(\mu+\lambda)(\mu+k_{1}\Lambda)}\\ &\qquad\quad+\dfrac{k_{1}a\beta(a^{2}(\mu+\lambda)^{2}+\lambda^{2}r^{2})(\mu+\lambda)+\beta(\mu+r)\sqrt{\Delta}}{\lambda^{2}k_{3}a(\mu+r)^{2}(\mu+\lambda)(\mu+k_{1}\Lambda)}\Bigg{)}=+\infty.\end{aligned}

(10)

Corollary 13.

The endemic equilibrium ${E^{*}}$ is globally asymptotically stable if $R_{0}>1$ and $\lambda$ is sufficiently small.

6. Numerical simulations

In this section we present some numerical simulations in order to illustrate our theoretical results. The system (1) is numerically integrated by using the fractional Euler’s method [37]. The approximate solutions of (2) with $0<\alpha\leq 1$ are displayed in Figures 1 and 2. The solutions converge to the equilibrium points. The parameter values used in the simulations are: $\Lambda=0.8$ , $\mu=0.1$ , $\lambda=0.5$ , $\beta=0.1$ , $r=0.5$ , $k_{1}=0.1$ , $k_{2}=0.02$ , and $k_{3}=0.003$ with initial conditions $S(0)=10.0$ , $I(0)=1.0$ , $R(0)=1.0$ . Using the MATLAB numerical computing environment, we get $R_{0}=0.7407$ . Hence, system (1) has a unique disease-free equilibrium $E_{0}=(8,0,0)$ . According to Theorem 11, $E_{0}$ is globally asymptotically stable (see Figure 1).

Refer to caption — Figure 1. Stability of the disease-free equilibrium $E_{0}$ .

Now, let us choose $\mu=0.02$ and keep the other parameter values. Then, $R_{0}=1.5385$ and $R^{*}=0.552\leq 21.8944$ . Hence, the condition (8) is satisfied, as well as conditions (9) and (10), and the model converges rapidly to its steady state for $\alpha=0.85$ , when compared to other fractional derivatives (see Figure 2).

7. Conclusion

The use of fractional order derivatives can help to reduce errors arising from the neglected parameters in modelling real life phenomena [18]. Here, we have studied a fractional-order SIRS epidemic model with a general incidence function. The stability of equilibrium points is investigated and numerical solutions are given. According to our theoretical analysis, the fractional order parameter $\alpha$ has no effect on the stability of free and endemic equilibria, but it can affect the time for arriving at the steady states. As future work, we plan to study the stability of a more general SIRS type model taking into account other parameters.

Acknowledgements

This work was partially supported by Fundação para a Ciência e a Tecnologia (FCT) within the R&D unit Centro de Investigação e Desenvolvimento em Matemática e Aplicações (CIDMA), project UIDB/04106/2020. The authors are very grateful to two anonymous reviewers, for several critical remarks and precious suggestions.

References

[1] R. P. Agarwal, D. Baleanu, J. J. Nieto, D. F. M. Torres and Y. Zhou, A survey on fuzzy fractional differential and optimal control nonlocal evolution equations, J. Comput. Appl. Math. 339 (2018), 3–29. arXiv:1709.07766
[2] N. Aguila-Camacho, M. A. Duarte-Mermoud and J. A. Gallegos, Lyapunov functions for fractional order systems, Commun. Nonlinear Sci. Numer. Simul. 19 (2014), no. 9, 2951–2957.
[3] R. Almeida, What is the best fractional derivative to fit data?, Appl. Anal. Discrete Math. 11 (2017), no. 2, 358–368.
[4] R. Almeida, Analysis of a fractional SEIR model with treatment, Appl. Math. Lett. 84 (2018), 56–62.
[5] R. Almeida, A. M. C. Brito da Cruz, N. Martins and M. T. T. Monteiro, An epidemiological MSEIR model described by the Caputo fractional derivative, Int. J. Dyn. Control 7 (2019), no. 2, 776–784.
[6] R. Almeida, S. Pooseh and D. F. M. Torres, Computational methods in the fractional calculus of variations, Imperial College Press, London, 2015.
[7] R. Almeida, D. Tavares and D. F. M. Torres, The variable-order fractional calculus of variations, SpringerBriefs in Applied Sciences and Technology, Springer, Cham, 2019. arXiv:1805.00720
[8] B. Bayour and D. F. M. Torres, Complex-valued fractional derivatives on time scales, in Differential and difference equations with applications, 79–87, Springer Proc. Math. Stat., 164, Springer, Cham, 2016. arXiv:1511.02153
[9] J. R. Beddington, Mutual interference between parasites or predators and its effect on searching efficiency, J. Anim. Ecol. 44 (1975), 331–340.
[10] A. R. M. Carvalho and C. M. A. Pinto, Non-integer order analysis of the impact of diabetes and resistant strains in a model for TB infection, Commun. Nonlinear Sci. Numer. Simul. 61 (2018), 104–126.
[11] A. R. M. Carvalho and C. M. A. Pinto, Immune response in HIV epidemics for distinct transmission rates and for saturated CTL response, Math. Model. Nat. Phenom. 14 (2019), no. 3, Art. 307, 13 pp.
[12] A. R. M. Carvalho, C. M. A. Pinto and D. Baleanu, HIV/HCV coinfection model: a fractional-order perspective for the effect of the HIV viral load, Adv. Difference Equ. 2018 (2018), Paper no. 2, 22 pp.
[13] P. H. Crowley and E. K. Martin, Functional responses and interference within and between year classes of a dragonfly population, J. North Am. Benthol. Soc. 8 (1989), 211–221.
[14] D. L. DeAngelis, R. A. Goldsten and R. V. O’Neill, A model for trophic interaction, Ecology 56 (1975), 881–892.
[15] A. Debbouche, J. J. Nieto and D. F. M. Torres, Optimal solutions to relaxation in multiple control problems of Sobolev type with nonlocal nonlinear fractional differential equations, J. Optim. Theory Appl. 174 (2017), no. 1, 7–31. arXiv:1504.05153
[16] K. Diethelm, N. J. Ford and A. D. Freed, Detailed error analysis for a fractional Adams method, Numer. Algorithms 36 (2004), no. 1, 31–52.
[17] E. F. Doungmo Goufo, R. Maritz and J. Munganga, Some properties of the Kermack-McKendrick epidemic model with fractional derivative and nonlinear incidence, Adv. Difference Equ. 2014 (2014), no. 278, 9 pp.
[18] H. A. A. El-Saka, The fractional-order SIR and SIRS epidemic models with variable population size, Math. Sci. Lett. 2 (2013), 195–200.
[19] M. El-Shahed and A. Alsaedi, The fractional SIRC model and influenza A, Math. Probl. Eng. 2011 (2011), Art. ID 480378, 9 pp.
[20] K. Hattaf, N. Yousfi and A. Tridane, Stability analysis of a virus dynamics model with general incidence rate and two delays, Appl. Math. Comput. 221 (2013), 514–521.
[21] J.-H. He, Variational iteration method for autonomous ordinary differential systems, Appl. Math. Comput. 114 (2000), 115–123.
[22] J.-H. He, S. K. Elagan and Z. B. Li, Geometrical explanation of the fractional complex transform and derivative chain rule for fractional calculus, Phys. Lett. A 376 (2012), no. 4, 257–259.
[23] J.-H. He and Y. Hu, On fractals Space-time and fractional calculus, Thermal Science 20 (2016), no. 3, 773–777.
[24] J.-H. He and Z. B. Li, Converting fractional differential equations into partial differential equations, Thermal Science 16 (2012), no. 2, 331–334.
[25] R. Hilfer, Applications of fractional calculus in physics, World Scientific Publishing Co., Inc., River Edge, NJ, 2000.
[26] J. Huo, H. Zhao and L. Zhu, The effect of vaccines on backward bifurcation in a fractional order HIV model, Nonlinear Anal. Real World Appl. 26 (2015), 289–305.
[27] W. O. Kermack and A. G. McKendrick, A contribution to the mathematical theory of epidemics, part I, Proc. R. Soc. A 115 (1927), 700–721.
[28] H. Khosravian-Arab and D. F. M. Torres, Uniform approximation of fractional derivatives and integrals with application to fractional differential equations, Nonlinear Stud. 20 (2013), 533–548. arXiv:1308.0451
[29] Z.-B. Li and J.-H. He, Fractional complex transform for fractional differential equations, Math. Comput. Appl. 15 (2010), no. 5, 970–973.
[30] W. Lin, Global existence theory and chaos control of fractional differential equations, J. Math. Anal. Appl. 332 (2007), no. 1, 709–726.
[31] H.-Y. Liu, J.-H. He and Z.-B. Li, Fractional calculus for nanoscale flow and heat transfer, Internat. J. Numer. Methods Heat Fluid Flow 24 (2014), no. 6, 1227–1250.
[32] E. M. Lotfi, M. Mahrouf, M. Maziane, C. J. Silva, D. F. M. Torres and N. Yousfi, A minimal HIV-AIDS infection model with general incidence rate and application to Morocco data, Stat. Optim. Inf. Comput. 7 (2019), no. 3, 588–603. arXiv:1812.06965
[33] A. B. Malinowska, T. Odzijewicz and D. F. M. Torres, Advanced methods in the fractional calculus of variations, SpringerBriefs in Applied Sciences and Technology, Springer, Cham, 2015.
[34] J. P. Mateus, P. Rebelo, S. Rosa, C. M. Silva and D. F. M. Torres, Optimal control of non-autonomous SEIRS models with vaccination and treatment, Discrete Contin. Dyn. Syst. Ser. S 11 (2018), no. 6, 1179–1199. arXiv:1706.06843
[35] Z. M. Odibat and S. Momani, Application of variational iteration method to nonlinear differential equations of fractional order, Int. J. Nonlinear Sci. Numer. Simul. 7 (2006), no. 1, 27–34.
[36] Z. M. Odibat and S. Momani, Approximate solutions for boundary value problems of time-fractional wave equation, Appl. Math. Comput. 181 (2006), no. 1, 767–774.
[37] Z. Odibat and S. Momani, An algorithm for the numerical solution of differential equations of fractional order, J. Appl. Math. Inform. 26 (2008), 15–27.
[38] Z. M. Odibat and N. T. Shawagfeh, Generalized Taylor’s formula, Appl. Math. Comput. 186 (2007), no. 1, 286–293.
[39] E. Okyere, F. T. Oduro, S. K. Amponsah, I. K. Dontwi and N. K. Frempong, Fractional order SIR model with constant population, British J. Math. Comput. Sci. 14 (2016), no. 2, 1–12.
[40] N. Özalp and E. Demirci, A fractional order SEIR model with vertical transmission, Math. Comput. Modelling 54 (2011), no. 1-2, 1–6.
[41] I. Petráš, Fractional-order nonlinear systems: modeling, analysis and simulation, Springer, Berlin, 2011.
[42] I. Podlubny, Fractional differential equations, Mathematics in Science and Engineering, 198, Academic Press, Inc., San Diego, CA, 1999.
[43] A. Rachah and D. F. M. Torres, Analysis, simulation and optimal control of a SEIR model for Ebola virus with demographic effects, Commun. Fac. Sci. Univ. Ank. Sér. A1 Math. Stat. 67 (2018), no. 1, 179–197. arXiv:1705.01079
[44] K. Razminia, A. Razminia and D. F. M. Torres, Pressure responses of a vertically hydraulic fractured well in a reservoir with fractal structure, Appl. Math. Comput. 257 (2015), 374–380. arXiv:1502.03944
[45] S. Rosa and D. F. M. Torres, Optimal control of a fractional order epidemic model with application to human respiratory syncytial virus infection, Chaos Solitons Fractals 117 (2018), 142–149. arXiv:1810.06900
[46] D. Rostamy and E. Mottaghi, Numerical solution and stability analysis of a nonlinear vaccination model with historical effects, J. Math. Stat. (2017), Conf. Proc., 17 pp. DOI: 10.15672/HJMS.20174720333
[47] A. B. Salati, M. Shamsi and D. F. M. Torres, Direct transcription methods based on fractional integral approximation formulas for solving nonlinear fractional optimal control problems, Commun. Nonlinear Sci. Numer. Simul. 67 (2019), 334–350. arXiv:1805.06537
[48] S. G. Samko, A. A. Kilbas and O. I. Marichev, Fractional integrals and derivatives, Gordon and Breach Science Publishers, Yverdon, 1993.
[49] D. Sierociuk, T. Skovranek, M. Macias, I. Podlubny, I. Petras, A. Dzielinski and P. Ziubinski, Diffusion process modeling by using fractional-order models, Appl. Math. Comput. 257 (2015), 2–11.
[50] A. A. Stanislavsky, Fractional oscillator, Phys. Rev. E 70 (2004), 051103, 6 pp.
[51] C. Vargas-De-León, Volterra-type Lyapunov functions for fractional-order epidemic systems, Commun. Nonlinear Sci. Numer. Simul. 24 (2015), no. 1-3, 75–85.
[52] X. Wang, Y. He and M. Wang, Chaos control of a fractional order modified coupled dynamos system, Nonlinear Anal. 71 (2009), no. 12, 6126–6134.
[53] K.-L. Wang and S. Y. Liu, He’s fractional derivative and its application for fractional Fornberg-Whitham equation, Thermal Science 21 (2017), no. 5, 2049–2055.
[54] B. J. West, M. Turalska and P. Grigolini, Fractional calculus ties the microscopic and macroscopic scales of complex network dynamics, New J. Phys. 17 (2015), no. 4, 045009, 13 pp.
[55] W. Wojtak, C. J. Silva and D. F. M. Torres, Uniform asymptotic stability of a fractional tuberculosis model, Math. Model. Nat. Phenom. 13 (2018), no. 1, Art. 9, 10 pp. arXiv:1801.07059
[56] X.-E. Wu and Y.-S. Liang, Relationship between fractal dimensions and fractional calculus, Nonlinear Sci. Lett. A 8 (2017), no. 1, 77–89.
[57] X. Wu, L. Tang and T. Zhong, An iteration algorithm for fractal dimensions of a self-similar set, Nonlinear Sci. Lett. A 8 (2017), no. 1, 117–120.
[58] X.-J. Yang, J. A. T. Machado and D. Baleanu, Exact traveling-wave solution for local fractional Boussinesq equation in fractal domain, Fractals 25 (2017), no. 4, 1740006, 7 pp.

	$\displaystyle D^{\alpha}L^{*}(t)$	$\displaystyle\leq\dfrac{1+k_{2}S^{}}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}S^{}\left(1-\dfrac{S^{*}}{S(t)}\right)\Bigg{(}\Lambda-\mu S(t)-\dfrac{\beta S(t)I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}$
		$\displaystyle\quad+\lambda R(t)\Bigg{)}+S^{}\left(1-\dfrac{I^{}}{I(t)}\right)\left(\dfrac{\beta S(t)I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-aI(t)\right)$
		$\displaystyle\quad+\dfrac{\lambda(1+k_{2}S^{})}{2\mu(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})}\Big{[}(S(t)-S^{})+(I(t)-I^{})+(R(t)-R^{})\Big{]}\Big{(}\Lambda-\mu N(t)\Big{)}$
		$\displaystyle\quad+\dfrac{\lambda(1+k_{2}S^{})}{r(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})}\left(R(t)-R^{}\right)\left(rI(t)-(\mu+\lambda)R(t)\right)$
		$\displaystyle\leq\dfrac{\mu S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{}}{S(t)}(S^{}-S(t))+aI^{}\dfrac{S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{}}{S(t)}$
		$\displaystyle\quad-\dfrac{S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{*}}{S(t)}\dfrac{\beta S(t)I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}$
		$\displaystyle\quad+\dfrac{\lambda S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{}}{S(t)}(R(t)-R^{})$
		$\displaystyle\quad+S^{}\left(1-\dfrac{I^{}}{I(t)}\right)\dfrac{\beta S(t)I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-aS^{}(I(t)-I^{})$
		$\displaystyle\quad+\dfrac{\lambda(1+k_{2}S^{})}{2\mu(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})}\Big{[}(S(t)-S^{})+(I(t)-I^{*})$
		$\displaystyle\quad+(R(t)-R^{})\Big{]}\Big{[}\mu(S^{}-S(t))+\mu(I^{}-I(t))+\mu(R^{}-R(t))\Big{]}$
		$\displaystyle\quad+\dfrac{\lambda(1+k_{2}S^{})}{r(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})}\left(R(t)-R^{}\right)\left(r(I(t)-I^{})-(\mu+\lambda)(R(t)-R^{})\right)$
		$\displaystyle\leq-\dfrac{\mu S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}+aI^{}\dfrac{S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{*}}{S(t)}$
		$\displaystyle\quad-\dfrac{a(1+k_{2}S^{})(S(t)-S^{})I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}+\left(1-\dfrac{I^{}}{I(t)}\right)\dfrac{a(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{*})S(t)I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}$
		$\displaystyle\quad-aS^{}I(t)+aS^{}I^{}+\dfrac{\lambda S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{}}{S(t)}(R(t)-R^{*})$
		$\displaystyle\quad-\dfrac{\lambda(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}(S(t)-S^{})(R(t)-R^{*})$
		$\displaystyle\leq-\dfrac{\mu S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}+aI^{}\dfrac{S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{*}}{S(t)}$
		$\displaystyle\quad+\dfrac{aS^{}(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{})I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-\dfrac{aI^{}(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})S(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-aS^{*}I(t)$
		$\displaystyle\quad+aS^{}I^{}+\dfrac{\lambda(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}(S(t)-S^{})(R(t)-R^{})\dfrac{S^{}-S(t)}{S(t)}$
		$\displaystyle\leq-\dfrac{\mu S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}+aI^{}\dfrac{S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{*}}{S(t)}$
		$\displaystyle\quad+\dfrac{aS^{}(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{})I(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-\dfrac{aI^{}(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})S(t)}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-aS^{*}I(t)$
		$\displaystyle\quad+aS^{}I^{}-\dfrac{\lambda(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}(R(t)-R^{*}))$
		$\displaystyle\leq-\dfrac{\mu S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{*})^{2}}{S(t)}$
		$\displaystyle\quad+aS^{}I^{}\left(\dfrac{(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{}}{S(t)}+\dfrac{(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{})I(t)}{(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t))I^{*}}\right)$
		$\displaystyle\quad+aS^{}I^{}\left(1-\frac{I(t)}{I^{}}-\dfrac{(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})S(t)}{(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t))S^{}}\right)$
		$\displaystyle\quad+R^{}\dfrac{\lambda(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{*})^{2}}{S(t)}$
		$\displaystyle\leq R^{}\dfrac{\lambda(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}-\dfrac{\mu S^{}(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}$
		$\displaystyle\quad+aS^{}I^{}\Bigg{(}\dfrac{(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{S(t)-S^{}}{S(t)}-1+1-\dfrac{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}{1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}}$
		$\displaystyle\quad-\dfrac{(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})S(t)}{(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t))S^{}}\Bigg{)}+aS^{}I^{}\Bigg{(}2-1-\frac{I(t)}{I^{}}$
		$\displaystyle\quad+\dfrac{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}{1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}}+\dfrac{(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{})I(t)}{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)I^{*}}\Bigg{)}$
		$\displaystyle\leq(\lambda R^{}-\mu S^{})\dfrac{(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}$
		$\displaystyle\quad+aS^{}I^{}\Bigg{(}3-\dfrac{(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{})S^{}}{(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})S(t)}-\dfrac{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}{1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{*}}$
		$\displaystyle\quad-\dfrac{(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})S(t)}{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)S^{}}\Bigg{)}+aS^{}I^{}\Bigg{(}-1-\frac{I(t)}{I^{}}$
		$\displaystyle\quad+\dfrac{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}{1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}}+\dfrac{(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{})I(t)}{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)I^{*}}\Bigg{)}$
		$\displaystyle\leq\left(\lambda R^{}-\mu S^{}\right)\dfrac{\left(1+k_{2}S^{}\right)}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}$
		$\displaystyle\quad+aS^{}I^{}\Bigg{[}\Bigg{(}3-\dfrac{\left(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}\right)S^{}}{\left(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}\right)S(t)}-\dfrac{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}{1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{*}}$
		$\displaystyle\quad-\dfrac{(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{})S(t)}{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)S^{}}\Bigg{)}+\ln\dfrac{\left(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}\right)S^{}}{\left(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}\right)S(t)}$
		$\displaystyle\quad+\ln\dfrac{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)S^{}}{\left(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}\right)S^{}}+\ln\dfrac{\left(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}\right)S(t)}{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)S^{*}}\Bigg{]}$
		$\displaystyle\quad+aS^{}I^{}\left[\dfrac{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}{1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}}-1\right]$
		$\displaystyle\quad+aS^{}I(t)\left[\dfrac{1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{*}}{1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)}-1\right]$
		$\displaystyle\leq(\lambda R^{}-\mu S^{})\dfrac{(1+k_{2}S^{})}{1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}}\dfrac{(S(t)-S^{})^{2}}{S(t)}$
		$\displaystyle\quad-aS^{}I^{}\Bigg{[}\Psi\left(\dfrac{\left(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}\right)S^{}}{\left(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}\right)S(t)}\right)+\Psi\left(\dfrac{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)S^{}}{\left(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}\right)S^{*}}\right)$
		$\displaystyle\quad+\Psi\left(\dfrac{\left(1+k_{1}S^{}+k_{2}I^{}+k_{3}S^{}I^{}\right)S(t)}{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)S^{*}}\right)\Bigg{]}$
		$\displaystyle\quad-\dfrac{aS^{}\left(k_{2}+k_{3}S(t)\right)\left(1+k_{1}S(t)\right)\left(I(t)-I^{}\right)^{2}}{\left(1+k_{1}S(t)+k_{2}I(t)+k_{3}S(t)I(t)\right)\left(1+k_{1}S(t)+k_{2}I^{}+k_{3}S(t)I^{}\right)}.$

Global stability of a Caputo fractional SIRS model with general incidence rate

Abstract

keywords:

1991 Mathematics Subject Classification:

1. Introduction

2. Basic results of fractional calculus

Definition 1 (See [42]).

Definition 2 (See [42]).

Definition 3 (See [42]).

Lemma 4 (See [30]).

3. Existence and uniqueness of positive solution

Lemma 5 (Generalized Mean Value Theorem [38]).

Corollary 6 (See [38]).

Theorem 7.

Proof.

4. Local stability

Theorem 8.

Proof.

Theorem 9.

Proof.

Theorem 10.

Proof.

5. Global stability

Theorem 11.

Proof.

Theorem 12.

Proof.

Corollary 13.

6. Numerical simulations

7. Conclusion

Acknowledgements

References

Global stability of a Caputo fractional SIRS model
with general incidence rate