Stability of a fractional HIV/AIDS model

TL;DR

This paper introduces a fractional order model for HIV/AIDS transmission, analyzing its stability properties both theoretically and through numerical simulations.

Contribution

It presents a novel fractional order HIV/AIDS model and studies its stability, which is a new approach compared to traditional models.

Findings

The model's local stability is established.

Numerical simulations confirm the theoretical stability results.

Abstract

We propose a fractional order model for HIV/AIDS transmission. Local and uniform stability of the fractional order model is studied. The theoretical results are illustrated through numerical simulations.

Equations

$${}_a^C{D}_t^{\alpha }f(t)=\frac{1}{\Gamma (n-\alpha )}{\int }_a^t\frac{f^{(n)}(\xi )}{(t-\xi {)}^{\alpha +1-n}}d\xi ,n-1<\alpha <n\in \mathbb{N}.$$

$${}_a^C{D}_t^{\alpha }(c_{1}f(t)+c_{2}g(t))=c_1{}_a^C{D}_t^{\alpha }f(t)+c_2{}_a^C{D}_t^{\alpha }g(t).$$

$${}_a^C{D}_t^{\alpha }c=0.$$

$${}_a^C{D}_t^{\alpha }x(t)=f(t,x(t)),\alpha \in (0,1),$$

$$|\arg (\lambda )|>\frac{\alpha \pi }{2}.$$

$$W_1(x)\le L(t,x(t))\le W_2(x)$$

$${}_a^C{D}_t^{\alpha }L(t,x(t))\le -W_3(x)$$

$${}_{t_0}^C{D}_t^{\alpha }\left[x(t)-x^{\ast }-x^{\ast }\ln \frac{x(t)}{x^{\ast }}\right]\le \left(1-\frac{x^{\ast }}{x(t)}\right){{}_{t_0}^{C}D}_t^{\alpha }x(t),x^{\ast }\in {\mathbb{R}}^{+},\forall \alpha \in (0,1).$$

$$\lambda =\beta \left(I+{\eta }_{C}C+{\eta }_{A}A\right),$$

