Global Stability for a HIV/AIDS Model

TL;DR

This paper analyzes the global stability of a HIV/AIDS population model with constant recruitment and mass action incidence, establishing conditions for disease-free and endemic equilibrium stability using Lyapunov methods.

Contribution

It provides new theoretical results on the existence, uniqueness, and global stability of equilibria in a HIV/AIDS model with variable population size.

Findings

Existence and uniqueness of disease-free and endemic equilibria.

Global stability of these equilibria established.

Conditions under which the disease persists or dies out.

Abstract

We investigate global stability properties of a HIV/AIDS population model with constant recruitment rate, mass action incidence, and variable population size. Existence and uniqueness results for disease-free and endemic equilibrium points are proved. Global stability of the equilibria is obtained through Lyapunov's direct method and LaSalle's invariance principle.

Equations

$$N(t)=S(t)+I(t)+C(t)+A(t).$$

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

$$\{\begin{array}{l}\dot{S}(t)=\Lambda -\beta \left(I(t)+{\eta }_{C}C(t)+{\eta }_{A}A(t)\right)S(t)-\mu S(t),\\ \dot{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)+\alpha A(t),\\ \dot{C}(t)=\varphi I(t)-(\omega +\mu )C(t),\\ \dot{A}(t)=\rho I(t)-(\alpha +\mu +d)A(t).\end{array}$$

$$\dot{N}(t)=\Lambda -\mu N(t)-dA(t).$$

$$\Omega =\left\{\left(S,I,C,A\right)\in {\mathbb{R}}_{+}^4:N\le \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{D}=\mu \left({\xi }_2\left(\rho +{\xi }_1\right)+\varphi {\xi }_1+\rho d\right)+\rho \omega d.$$

$$\begin{array}{rl}V& =\left({\xi }_1{\xi }_2+{\xi }_1\varphi {\eta }_C+{\xi }_2\rho {\eta }_A\right)I+\left({\xi }_1\omega +{\xi }_1{\xi }_3{\eta }_C+\rho {\eta }_A\omega -{\eta }_C\rho \alpha \right)C\\ & +\left(\alpha {\xi }_2+{\xi }_2{\xi }_3{\eta }_A+\varphi {\eta }_C\alpha -\varphi {\eta }_A\omega \right)A.\end{array}$$

$$\begin{array}{rl}\dot{V}& =\left({\xi }_1{\xi }_2+{\xi }_1\varphi {\eta }_C+{\xi }_2\rho {\eta }_A\right)\dot{I}+\left({\xi }_1\omega +{\xi }_1{\xi }_3{\eta }_C+\rho {\eta }_A\omega -{\eta }_C\rho \alpha \right)\dot{C}\\ & +\left(\alpha {\xi }_2+{\xi }_2{\xi }_3{\eta }_A+\varphi {\eta }_C\alpha -\varphi {\eta }_A\omega \right)\dot{A}\\ & =\left({\xi }_1{\xi }_2+{\xi }_1\varphi {\eta }_C+{\xi }_2\rho {\eta }_A\right)\left(\beta \left(I+{\eta }_{C}C+{\eta }_{A}A\right)S-{\xi }_{3}I+\alpha A+\omega C\right)\\ & +\left({\xi }_1\omega +{\xi }_1{\xi }_3{\eta }_C+\rho {\eta }_A\omega -{\eta }_C\rho \alpha \right)\left(\varphi I-{\xi }_{2}C\right)\\ & +\left(\alpha {\xi }_2+{\xi }_2{\xi }_3{\eta }_A+\varphi {\eta }_C\alpha -\varphi {\eta }_A\omega \right)\left(\rho I-{\xi }_{1}A\right),\end{array}$$

$$\begin{array}{rl}\dot{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 +\alpha {\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 \alpha {\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 \alpha )A.\end{array}$$

$$\begin{array}{rl}\dot{V}& \le ({\xi }_1{\xi }_2\beta +{\xi }_1\varphi {\eta }_C\beta +{\xi }_2\rho {\eta }_A\beta )I{S}^0+(-{\xi }_1{\xi }_2{\xi }_3+{\xi }_1\omega \varphi +\alpha {\xi }_2\rho )I\\ & +{\eta }_C({\xi }_1{\xi }_2\beta +{\xi }_1\varphi {\eta }_C\beta +{\xi }_2\rho {\eta }_A\beta )C{S}^0+{\eta }_C\left(-{\xi }_1{\xi }_3{\xi }_2+{\xi }_1\varphi \omega +\rho \alpha {\xi }_2\right)C\\ & +{\eta }_A({\xi }_1{\xi }_2\beta +{\xi }_1\varphi {\eta }_C\beta +{\xi }_2\rho {\eta }_A\beta )A{S}^0+{\eta }_A\left(-{\xi }_2{\xi }_3{\xi }_1+\varphi \omega {\xi }_1+{\xi }_2\rho \alpha \right)A.\end{array}$$

$$\begin{array}{rl}\dot{V}& \le \mathcal{N}I-\mathcal{D}I+{\eta }_C(\mathcal{N}C-\mathcal{D}C)+{\eta }_A(\mathcal{N}A-\mathcal{D}A)\\ & =\mathcal{D}I\left(R_0-1\right)+{\eta }_C\mathcal{D}C\left(R_0-1\right)+{\eta }_A\mathcal{D}A\left(R_0-1\right)\\ & \le 0\text{for}{R}_0<1.\end{array}$$

$${\Sigma }_{+}=(S^{\ast },I^{\ast },C^{\ast },A^{\ast })$$

$$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(\Lambda \mathcal{N}-\mu \mathcal{D})}{\mathcal{DN}}.$$

$$V=\left(S-S^{\ast }\ln (S)\right)+\left(I-I^{\ast }\ln (I)\right)+\frac{\omega }{{\xi }_2}\left(C-C^{\ast }\ln (C)\right)+\frac{\alpha }{{\xi }_1}\left(A-A^{\ast }\ln (A)\right).$$

$$\dot{V}=\left(1-\frac{S^{\ast }}{S}\right)\dot{S}+\left(1-\frac{I^{\ast }}{I}\right)\dot{I}+\frac{\omega }{{\xi }_2}\left(1-\frac{C^{\ast }}{C}\right)\dot{C}+\frac{\alpha }{{\xi }_1}\left(1-\frac{A^{\ast }}{A}\right)\dot{A}.$$

$$\begin{array}{rl}\dot{V}& =\left(1-\frac{S^{\ast }}{S}\right)\left[\Lambda -\beta \left(I+{\eta }_{C}C+{\eta }_{A}A\right)S-\mu S\right]\\ & +\left(1-\frac{I^{\ast }}{I}\right)\left[\beta \left(I+{\eta }_{C}C+{\eta }_{A}A\right)S-{\xi }_{3}I+\alpha A+\omega C\right]\\ & +\frac{\omega }{{\xi }_2}\left(1-\frac{C^{\ast }}{C}\right)\left[\varphi I-{\xi }_{2}C\right]+\frac{\alpha }{{\xi }_1}\left(1-\frac{A^{\ast }}{A}\right)\left[\rho I-{\xi }_{1}A\right].\end{array}$$

$$\begin{array}{rl}\dot{V}& =\left(1-\frac{S^{\ast }}{S}\right)\left[\beta \left(I^{\ast }+{\eta }_C{C}^{\ast }+{\eta }_A{A}^{\ast }\right)S^{\ast }+\mu S^{\ast }-\beta \left(I+{\eta }_{C}C+{\eta }_{A}A\right)S-\mu S\right]\\ & +\left(1-\frac{I^{\ast }}{I}\right)\left[\beta \left(I+{\eta }_{C}C+{\eta }_{A}A\right)S-{\xi }_{3}I+\alpha A+\omega C\right]\\ & +\frac{\omega }{{\xi }_2}\left(1-\frac{C^{\ast }}{C}\right)\left[\varphi I-{\xi }_{2}C\right]+\frac{\alpha }{{\xi }_1}\left(1-\frac{A^{\ast }}{A}\right)\left[\rho I-{\xi }_{1}A\right],\end{array}$$

$$\begin{array}{rl}\dot{V}& =\left(1-\frac{S^{\ast }}{S}\right)\beta I^{\ast }{S}^{\ast }+\mu S^{\ast }\left(2-\frac{S}{S^{\ast }}-\frac{S^{\ast }}{S}\right)-\beta IS+\beta I{S}^{\ast }\\ & +\beta ({\eta }_C{C}^{\ast }+{\eta }_A{A}^{\ast })S^{\ast }-\beta ({\eta }_{C}C+{\eta }_{A}A)S-\frac{S^{\ast }}{S}\beta ({\eta }_C{C}^{\ast }+{\eta }_A{A}^{\ast })S^{\ast }\\ & +S^{\ast }\beta ({\eta }_{C}C+{\eta }_{A}A)+\left(1-\frac{I^{\ast }}{I}\right)\left[\beta \left(I+{\eta }_{C}C+{\eta }_{A}A\right)S-{\xi }_{3}I+\alpha A+\omega C\right]\\ & +\frac{\omega }{{\xi }_2}\left(1-\frac{C^{\ast }}{C}\right)\left[\varphi I-{\xi }_{2}C\right]+\frac{\alpha }{{\xi }_1}\left(1-\frac{A^{\ast }}{A}\right)\left[\rho I-{\xi }_{1}A\right].\end{array}$$

$${\xi }_3{I}^{\ast }=\beta (I^{\ast }+{\eta }_C{C}^{\ast }+{\eta }_A{A}^{\ast })S^{\ast }+\alpha A^{\ast }+\omega C^{\ast },{\xi }_2{C}^{\ast }=\varphi I^{\ast },{\xi }_1{A}^{\ast }=\rho I^{\ast },$$

$$\begin{array}{rl}\dot{V}& =\left(\beta I^{\ast }{S}^{\ast }+\mu S^{\ast }\right)\left(2-\frac{S}{S^{\ast }}-\frac{S^{\ast }}{S}\right)+\beta S^{\ast }\left({\eta }_C{C}^{\ast }+{\eta }_A{A}^{\ast }\right)\left(2-\frac{S^{\ast }}{S}-\frac{I}{I^{\ast }}\right)\\ & +\beta S^{\ast }\left({\eta }_{C}C+{\eta }_{A}A\right)\left(1-\frac{I^{\ast }}{I}\frac{S}{S^{\ast }}\right)+\alpha A^{\ast }\left(1-\frac{A}{A^{\ast }}\frac{I^{\ast }}{I}\right)+\omega C^{\ast }\left(1-\frac{C}{C^{\ast }}\frac{I^{\ast }}{I}\right)\\ & +\frac{\omega \varphi }{{\xi }_2}{I}^{\ast }\left(1-\frac{I}{I^{\ast }}\frac{C^{\ast }}{C}\right)+\frac{\alpha \rho }{{\xi }_1}{I}^{\ast }\left(1-\frac{I}{I^{\ast }}\frac{A^{\ast }}{A}\right).\end{array}$$

Full text

1. Introduction

Mathematical models may represent a useful tool in the development of public health policies [8, 20, 21]. Although it is unlikely that a mathematical model will provide accurate long-term predictions on the number of AIDS cases, one such model, based on interactions that lead to disease transmission, could eventually allow researchers to answer many useful questions [13]. As a result, several mathematical models have been proposed in the last decades for HIV/AIDS transmission dynamics: see, e.g., [1, 2, 3, 4, 11, 17, 19, 22] and references cited therein.

Global stability of equilibrium points for mathematical models of HIV/AIDS transmission dynamics has been studied by different authors: see, e.g., [5, 9, 12]. In [12], the authors consider different latent stages depending on other chronic diseases that each individual may have. The epidemic model in [16] considers a latent stage and vaccination of newborns and susceptible. In [18], it is assumed that the HIV epidemic spreads both through horizontal and vertical transmission; in [23], the immigration of infective individuals is considered, both models with a variable size population. The effect of screening unaware infective individuals on the spread of HIV, in a constant population, is considered in the mathematical model proposed in [25]. In [5], the global stability is studied for a HIV/AIDS model with two infective stages and where a discrete time delay is introduced, describing the time from start of treatment in the symptomatic stage until treatment effects become visible.

Motivated by the results of [24], in this paper we propose a mathematical model for HIV/AIDS transmission with varying population size in a homogeneously mixing population. Differently from [24], here we consider a mass action hypothesis for the transmission rate. We assume that the rate at which susceptible are infected by individuals with AIDS symptoms is bigger or equal than the rate of infection by contact with HIV-infected individuals (pre-AIDS). This is justifiable because individuals with AIDS symptoms have a higher viral load and it is known that there exists a positive correlation between viral load and infectiousness [6]. On the other hand, individuals with HIV-infection under anti-retroviral treatment (ART) suffer a partial restoration of the immune function and, therefore, we assume that the rate of infection by contact with individuals under ART is smaller or equal than the rate of infection by contact with HIV-infected individuals (pre-AIDS), which are not under ART (see, e.g., [7]). We prove the global stability of the disease free equilibrium whenever the basic reproduction number $R_{0}$ is less than one; and the global stability of the unique endemic equilibrium when $R_{0}$ is greater than one. The global stability analysis is done through Lyapunov’s direct method combined with LaSalle’s invariance principle.

The paper is organized as follows. In Section 2, we describe the mathematical model for HIV/AIDS transmission. Then, in Section 3, we prove existence and global stability of the disease free equilibrium. The existence and global stability of the unique endemic equilibrium point is proved in Section 4. The stability results are then illustrated through numerical simulations in Section 5. We finish the paper with Section 6 of concluding remarks.

2. Model for HIV/AIDS transmission

In this paper, we propose and analyze a mathematical model for HIV/AIDS transmission with varying population size in a homogeneously mixing population. The model is based on that of [24], and subdivides the 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 other individuals ($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

[TABLE]

The effective contact with people infected with HIV is at a rate $\lambda$, given by

[TABLE]

where $\beta$ is the 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 and 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 [6]. On the other hand, $\eta_{C}\leq 1$ translates the partial restoration of the immune function of individuals with HIV infection that are correctly treated under ART [7].

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 a rate $\alpha$. We assume that HIV-infected individuals with AIDS symptoms, $A$, that start treatment, move to the class of HIV-infected individuals, $I$, and will move 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$. We assume that only HIV-infected individuals with AIDS symptoms, $A$, suffer from an AIDS induced death, at a rate $d$. These assumptions are translated into the following mathematical model:

[TABLE]

From $N(t)=S(t)+I(t)+C(t)+A(t)$ and (2.1), it follows that

[TABLE]

Thus, the total population size $N$ may vary in time. Let $\Omega$ denote the biologically feasible region

[TABLE]

Using a standard comparison theorem (see [14]), one can easily show that $N(t)\leq\frac{\Lambda}{\mu}$ if $N(0)\leq\frac{\Lambda}{\mu}$. Thus, the region $\Omega$ is positively invariant. Hence, it is sufficient to consider the dynamics of the flow generated by (2.1) in $\Omega$. In this region, the model is epidemiologically and mathematically well posed in the sense of [10]. In other words, every solution of the model (2.1) with initial conditions in $\Omega$ remains in $\Omega$ for all $t>0$. Therefore, the dynamics of our model will be considered in $\Omega$.

3. Existence and global stability of the DFE

Model (2.1) has a disease-free equilibrium (DFE) given by

[TABLE]

Following [26], the basic reproduction number $R_{0}$ for (2.1), which represents the expected average number of new HIV infections produced by a single HIV-infected individual when in contact with a completely susceptible population, is given by

[TABLE]

where $\xi_{1}=\alpha+\mu+d$, $\xi_{2}=\omega+\mu$, $\mathcal{N}=\beta\,\left(\xi_{2}\left(\xi_{1}+\rho\,\eta_{A}\right)+\eta_{C}\,\phi\,\xi_{1}\right)$, and

[TABLE]

The following local stability result follows easily from Theorem 2 of [26].

Lemma 1**.**

The disease free equilibrium $\Sigma_{0}$ is locally asymptotically stable if $R_{0}<1$ and unstable if $R_{0}>1$.

Now we prove the global stability of the disease free equilibrium (3.1).

Theorem 1**.**

The disease free equilibrium $\Sigma_{0}$ is globally asymptotically stable for $R_{0}<1$.

Proof.

Let $\xi_{3}=\rho+\phi+\mu$. Consider the following Lyapunov function:

[TABLE]

Note that $\xi_{1}\omega+\xi_{1}\xi_{3}\eta_{C}+\rho\eta_{A}\omega-\eta_{C}\rho\alpha=\xi_{1}\omega+\alpha(\phi+\mu)\eta_{C}+(\mu+d)\xi_{3}\eta_{C}+\rho\eta_{A}\omega>0$ and $\alpha\xi_{2}+\xi_{2}\xi_{3}\eta_{A}+\phi\eta_{C}\alpha-\phi\eta_{A}\omega=\alpha\xi_{2}+\omega(\rho+\mu)\eta_{A}+\mu\xi_{3}\eta_{A}+\phi\eta_{C}\alpha>0$. The time derivative of $V$ computed along the solutions of (2.1) is given by

[TABLE]

which can be further simplified to

[TABLE]

As $S\leq S^{0}$, the following inequality holds:

[TABLE]

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+\alpha\xi_{2}\rho=-\mathcal{D}$,

[TABLE]

Because all model parameters are nonnegative, it follows that $\dot{V}\leq 0$, for $R_{0}<1$ with $\dot{V}=0$, if and only if $I=C=A=0$. Substituting $(I,C,A)=(0,0,0)$ into (2.1) shows that $S\to S^{0}=\frac{\Lambda}{\mu}$ as $t\to\infty$. Hence, $V$ is a Lyapunov function on $\Omega$ and the largest compact invariant set in $\{(S,I,C,A)\in\Omega\,:\,\dot{V}=0\}$ is the singleton $\{\Sigma_{0}\}$. Thus, by LaSalle’s invariance principle [15], every solution of (2.1), with initial conditions in $\Omega$, approaches $\Sigma_{0}$ as $t\to\infty$, whenever $R_{0}<1$. ∎

4. Existence and global stability of the endemic equilibrium

It is easy to show that model (2.1) has a unique endemic equilibrium

[TABLE]

whenever $R_{0}>1$. This is precisely stated in Lemma 2.

Lemma 2**.**

The model (2.1) has a unique endemic equilibrium $\Sigma_{+}=(S^{*},I^{*},C^{*},A^{*})$ whenever $R_{0}>1$, which is given by

[TABLE]

We now prove the global stability of the endemic equilibrium $\Sigma_{+}$.

Theorem 2**.**

The endemic equilibrium $\Sigma_{+}$ of model (2.1) is globally asymptotically stable for $R_{0}>1$.

Proof.

We start by defining the region $\Omega_{0}=\left\{(S,I,C,A)\in\Omega\,|\,I=C=A=0\right\}$. Consider the following Lyapunov function:

[TABLE]

Differentiating $V$ with respect to time gives

[TABLE]

Substituting the expressions for the derivatives in $\dot{V}$, it follows from (2.1) 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 (2.1) at steady-state that (4.1) 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]

Because the geometric mean is less or equal than the arithmetic mean, it follows that the terms between the larger brackets are less or equal than zero and $\dot{V}=0$ holds if and only if $(S,I,C,A)$ take the equilibrium values $(S^{*},I^{*},C^{*},A^{*})$. Thus, by LaSalle’s invariance principle, the endemic equilibrium $\Sigma_{+}$ is globally asymptotically stable. ∎

5. Numerical simulations

In this section, we provide some numerical simulations that illustrate the analytic results proved in Sections 3 and 4. Consider the parameter values $\mu=1/70$, $\Lambda=2$, $\beta=0.001$, $\eta_{C}=0.04$, $\eta_{A}=1.3$, $\omega=0.09$, $\rho=0.1$, $\phi=1$, $\alpha=0.33$ and $d=1$. The corresponding basic reproduction number is equal to $R_{0}=0.9141$. The disease free equilibrium is given by $\left(S^{0},I^{0},C^{0},A^{0}\right)=(140,0,0,0)$. Figure 1 illustrates the stability of the disease free equilibrium proved in Theorem 1.

In Figure 2, we can observe the stability of the endemic equilibrium proved in Theorem 2 for the paremeter values $\mu=1/70$, $\Lambda=2$, $\beta=0.002$, $\eta_{C}=0.04$, $\eta_{A}=1.3$, $\omega=0.09$, $\rho=0.1$, $\phi=1$, $\alpha=0.33$ and $d=1$, which corresponds to a basic reproduction number equal to $R_{0}=1.8281$ and where the unique endemic equilibrium is given by $\Sigma_{+}=(S^{*},I^{*},C^{*},A^{*})=(76.5820,3.9959,38.3171,0.2973)$.

6. Conclusion

We proposed a mathematical model for HIV/AIDS transmission with variable total population size and different transmission rates depending on the viral load of HIV infected individuals. We proved existence of a disease free equilibrium and computed the basic reproduction number $R_{0}$ using the method in [26]. Existence of an endemic equilibrium is proved for $R_{0}>1$. We also proved the global stability of the disease free equilibrium when $R_{0}<1$ and the global stability of the endemic equilibrium for $R_{0}>1$. The proofs of global stability are carried out through Lyapunov’s direct method combined with LaSalle’s invariance principle. The numerical simulations provided in Section 5 illustrate the obtained stability results.

Acknowledgments

This research was partially supported by FCT (The Portuguese Foundation for Science and Technology) within projects UID/MAT/04106/2013, CIDMA, and 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 grateful to the FCT post-doc fellowship SFRH/BPD/72061/2010.