$$\{\begin{array}{l}{}_{t_0}^C{D}_t^{\alpha }S(t)=\Lambda -\beta \left(I(t)+{\eta }_{C}C(t)+{\eta }_{A}A(t)\right)S(t)-\mu S(t),\\ {}_{t_0}^C{D}_t^{\alpha }I(t)=\beta \left(I(t)+{\eta }_{C}C(t)+{\eta }_{A}A(t)\right)S(t)-\left(\rho +\varphi +\mu \right)I(t)+\omega C(t)+\gamma A(t),\\ {}_{t_0}^C{D}_t^{\alpha }C(t)=\varphi I(t)-(\omega +\mu )C(t),\\ {}_{t_0}^C{D}_t^{\alpha }A(t)=\rho I(t)-(\gamma +\mu +d)A(t).\end{array}$$

$$\Omega =\left\{\left(S,I,C,A\right)\in {\mathbb{R}}_{+}^4:N\le \frac{\Lambda }{\mu }\right\}.$$

$${\Sigma }_0=\left(S^0,I^0,C^0,A^0\right)=\left(\frac{\Lambda }{\mu },0,0,0\right).$$

$$R_0=\frac{S^0\beta \left({\xi }_2\left({\xi }_1+\rho {\eta }_A\right)+{\eta }_C\varphi {\xi }_1\right)}{\mu \left({\xi }_2\left(\rho +{\xi }_1\right)+\varphi {\xi }_1+\rho d\right)+\rho \omega d}=\frac{S^0\mathcal{N}}{\mathcal{D}},$$

$$\mathcal{N}=\beta \left({\xi }_2\left({\xi }_1+\rho {\eta }_A\right)+{\eta }_C\varphi {\xi }_1\right)$$

$$\mathcal{D}=\mu \left({\xi }_2\left(\rho +{\xi }_1\right)+\varphi {\xi }_1+\rho d\right)+\rho \omega d.$$

$$S^{\ast }=\frac{\mathcal{D}}{\mathcal{N}},I^{\ast }=\frac{{\xi }_1{\xi }_2(\Lambda \mathcal{N}-\mu \mathcal{D})}{\mathcal{DN}},C^{\ast }=\frac{\varphi {\xi }_1(\Lambda \mathcal{N}-\mu \mathcal{D})}{\mathcal{DN}},A^{\ast }=\frac{\rho {\xi }_2\left(\Lambda \mathcal{N}-\mu \mathcal{D}\right)}{\mathcal{DN}}.$$

$$J(\Sigma_{0})=\left[\begin{array}[]{cccc}-\mu&-{\frac{\beta\,\Lambda}{\mu}}&-{\frac{\beta\,\Lambda\,\eta_{C}}{\mu}}&-{\frac{\Lambda\,\beta\,\eta_{A}}{\mu}}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&{\frac{\Lambda\,\beta}{\mu}}-\mu-\phi-\rho&{\frac{\Lambda\,\beta\,\eta_{C}}{\mu}}+\omega&{\frac{\Lambda\,\beta\,\eta_{A}}{\mu}}+\alpha\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&\phi&-\omega-\mu&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&\rho&0&-\gamma-\mu-d\end{array}\right].$$

$$|\arg (\lambda )|>\frac{\alpha \pi }{2}.$$

$$qp=0$$

$$q=(\lambda +\mu )$$

$$p={\lambda }^3+b_1{\lambda }^2+b_2\lambda +b_3,$$

$$\begin{array}{rl}{b}_1& =-\frac{\Lambda \beta -\mu ({\xi }_1+{\xi }_2+{\xi }_3)}{\mu },\\ b_2& =-\frac{1}{\mu }\left(\Lambda \beta \left({\eta }_A\rho +{\eta }_C\varphi +{\xi }_1+{\xi }_2\right)-\mu \left(d({\xi }_2+{\xi }_3)+\gamma (\mu +{\xi }_2+\varphi )+\mu ({\xi }_2+\omega +2{\xi }_3)+\omega \rho \right)\right),\\ b_3& =-\frac{1}{\mu }\left(\Lambda \mathcal{N}-\mu \mathcal{D}\right).\end{array}$$

$$D(p)=-\left|\begin{array}[]{c c c c c}1&b_{1}&b_{2}&b_{3}&0\\ 0&1&b_{1}&b_{2}&b_{3}\\ 3&2b_{1}&b_{2}&0&0\\ 0&3&2b_{1}&b_{2}&0\\ 0&0&3&2b_{1}&b_{2}\end{array}\right|=18b_{1}b_{2}b_{3}+(b_{1}b_{2})^{2}-4b_{3}b_{1}^{3}-4b_{2}^{3}-27b_{3}^{3}.$$

$$V(t)=c_{1}I(t)+c_{2}C(t)+c_{3}A(t),$$

$$\begin{array}{rl}{c}_1& ={\xi }_1{\xi }_2+{\xi }_1\varphi {\eta }_C+{\xi }_2\rho {\eta }_A,\\ c_2& ={\xi }_1\omega +{\xi }_1{\xi }_3{\eta }_C+\rho {\eta }_A\omega -{\eta }_C\rho \gamma ,\\ c_3& =\gamma {\xi }_2+{\xi }_2{\xi }_3{\eta }_A+\varphi {\eta }_C\gamma -\varphi {\eta }_A\omega .\end{array}$$

$${}_{t_0}^C{D}_t^{\alpha }V=c_1{{}_a^{C}D}_t^{\gamma }I+c_2{}_a^C{D}_t^{\gamma }C+c_3{}_a^C{D}_t^{\gamma }A.$$

$${}_{t_0}^C{D}_t^{\alpha }V=c_1\left(\beta \left(I+{\eta }_{C}C+{\eta }_{A}A\right)S-{\xi }_{3}I+\gamma A+\omega C\right)+c_2\left(\varphi I-{\xi }_{2}C\right)+c_3\left(\rho I-{\xi }_{1}A\right).$$

$${\xi }_1\omega +{\xi }_1{\xi }_3{\eta }_C+\rho {\eta }_A\omega -{\eta }_C\rho \gamma ={\xi }_1\omega +\gamma (\varphi +\mu ){\eta }_C+(\mu +d){\xi }_3{\eta }_C+\rho {\eta }_A\omega >0$$

$$\gamma {\xi }_2+{\xi }_2{\xi }_3{\eta }_A+\varphi {\eta }_C\gamma -\varphi {\eta }_A\omega =\gamma {\xi }_2+\omega (\rho +\mu ){\eta }_A+\mu {\xi }_3{\eta }_A+\varphi {\eta }_C\gamma >0.$$

$$\begin{array}{rl}{}_{t_0}^C{D}_t^{\alpha }V& =({\xi }_1{\xi }_2\beta +{\xi }_1\varphi {\eta }_C\beta +{\xi }_2\rho {\eta }_A\beta )IS+(-{\xi }_1{\xi }_2{\xi }_3+{\xi }_1\omega \varphi +\gamma {\xi }_2\rho )I\\ & +{\eta }_C({\xi }_1{\xi }_2\beta +{\xi }_1\varphi {\eta }_C\beta +{\xi }_2\rho {\eta }_A\beta )CS+{\eta }_C(-{\xi }_1{\xi }_3{\xi }_2+{\xi }_1\varphi \omega +\rho \gamma {\xi }_2)C\\ & +{\eta }_A({\xi }_1{\xi }_2\beta +{\xi }_1\varphi {\eta }_C\beta +{\xi }_2\rho {\eta }_A\beta )AS+{\eta }_A(-{\xi }_2{\xi }_3{\xi }_1+\varphi \omega {\xi }_1+{\xi }_2\rho \gamma )A.\end{array}$$

Full text

††journal: Mathematics and Computers in Simulation

1 Introduction

Fractional differential equations (FDEs), also known in the literature as extraordinary differential equations, are a generalization of differential equations through the application of fractional calculus, that is, the branch of mathematical analysis that studies different possibilities of defining differentiation operators of noninteger order [2, 3]. FDEs are naturally related to systems with memory, which explains their usefulness in most biological systems [4]. Indeed, FDEs have been considered in many epidemiological models. In [5], a fractional order model for nonlocal epidemics is considered, and the results are expected to be relevant to foot-and-mouth disease, SARS, and avian flu. Some necessary and sufficient conditions for local stability of fractional order differential systems are provided [5]. In [6], a fractional order SEIR model with vertical transmission within a nonconstant population is considered, and the asymptotic stability of the disease free and endemic equilibria are analyzed. The stability of the endemic equilibrium of a fractional order SIR model is studied in [7]. A fractional order model of HIV infection of CD4+ T-cells is analyzed in [8]. A fractional order predator prey model and a fractional order rabies model are proposed in [9]. The stability of equilibrium points are studied, and an example is given where the equilibrium point is a centre for the integer order system but locally asymptotically stable for its fractional-order counterpart [9]. A fractional control model for malaria transmission is proposed and studied numerically in [10].

The question of stability for FDEs is crucial: see, e.g., [11, 12] for good overviews on stability of linear/nonlinear, positive, with delay, distributed, and continuous/discrete fractional order systems. In [13], an extension of the Lyapunov direct method for fractional-order systems using Bihari’s and Bellman–Gronwall’s inequality, and a proof of a comparison theorem for fractional-order systems, are obtained. A new lemma for Caputo fractional derivatives, when $0<\alpha<1$, is proposed in [14], which allows to find Lyapunov candidate functions for proving the stability of many fractional order systems, using the fractional-order extension of the Lyapunov direct method. Motivated by the work [14], the authors of [15] extended the Volterra-type Lyapunov function to fractional-order biological systems through an inequality to estimate the Caputo fractional derivatives of order $\alpha\in(0,1)$. Using this result, the uniform asymptotic stability of some Caputo-type epidemic systems with a pair of fractional-order differential equations is proved. Such systems are the basic models of infectious disease dynamics (SIS, SIR and SIRS models) and Ross–Macdonald model for vector-borne diseases. For more on the subject see [16], where the problem of output feedback stabilization for fractional order linear time-invariant systems with fractional commensurate order is investigated, and [17], where the stability of a special observer with a nonlinear weighted function and a transient dynamics function is rigorously analyzed for slowly varying disturbances and higher-order disturbances of fractional-order systems.

Here we propose a Caputo fractional order SICA epidemiological model with constant recruitment rate, mass action incidence and variable population size, for HIV/AIDS transmission. The model is based on an integer-order HIV/AIDS model without memory effects firstly proposed in [18] and later modified in [19, 20]. The model for $\alpha=1$ describes well the clinical reality given by the data of HIV/AIDS infection in Cape Verde from 1987 to 2014 [19]. In the present work, we extend the model by considering fractional differentiation, in order to capture memory effects, long-rage interactions, and hereditary properties, which exist in the process of HIV/AIDS transmission but are neglected in the case $\alpha=1$, that is, for integer-order differentiation [21, 22]. Using the results from [23] and [5], we prove the local asymptotic stability of the disease free equilibrium. Then, we extend the results of [13] and [15] and prove the uniform asymptotic stability of the disease free and endemic equilibrium points. For the numerical implementation of the fractional derivatives, we have used the Adams–Bashforth–Moulton scheme, which has been implemented in the fde12 Matlab routine by Garrappa [24]. The software code implements a predictor-corrector PECE method, as described in [25].

The paper is organized as follows. In Section 2, we present basic definitions and recall necessary results on Caputo fractional calculus and local and uniform asymptotic stability and Volterra-type Lyapunov functions for fractional-order systems. The original results appear in Section 3: we introduce our Caputo fractional-order HIV/AIDS model and study the existence of equilibrium points. More precisely, in Section 3.1 we prove local asymptotic stability of the disease free equilibrium, while in Sections 3.2 and 3.3 we prove uniform asymptotic stability of the disease free and endemic equilibrium points, respectively. We end with Section 4 of numerical simulations, which illustrate the stability results proved in Sections 3.1–3.3.

2 Preliminaries on the Caputo fractional calculus

We begin by introducing the definition of Caputo fractional derivative and recalling its main properties.

Definition 2.1** (See [26]).**

Let $a>0$, $t>a$, and $\alpha,a,t\in\mathbb{R}$. The Caputo fractional derivative of order $\alpha$ of a function $f\in C^{n}$ is given by

[TABLE]

Property 2.1** (Linearity; see, e.g., [27]).**

Let $f,g:[a,b]\rightarrow\mathbb{R}$ be such that ${}_{a}^{C}D_{t}^{\alpha}f(t)$ and ${}_{a}^{C}D_{t}^{\alpha}g(t)$ exist almost everywhere and let $c_{1},c_{2}\in\mathbb{R}$. Then, ${}_{a}^{C}D_{t}^{\alpha}(c_{1}f(t)+c_{2}g(t))$ exists almost everywhere with

[TABLE]

Property 2.2** (Caputo derivative of a constant; see, e.g., [28]).**

The fractional derivative of a constant function $f(t)\equiv c$ is zero:

[TABLE]

Let us consider the following general fractional differential equation involving the Caputo derivative:

[TABLE]

subject to a given initial condition $x_{0}=x(t_{0})$.

Definition 2.2** (See, e.g., [29]).**

The constant $x^{*}$ is an equilibrium point of the Caputo fractional dynamic system (1) if, and only if, $f(t,x^{*})=0$.

Following [23], an equilibrium point $x^{*}$ of the Caputo fractional dynamic system (1) is locally asymptotically stable if all the eigenvalues $\lambda$ of the Jacobian matrix of system (1), evaluated at the equilibrium point $x^{*}$, satisfies the following condition:

[TABLE]

Next theorem gives an extension of the celebrated Lyapunov direct method for Caputo type fractional order nonlinear systems [13].

Theorem 2.3** (Uniform Asymptotic Stability [13]).**

Let $x^{*}$ be an equilibrium point for the nonautonomous fractional order system (1) and $\Omega\subset\mathbb{R}^{n}$ be a domain containing $x^{*}$. Let $L:[0,\infty)\times\Omega\rightarrow\mathbb{R}$ be a continuously differentiable function such that

[TABLE]

and

[TABLE]

for all $\alpha\in(0,1)$ and all $x\in\Omega$, where $W_{1}(\cdot)$, $W_{2}(\cdot)$ and $W_{3}(\cdot)$ are continuous positive definite functions on $\Omega$. Then the equilibrium point $x^{*}$ of system (1) is uniformly asymptotically stable.

In what follows, we recall a lemma proved in [15], where a Volterra-type Lyapunov function is obtained for fractional-order epidemic systems.

Lemma 2.4** (See [15]).**

Let $x(\cdot)$ be a continuous and differentiable function with $x(t)\in\mathbb{R_{+}}$. Then, for any time instant $t\geq t_{0}$, one has

[TABLE]

3 The fractional HIV/AIDS model

In this section we propose a Caputo fractional-order model for HIV/AIDS with memory effects. Our population model assumes a constant recruitment rate, mass action incidence, and variable population size.

The model subdivides human population into four mutually-exclusive compartments: susceptible individuals ($S$); HIV-infected individuals with no clinical symptoms of AIDS (the virus is living or developing in the individuals but without producing symptoms or only mild ones) but able to transmit HIV to others ($I$); HIV-infected individuals under ART treatment (the so called chronic stage) with a viral load remaining low ($C$); and HIV-infected individuals with AIDS clinical symptoms ($A$). The total population at time $t$, denoted by $N(t)$, is given by $N(t)=S(t)+I(t)+C(t)+A(t)$. Effective contact with people infected with HIV is at a rate $\lambda$, given by

[TABLE]

where $\beta$ is the effective contact rate for HIV transmission. The modification parameter $\eta_{A}\geq 1$ accounts for the relative infectiousness of individuals with AIDS symptoms, in comparison to those infected with HIV with no AIDS symptoms. Individuals with AIDS symptoms are more infectious than HIV-infected individuals (pre-AIDS) because they have a higher viral load and there is a positive correlation between viral load and infectiousness. On the other hand, $\eta_{C}\leq 1$ translates the partial restoration of immune function of individuals with HIV infection that use ART correctly. All individuals suffer from natural death, at a constant rate $\mu$. We assume that HIV-infected individuals with and without AIDS symptoms have access to ART treatment. HIV-infected individuals with no AIDS symptoms $I$ progress to the class of individuals with HIV infection under ART treatment $C$ at a rate $\phi$, and HIV-infected individuals with AIDS symptoms are treated for HIV at rate $\gamma$. Individuals in the class $C$ leave to the class $I$ at a rate $\omega$. We also assume that an HIV-infected individual with AIDS symptoms $A$ that starts treatment moves to the class of HIV-infected individuals $I$, moving to the chronic class $C$ only if the treatment is maintained. HIV-infected individuals with no AIDS symptoms $I$ that do not take ART treatment progress to the AIDS class $A$ at rate $\rho$. Note that only HIV-infected individuals with AIDS symptoms $A$ suffer from an AIDS induced death, at a rate $d$. The Caputo fractional-order system that describes the previous assumptions is:

[TABLE]

The biologically feasible region of system (3) is given by

[TABLE]

The model (3) has a disease free equilibrium given by

[TABLE]

Let

[TABLE]

where $\xi_{1}=\gamma+\mu+d$, $\xi_{2}=\omega+\mu$,

[TABLE]

and

[TABLE]

Whenever $R_{0}>1$, the model (3) has a unique endemic equilibrium $\Sigma_{*}=\left(S^{*},I^{*},C^{*},A^{*}\right)$ given by

[TABLE]

3.1 Local asymptotic stability of the disease free equilibrium $\Sigma_{0}$

As firstly proved in [23], stability is guaranteed if and only if the roots of some polynomial (the eigenvalues of the matrix of dynamics or the poles of the transfer matrix) lie outside the closed angular sector $|\arg(\lambda)|\leq\frac{\alpha\pi}{2}$. In our case, the Jacobian matrix $J(\Sigma_{0})$ for system (3) evaluated at the uninfected steady state $\Sigma_{0}$ (5) is given by

[TABLE]

The uninfected steady state is asymptotically stable if all of the eigenvalues $\lambda$ of the Jacobian matrix $J(\Sigma_{0})$ satisfy the following condition (see, e.g., [23]):

[TABLE]

Let $\xi_{3}=\rho+\phi+\mu$. The eigenvalues are determined by solving the characteristic equation $\det(J(\Sigma_{0})−\lambda I)=0$. For $J(\Sigma_{0})$ as in (8), the characteristic equation is given by

[TABLE]

with

[TABLE]

and

[TABLE]

where

[TABLE]

From (9) we have that the eigenvalue $\lambda_{1}=-\mu$ satisfies $|\arg(\lambda_{1})|>\frac{\alpha\pi}{2}$ for all $\alpha\in(0,1)$. The discriminant $D(p)$ of the polynomial (10) is given (see [5]) by

[TABLE]

Following [5], all roots of the polynomial (10) satisfy condition (2) if the following conditions hold:

• (i)

if $D(p)>0$, then the Routh–Hurwitz conditions are a necessary and sufficient condition for the equilibrium point $\Sigma_{0}$ to be locally asymptotically stable, i.e., $b_{1}>0$, $b_{3}>0$ and $b_{1}b_{2}−b_{3}>0$;

• (ii)

if $D(p)<0$, $b_{1}\geq 0$, $b_{2}\geq 0$, $b_{3}>0$, and $\alpha<2/3$, then $\Sigma_{0}$ is locally asymptotically stable;

• (iii)

if $D(p)<0$, $b_{1}<0$, $b_{2}<0$, and $\alpha>2/3$, then $\Sigma_{0}$ is unstable;

• (iv)

if $D(p)<0$, $b_{1}>0$, $b_{2}>0$, and $b_{1}b_{2}−b_{3}=0$, then $\Sigma_{0}$ is locally asymptotically stable for all $\alpha\in[0,1)$;

• (v)

$b_{3}>0$ is a necessary condition for local asymptotic stability of $\Sigma_{0}$.

3.2 Uniform asymptotic stability of the disease free equilibrium $\Sigma_{0}$

In this section, we prove the uniform asymptotic stability of the disease free equilibrium $\Sigma_{0}$ (5) of the fractional order system (3).

Theorem 3.1**.**

Let $\alpha\in(0,1)$. The disease free equilibrium $\Sigma_{0}$ (5), of the fractional system (3), is uniformly asymptotically stable in $\Omega$ (4), whenever (6) satisfies $R_{0}<1$.

Proof.

Consider the following Lyapunov function:

[TABLE]

where

[TABLE]

Function $V$ is defined, continuous and positive definite for all $I(t)>0$, $C(t)>0$ and $A(t)>0$. By Property 2.1, we have

[TABLE]

From (3) we have

[TABLE]

Note that

[TABLE]

and

[TABLE]

Therefore, we have

[TABLE]

As $S\leq S^{0}$,

[TABLE]

holds. From $S^{0}\left(\xi_{1}\xi_{2}\beta+\xi_{1}\phi\eta_{C}\beta+\xi_{2}\rho\eta_{A}\beta\right)=\mathcal{N}$ and $-\xi_{1}\xi_{2}\xi_{3}+\xi_{1}\omega\phi+\gamma\xi_{2}\rho=-\mathcal{D}$, one has

[TABLE]

Because all the model parameters are nonnegative, it follows that ${}_{t_{0}}^{C}D_{t}^{\alpha}V\leq 0$ for $R_{0}<1$ with ${}_{t_{0}}^{C}D_{t}^{\alpha}V=0$ if, and only if, $I=C=A=0$. Substituting $(I,C,A)=(0,0,0)$ in (3) shows that $S\to S^{0}=\frac{\Lambda}{\mu}$ as $t\to\infty$. Hence, by Theorem 2.3, the equilibrium point $\Sigma_{0}$ of system (3) is uniformly asymptotically stable in $\Omega$, whenever $R_{0}<1$. ∎

3.3 Uniform asymptotic stability of the endemic equilibrium $\Sigma_{*}$

In this section we prove uniform asymptotic stability of the endemic equilibrium $\Sigma_{*}$ (7) of the fractional order system (3).

Theorem 3.2**.**

Let $\alpha\in(0,1)$ and (6) be such that $R_{0}>1$. Then the unique endemic equilibrium $\Sigma_{*}$ (7) of the fractional order system (3) is uniformly asymptotically stable in the interior of $\Omega$ (4).

Proof.

Consider the following function:

[TABLE]

where

[TABLE]

Function $V$ is a Lyapunov function because it is defined, continuous, and positive definite for all $S(t)>0$, $I(t)>0$, $C(t)>0$ and $A(t)>0$. By Lemma 2.4, we have

[TABLE]

It follows from (3) that

[TABLE]

Using the relation $\Lambda=\beta\left(I^{*}+\eta_{C}\,C^{*}+\eta_{A}A^{*}\right)S^{*}+\mu S^{*}$, we have from the first equation of system (3) at steady-state that (11) can be written as

[TABLE]

which can then be simplified to

[TABLE]

Using the relations at the steady state,

[TABLE]

and, after some simplifications, we have

[TABLE]

The terms between the larger brackets are less than or equal to zero by the well-known inequality that asserts the geometric mean to be less than or equal to the arithmetic mean. Therefore, ${}_{t_{0}}^{C}D_{t}^{\alpha}V(S,I,C,A)$ is negative definite when $0<\alpha<1$. By Theorem 2.3 (the uniform asymptotic stability theorem), the endemic equilibrium $\Sigma_{*}$ is uniformly asymptotically stable in the interior of $\Omega$, whenever $R_{0}>1$. ∎

Note that the fractional model (3) is stable independently of the parameter values. Indeed, the values of the parameters determine the value of $R_{0}$ and, for $R_{0}<1$, the stability of the system is, according with Theorem 3.1, “around” the disease free equilibrium $\Sigma_{0}$; for $R_{0}>1$, the stability of the system is, in agreement with Theorem 3.2, “around” the endemic equilibrium $\Sigma_{*}$.

4 Numerical simulations

In this section we study the dynamical behavior of our model (3), by variation of the noninteger order derivative $\alpha$.

4.1 Local asymptotic stability of the disease free equilibrium $\Sigma_{0}$

Consider the parameter values of Table 1 and $\beta=0.001$. The basic reproduction number (6) is

[TABLE]

while the disease free equilibrium (5) takes the value

[TABLE]

On the other hand, the discriminant $D(p)$ of the polynomial $p$ (10) is given by $D(p)=0.51045>0$, $b_{1}=2.41711>0$, $b_{3}=0.00652>0$ and $b_{1}b_{2}b_{3}=0.02205>0$. Therefore, the Routh–Hurwitz conditions are a necessary and sufficient condition for the equilibrium point $\Sigma_{0}$ to be locally asymptotically stable (see Section 3.1). The stability of the disease free equilibrium $\Sigma_{0}$ is illustrated in Figure 1, where we considered the initial conditions

[TABLE]

and a fixed time step size of $h=2^{-6}$.

For the numerical implementation of the fractional derivatives, we have used the Adams–Bashforth–Moulton scheme, which has been implemented in the Matlab code fde12 by Garrappa [24]. This code implements a predictor-corrector PECE method of Adams–Bashforth–Moulton type, as described in [25].

Regarding convergence and accuracy of the numerical method, we refer to [31]. The stability properties of the method implemented by fde12 have been studied in [32]. Here we considered, without loss of generality, the fractional-order derivatives $\alpha=1.0,0.9,0.8$ and $0.7$.

4.2 Stability of the endemic equilibrium $\Sigma_{*}$

For the numerical study of the stability of the endemic equilibrium $\Sigma_{*}$ (7), we consider the parameter values from Table 1 and $\beta=0.01$. The basic reproduction number (6) takes the value $R_{0}=7.95871$. The concrete value of the endemic equilibrium (7) is $\Sigma_{*}=\left(18.3490,8.0673,77.2881,0.6001\right)$. Figure 2 illustrates the stability of the endemic equilibrium for the initial conditions

[TABLE]

where a fixed time step size of $h=2^{-6}$ has been used.

Our results show that the smaller the order $\alpha$ of the fractional derivative, the slower the convergence to the equilibrium point.

Acknowledgments

This research was partially supported by the Portuguese Foundation for Science and Technology (FCT) through the R&D unit CIDMA, reference UID/MAT/04106/2019, and by project PTDC/EEI-AUT/2933/2014 (TOCCATA), funded by FEDER funds through COMPETE 2020 – Programa Operacional Competitividade e Internacionalização (POCI) and by national funds through FCT. Silva is also supported by national funds (OE), through FCT, I.P., in the scope of the framework contract foreseen in the numbers 4, 5 and 6 of the article 23, of the Decree-Law 57/2016, of August 29, changed by Law 57/2017, of July 19. The authors are grateful to three reviewers for their critical remarks and precious suggestions, which helped them to improve the quality and clarity of the manuscript.