Modeling and optimal control of HIV/AIDS prevention through PrEP

TL;DR

This paper develops a new epidemiological model for HIV/AIDS incorporating PrEP, proves its mathematical properties, calibrates it with real data from Cape Verde, and determines optimal PrEP strategies to reduce HIV transmission effectively.

Contribution

It introduces a novel HIV/AIDS model with PrEP, proves its stability, calibrates it with real data, and derives optimal control strategies considering cost and coverage constraints.

Findings

PrEP significantly reduces HIV transmission in simulations.

The model accurately predicts HIV/AIDS cases in Cape Verde.

Optimal strategies balance PrEP costs and epidemic control.

Abstract

Pre-exposure prophylaxis (PrEP) consists in the use of an antiretroviral medication to prevent the acquisition of HIV infection by uninfected individuals and has recently demonstrated to be highly efficacious for HIV prevention. We propose a new epidemiological model for HIV/AIDS transmission including PrEP. Existence, uniqueness and global stability of the disease free and endemic equilibriums are proved. The model with no PrEP is calibrated with the cumulative cases of infection by HIV and AIDS reported in Cape Verde from 1987 to 2014, showing that it predicts well such reality. An optimal control problem with a mixed state control constraint is then proposed and analyzed, where the control function represents the PrEP strategy and the mixed constraint models the fact that, due to PrEP costs, epidemic context and program coverage, the number of individuals under PrEP is limited at each instant of time. The objective is to determine the PrEP strategy that satisfies the mixed state control constraint and minimizes the number of individuals with pre-AIDS HIV-infection as well as the costs associated with PrEP. The optimal control problem is studied analytically. Through numerical simulations, we demonstrate that PrEP reduces HIV transmission significantly.

Tables

Table 1. Table 1. Cumulative cases of infection by HIV/AIDS and total population in Cape Verde in the period 1987–2014 [ 25 , 42 ] .

Year	1987	1988	1989	1990	1991	1992	1993
HIV/AIDS	61	107	160	211	244	303	337
Population	323972	328861	334473	341256	349326	358473	368423
Year	1994	1995	1996	1997	1998	1999	2000
HIV/AIDS	358	395	432	471	560	660	779
Population	378763	389156	399508	409805	419884	429576	438737
Year	2001	2002	2003	2004	2005	2006	2007
HIV/AIDS	913	1064	1233	1493	1716	2015	2334
Population	447357	455396	462675	468985	474224	478265	481278
Year	2008	2009	2010	2011	2012	2013	2014
HIV/AIDS	2610	2929	3340	3739	4090	4537	4946
Population	483824	486673	490379	495159	500870	507258	513906

Table 2. Table 2. Parameters of the HIV/AIDS model ( 1 ) for Cape Verde.

Symbol	Description	Value	Reference
$N\left(0\right)$	Initial population	$323972$	[38]
$\Lambda$	Recruitment rate	$13045$	[38]
$\mu$	Natural death rate	$1/69.54$	[38]
$\beta$	HIV transmission rate	$0.752$	Estimated
${\eta }_C$	Modification parameter	$0.015$, $0.04$	Assumed
${\eta }_A$	Modification parameter	$1.3$, $1.35$	Assumed
$\varphi$	HIV treatment rate for $I$ individuals	$1$	[30]
$\rho$	Default treatment rate for $I$ individuals	$0.1$	[30]
$\alpha$	AIDS treatment rate	$0.33$	[30]
$\omega$	Default treatment rate for $C$ individuals	$0.09$	[30]
$d$	AIDS induced death rate	$1$	[39]

Equations

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

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

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

$$\underset{t\to \infty }{\mathrm{liminf}}N(t)\ge \epsilon$$

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

$$R_0=\frac{\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{\mathcal{N}}{\mathcal{D}},$$

$${\lambda }^{\ast }=\frac{\beta \left(I^{\ast }+{\eta }_C{C}^{\ast }+{\eta }_A{A}^{\ast }\right)}{N^{\ast }}.$$

$$S^{\ast }=\frac{\Lambda }{{\lambda }^{\ast }+\mu },I^{\ast }=-\frac{{\lambda }^{\ast }\Lambda {\xi }_1{\xi }_2}{D},C^{\ast }=-\frac{\varphi {\lambda }^{\ast }\Lambda {\xi }_1}{D},A^{\ast }=-\frac{{\rho }_1{\lambda }^{\ast }\Lambda {\xi }_2}{D}$$

$${\lambda }^{\ast }=-\mu (1-R_0).$$

$$\begin{array}{rl}{S}^{\ast }& =\frac{\Lambda (\rho d{\xi }_2-\mathcal{D})}{\mu (\rho d{\xi }_2-\mathcal{N})},I^{\ast }=\frac{\Lambda {\xi }_1{\xi }_2(\mathcal{D}-\mathcal{N})}{\mathcal{D}(\rho d{\xi }_2-\mathcal{N})},\\ C^{\ast }& =\frac{\Lambda \varphi {\xi }_1(\mathcal{D}-\mathcal{N})}{\mathcal{D}(\rho d{\xi }_2-\mathcal{N})},A^{\ast }=\frac{\Lambda \rho {\xi }_2(\mathcal{D}-\mathcal{N})}{\mathcal{D}(\rho d{\xi }_2-\mathcal{N})}.\end{array}$$

$$\lambda ={\beta }_1\left(I+{\eta }_{C}C+{\eta }_{A}A\right),\text{where}{\beta }_1=\frac{\beta \mu }{\Lambda }.$$

$$\{\begin{array}{l}\dot{S}(t)=\Lambda -{\beta }_1\left(I(t)+{\eta }_{C}C(t)+{\eta }_{A}A(t)\right)S(t)-\mu S(t),\\ \dot{I}(t)={\beta }_1\left(I(t)+{\eta }_{C}C(t)+{\eta }_{A}A(t)\right)S(t)-{\xi }_{3}I(t)+\alpha A(t)+\omega C(t),\\ \dot{C}(t)=\varphi I(t)-{\xi }_{2}C(t),\\ \dot{A}(t)=\rho I(t)-{\xi }_{1}A(t),\end{array}$$

$${\Omega }_0=\{\left(S,I,C,A\right)\in \Omega :I=C=A=0\}.$$

$$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 }_1\left(I+{\eta }_{C}C+{\eta }_{A}A\right)S-\mu S\right]\\ & +\left(1-\frac{I^{\ast }}{I}\right)\left[{\beta }_1\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 }_1\left(I^{\ast }+{\eta }_C{C}^{\ast }+{\eta }_A{A}^{\ast }\right)S^{\ast }+\mu S^{\ast }-{\beta }_1\left(I+{\eta }_{C}C+{\eta }_{A}A\right)S-\mu S\right]\\ & +\left(1-\frac{I^{\ast }}{I}\right)\left[{\beta }_1\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 }_1{I}^{\ast }{S}^{\ast }+\mu S^{\ast }\left(2-\frac{S}{S^{\ast }}-\frac{S^{\ast }}{S}\right)-{\beta }_{1}IS+{\beta }_{1}I{S}^{\ast }\\ & +{\beta }_1({\eta }_C{C}^{\ast }+{\eta }_A{A}^{\ast })S^{\ast }-{\beta }_1({\eta }_{C}C+{\eta }_{A}A)S-\frac{S^{\ast }}{S}{\beta }_1({\eta }_C{C}^{\ast }+{\eta }_A{A}^{\ast })S^{\ast }\\ & +S^{\ast }{\beta }_1({\eta }_{C}C+{\eta }_{A}A)+\left(1-\frac{I^{\ast }}{I}\right)\left[{\beta }_1\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 }_1(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 }_1{I}^{\ast }{S}^{\ast }+\mu S^{\ast }\right)\left(2-\frac{S}{S^{\ast }}-\frac{S^{\ast }}{S}\right)+{\beta }_1{S}^{\ast }\left({\eta }_C{C}^{\ast }+{\eta }_A{A}^{\ast }\right)\left(1-\frac{S^{\ast }}{S}\right)\\ & +{\beta }_1{S}^{\ast }\left({\eta }_C{C}^{\ast }+{\eta }_A{A}^{\ast }\right)\left(1-\frac{I}{I^{\ast }}\right)+{\beta }_1{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}$$

$$S_0=S(0)=323911,I_0=I(0)=61,C_0=C(0)=0,A_0=A(0)=0.$$

$$\{\begin{array}{l}\dot{S}(t)=\Lambda -\frac{\beta \left(I(t)+{\eta }_{C}C(t)+{\eta }_{A}A(t)\right)}{N(t)}S(t)-\mu S(t)-\psi S(t)+\theta E(t),\\ \dot{I}(t)=\frac{\beta \left(I(t)+{\eta }_{C}C(t)+{\eta }_{A}A(t)\right)}{N(t)}S(t)-(\rho +\varphi +\mu )I(t)+\alpha A(t)+\omega C(t),\\ \dot{C}(t)=\varphi I(t)-(\omega +\mu )C(t),\\ \dot{A}(t)=\rho I(t)-(\alpha +\mu +d)A(t),\\ \dot{E}(t)=\psi S(t)-(\mu +\theta )E(t).\end{array}$$

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

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

$$R_0=\frac{\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{\mathcal{N}}{\mathcal{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(\frac{\beta }{N}\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+{\xi }_1\varphi {\eta }_C+{\xi }_2\rho {\eta }_A)\frac{\beta IS}{N}+(-{\xi }_1{\xi }_2{\xi }_3+{\xi }_1\omega \varphi +\alpha {\xi }_2\rho )I\\ & +{\eta }_C({\xi }_1{\xi }_2+{\xi }_1\varphi {\eta }_C+{\xi }_2\rho {\eta }_A)\frac{\beta CS}{N}+{\eta }_C(-{\xi }_1{\xi }_3{\xi }_2+{\xi }_1\varphi \omega +\rho \alpha {\xi }_2)C\\ & +{\eta }_A({\xi }_1{\xi }_2+{\xi }_1\varphi {\eta }_C+{\xi }_2\rho {\eta }_A)\frac{\beta AS}{N}+{\eta }_A(-{\xi }_2{\xi }_3{\xi }_1+\varphi \omega {\xi }_1+{\xi }_2\rho \alpha )A\\ & =\mathcal{D}\left(R_0\frac{S}{N}-1\right)I+{\eta }_C\mathcal{D}\left(R_0\frac{S}{N}-1\right)C+{\eta }_A\mathcal{D}\left(R_0\frac{S}{N}-1\right)A\\ & \le \mathcal{D}(R_0-1)I+{\eta }_C\mathcal{D}(R_0-1)C+{\eta }_A\mathcal{D}(R_0-1)A(\text{because}S\le N\text{ in }\Omega )\\ & \le 0\text{for}{R}_0<1.\end{array}$$

$$\begin{array}{rl}{S}^{\ast }& =\frac{\Lambda {\xi }_4\left({\xi }_1(\varphi +{\xi }_2)+\rho {\xi }_2\right)}{\mathcal{F}},I^{\ast }=\frac{-{\xi }_4\Lambda {\xi }_1{\xi }_2(\mathcal{D}-\mathcal{N})}{\mathcal{D}\mathcal{F}},\\ C^{\ast }& =\frac{-{\xi }_4\Lambda \varphi {\xi }_1(\mathcal{D}-\mathcal{N})}{\mathcal{D}\mathcal{F}},\\ A^{\ast }& =\frac{-{\xi }_4\Lambda \rho {\xi }_2(\mathcal{D}-\mathcal{N})}{\mathcal{D}\mathcal{F}},E^{\ast }=\frac{\psi \Lambda \left({\xi }_1(\varphi +{\xi }_2)+\rho {\xi }_2\right)}{\mathcal{F}},\end{array}$$

$$\lambda ={\beta }_1\left(I+{\eta }_{C}C+{\eta }_{A}A\right),\text{where}{\beta }_1=\frac{\beta \mu }{\Lambda }.$$

Full text

∗Corresponding author: Cristiana J. Silva (cjoaosilva@ua.pt)

Cristiana J. Silva∗ and Delfim F. M. Torres

Center for Research & Development in Mathematics and Applications (CIDMA)

Department of Mathematics, University of Aveiro

3810–193 Aveiro, Portugal

1. Introduction

The human immunodeficiency virus (HIV) is a retrovirus that causes HIV infection and, over time, acquired immunodeficiency syndrome (AIDS) [35]. The most significant advance in medical management of HIV infection has been the treatment of patients with antiviral drugs, which can suppress HIV replication to undetectable levels. Before 1996, few antiretroviral (ART) treatment options for HIV infection existed. The treatment of HIV infection was revolutionized in the mid 1990s with the introduction of drug regimens that combine inhibitors of the reverse transcriptase and protease and two of three essential enzymes of HIV. Combination of antiretroviral drugs, dramatically suppresses viral replication and reduces the HIV viral load to levels below the limits of detection of the most sensitive clinical assays, resulting in a significant reconstitution of the immune system [4].

The Global AIDS Update 2016 of the Joint United Nations Programme on HIV/AIDS, reports that the global coverage of ART therapy reached approximately 46% at the end of 2015. The gains in treatment are largely responsible for a 26% decline in AIDS-related deaths globally since 2010, from an estimated 1.5 million in 2010, to 1.1 million in 2015. Despite this significant achievement, there has not been a decrease in new infections since 2010, with more than 2 million new infections reported in 2015 [33]. The World Health Organization’s (WHO) Global Health Sector Strategy on HIV embraces innovation in the HIV response, recommending, for example, that people at substantial risk of HIV infection should be offered pre-exposure prophylaxis (PrEP) as an additional prevention choice, as part of comprehensive prevention. PrEP is the use of an antiretroviral medication to prevent the acquisition of HIV infection by uninfected individuals. Several trials among men who have sex with men, people who inject drugs, transgender people, women and serodiscordant couples (one partner is HIV-positive and the other is HIV-negative) have shown that when PrEP is taken, it is an effective and safe mechanism for preventing HIV-infection [36].

In [1], it is concluded that PrEP could prevent 2.7 to 3.2 million new cases of HIV in sub-Saharan Africa over 10 years, if it is targeted to the highest risk groups, and disinhibition could be prevented. In 2008, the authors of [14] claimed that PrEP represents the most powerful available biologic intervention for HIV prevention. In 2016, WHO has welcomed a plan by the South African Ministry of Health to provide immediate antiretroviral treatment to all sex workers with HIV, and to offer daily oral PrEP to HIV-negative sex workers to prevent them from acquiring the infection [41]. Following WHO, making PrEP drugs available for safe, effective prevention outside the clinical trial setting is the current challenge.

Substantial gaps remain in understanding the trade-offs between the costs and benefits of choosing alternative HIV prevention strategies, such as the initiation of PrEP by high risk uninfected individuals [13]. Mathematical models of HIV that include PrEP are scarce [6]. In [17], a mathematical model is used to estimate the effects of early diagnosis, early treatment and PrEP, on the HIV epidemic in South Korea over the next 40 years, as compared with the current situation. In [13], the authors develop a mathematical model to simulate HIV incidence among men residing in Los Angeles County, CA, aged between 15 to 65 years old, who have sex with men. They claim that PrEP and Test-and-Treat yield the largest reductions in HIV incidence, and are highly cost-effective. Another cost-effectiveness study was done in [3] and the results showed that PrEP can be cost-saving if delivered to individuals at increased risk of infection. In this paper, we propose a new mathematical epidemiological model for HIV/AIDS transmission including PrEP, which generalizes the HIV/AIDS sub-model recently proposed in [30].

First, we consider the mathematical model with no PrEP, calibrate the model to the cumulative cases of infection by HIV and AIDS from 1987 to 2014 reported in [25], and we show that the model predicts well the reality given in [25]. In this model, the effective contact with people infected with HIV includes two modification parameters that account for the relative infectiousness of individuals with AIDS symptoms, in comparison to those infected with HIV with no AIDS symptoms [37], and for partial restoration of the immune function of individuals with HIV infection that use ART correctly [10]. It should be noted that the infectiousness of the HIV-infected individuals under ART treatment is a controversial subject. HIV treatment reduces the viral load in the blood, and also in other body fluids, such as semen and vaginal fluid. However, not all people living with HIV who take HIV treatment and have an undetectable viral load in the blood also have an undetectable viral load in their other bodily fluids: see, e.g., [22] and references therein. In our paper, the values considered for the modification parameters are based on two different research studies: the first one is known as HPTN 052, where it was found that the risk of HIV transmission among heterosexual serodiscordant couples is 96% lower when the HIV-positive partner is on treatment [9]; and the other one, where it was proved that HIV-infected individuals under ART treatment have a very low probability (assumed inferior to $0.04$) of transmitting HIV [11]. We prove the global stability of the endemic equilibrium and we provide numerical simulations, illustrating the calibration of the model to the HIV/AIDS situation in Cape Verde.

In [23], the authors evaluate the effect of early HIV treatment and optimization of care, HIV testing, condom distribution, and substance abuse treatment on HIV incidence from 2011 to 2021, using Cape Verde as an example. However, they did not include PrEP for HIV negative groups at risk as a possible prevention measure. Here, we show that inclusion of PrEP can reduce significantly the HIV incidence. We start by proving existence and global stability of the disease free equilibrium of the HIV/AIDS-PrEP model. Moreover, we also prove existence of an unique endemic equilibrium and the global stability for some specific relevant cases. It is important, however, to highlight that PrEP is not for everyone [40]. Only people who are HIV-negative and at very high risk for HIV infection should take PrEP. Moreover, PrEP is highly expensive and it is still not approved in many countries, e.g., by the European Medicines Agency (EMEA) [32]. Therefore, the number of individuals that should take PrEP is limited at each instant of time. In order to study this health public problem, we formulate an optimal control problem with a mixed state control constraint.

Optimal control theory has been successfully applied to several epidemiological models, e.g., dengue [27, 28], tuberculosis [12, 29], Ebola [16, 24], cholera [20], and HIV/AIDS [26, 31]. However, we claim our work to be the first to apply optimal control to an HIV/AIDS model with PrEP. More precisely, we consider the HIV/AIDS-PrEP model and formulate an optimal control problem with the aim to determine the PrEP strategy that satisfies the mixed state control constraint and minimizes the number of individuals with pre-AIDS HIV-infection as well as the costs associated with PrEP. We solve the optimal control problem and provide numerical simulations, which show that it is possible to reduce the HIV incidence through PrEP and having into consideration the limitations related to the implementation of PrEP (cost, epidemic context, program coverage and individual-level adherence).

The paper is organized as follows. In Section 2, we consider the SICA model for HIV/AIDS transmission proposed in [30] and prove the global stability of the unique endemic equilibrium for the case when the associated AIDS-induced mortality is negligible. We calibrate the SICA model to the cumulative cases of infection by HIV and AIDS from 1987 to 2014 in Cape Verde and show that it predicts well this reality. In Section 3, we generalize the SICA model by including the possibility of providing PrEP to susceptible individuals. We prove the existence and global stability of the disease-free equilibrium for $R_{0}<1$, where $R_{0}$ denotes the basic reproduction number, which is computed following [34]. We prove existence of an unique endemic equilibrium point for $R_{0}>1$ and its global stability for a negligible AIDS-induced death rate and strict adherence to PrEP. Through numerical simulations, we investigate the impact of PrEP in the reduction of HIV transmission. In Section 4, we propose and analyze an optimal control problem with a mixed state control constraint. Numerical simulations show that the extremal solutions combine a reduction of HIV transmission with limited number of individuals under PrEP at each instant of time. We end with Section 5 of conclusions and future work.

2. The SICA model for HIV/AIDS transmission

In this section, we analyze a mathematical model for HIV/AIDS transmission with varying population size in a homogeneously mixing population, first proposed in [30], and prove the global stability of the unique endemic equilibrium for the case when the associated AIDS-induced mortality is negligible. 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 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 $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 [37]. On the other hand, $\eta_{C}\leq 1$ translates the partial restoration of immune function of individuals with HIV infection that use ART correctly [10]. 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 $\alpha$. 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$. These assumptions are translated in the following mathematical model:

[TABLE]

We consider the biologically feasible region

[TABLE]

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

Theorem 2.1** (See [30]).**

The population $N(t)$ is uniformly persistent, that is,

[TABLE]

with $\varepsilon>0$ not depending on the initial data.

Model (1) has a disease-free equilibrium, given by

[TABLE]

Following [34], the basic reproduction number $R_{0}$ for model (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$ and $\xi_{3}=\rho+\phi+\mu$.

Lemma 2.2** (See [30]).**

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

Theorem 2.3** (See [30]).**

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

To find conditions for the existence of an equilibrium for which HIV is endemic in the population (i.e., at least one of $I^{*}$, $C^{*}$ or $A^{*}$ is nonzero), denoted by $\Sigma_{+}=\left(S^{*},I^{*},C^{*},A^{*}\right)$, the equations in (1) are solved in terms of the force of infection at the steady-state $\lambda^{*}$, given by

[TABLE]

Setting the right hand side of the equations of the model to zero, and noting that $\lambda=\lambda^{*}$ at equilibrium gives

[TABLE]

with $D=-(\lambda^{*}+\mu)(\mu\left(\xi_{2}(\rho+\xi_{1})+\xi_{1}\phi+\rho d\right)+\rho\omega d)$, we use (5) in the expression for $\lambda^{*}$ in (4) to show that the nonzero (endemic) equilibrium of the model satisfies

[TABLE]

The force of infection at the steady-state $\lambda^{*}$ is positive only if $R_{0}>1$. Thus, the existence and uniqueness of the endemic equilibrium follows.

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

The model (1) has a unique endemic equilibrium whenever $R_{0}>1$.

Remark 1**.**

The expression of the endemic equilibrium of model (1) is given by

[TABLE]

Theorem 2.5** (See [30]).**

The endemic equilibrium $\Sigma_{+}$ is locally asymptotically stable for $R_{0}$ near 1.

2.1. Global stability of the endemic equilibrium for negligible AIDS-induced death rate ($d=0$)

In this section, we investigate the global stability of the endemic equilibrium of model (1) for the case when the associated AIDS-induced mortality is negligible ($d=0$). Adding the equations of the model (1), with $d=0$, gives $\dot{N}=\Lambda-\mu N$, so that $N\to\frac{\Lambda}{\mu}$ as $t\to\infty$. Thus, $\frac{\Lambda}{\mu}$ is an upper bound of $N(t)$ provided that $N(0)\leq\frac{\Lambda}{\mu}$. Further, if $N(0)>\frac{\Lambda}{\mu}$, then $N(t)$ decreases to this level. Using $N=\frac{\Lambda}{\mu}$ in the force of infection $\lambda=\frac{\beta}{N}\left(I+\eta_{C}\,C+\eta_{A}A\right)$ gives a limiting (mass action) system (see, e.g., [2]). Then, the force of infection becomes

[TABLE]

Therefore, we consider the model

[TABLE]

where $\xi_{1}=\alpha+\mu$. The model (6) has a unique endemic equilibrium given by $\tilde{\Sigma}_{+}=\Sigma_{+}|_{d=0}$, whenever $\tilde{R}_{0}=R_{0}|_{d=0}>1$. Let us define

[TABLE]

Theorem 2.6**.**

The endemic equilibrium $\tilde{\Sigma}_{+}$ of model (6) is globally asymptotically stable in $\Omega\backslash\Omega_{0}$ whenever $\tilde{R}_{0}>1$.

Proof.

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 (1) that

[TABLE]

Using relation $\Lambda=\beta_{1}\left(I^{*}+\eta_{C}\,C^{*}+\eta_{A}A^{*}\right)S^{*}+\mu S^{*}$, we have from the first equation of system (6) at steady-state that (8) 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 a well-known inequality: the geometric mean is less than or equal to the arithmetic mean. The equality $\frac{dV}{dt}=0$ holds if and only if $(S,I,C,A)$ take the equilibrium values $(S^{*},I^{*},C^{*},A^{*})$. Therefore, by LaSalle’s Invariance Principle [19], the endemic equilibrium $\Sigma_{+}$ is globally asymptotically stable. ∎

We conjecture that the endemic equilibrium for positive AIDS-induced death rate ($d>0$) is globally asymptotically stable. This remains an open question.

2.2. Case study for Cape Verde

In this section, we calibrate model (1) to the cumulative cases of infection by HIV and AIDS in Cape Verde from 1987 to 2014. We show that our model predicts well this reality. In Table 1, the cumulative cases of infection by HIV and AIDS in Cape Verde are depicted for the years 1987–2014 [25].

We consider the initial conditions (9) based on [25, 38]:

[TABLE]

We borrow the parameter values $\phi=1$, $\rho=0.1$, $\alpha=0.33$ and $\omega=0.09$ from [30]. Following the World Bank data [38, 42], the natural death rate is assumed to take the value $\mu=1/69.54$. The recruitment rate $\Lambda=13045$ was estimated in order to approximate the values of the total population of Cape Verde given in Table 1. See Figure 1, were we can observe that model (1) fits well the total population of Cape Verde.

The AIDS induced death rate is assume to be $d=1$ based on [39]. Two cases are considered: $\eta_{C}=0.04$, based on a research study known as HPTN 052, where it is found that the risk of HIV transmission among heterosexual serodiscordant couples is 96% lower when the HIV-positive partner is on treatment [9]; and $\eta_{C}=0.015$, which means that HIV-infected individuals under ART treatment have a very low probability of transmitting HIV, based on [11]. For the modification parameter $\eta_{A}\geq 1$ that accounts for the relative infectiousness of individuals with AIDS symptoms, in comparison to those infected with HIV with no AIDS symptoms, we assume $\eta_{A}=1.3$ and $\eta_{A}=1.35$, based in [37]. We estimated the value of the HIV transmission rate $\beta$ for $(\eta_{C},\eta_{A})=(0.04,1.35)$ equal to $0.695$ and for $(\eta_{C},\eta_{A})=(0.015,1.3)$ equal to $0.752$, and show that the model (1) predicts well the reality of Cape Verde for these parameter values: see Figure 2. All the considered parameter values are resumed in Table 2.

For the triplets $(\beta,\eta_{C},\eta_{A})=(0.752,0.015,1.3)$ and $(\beta,\eta_{C},\eta_{A})=(0.695,0.04,1.35)$, and the other parameter values from Table 2, we have that the basic reproduction number is given by $R_{0}=4.0983$ and $R_{0}=4.5304$, respectively.

3. The SICAE model

Now we generalize the model proposed in Section 2 by adding the possibility of providing PrEP to susceptible individuals. We add a class of individuals to the total population $N$, denoted by $E$, which represents the individuals that are under PrEP. The proportion of susceptible individuals that takes PrEP is denoted by $\psi$. We assume that PrEP is effective so that all susceptible individuals under PrEP treatment are transferred to class $E$. The individuals that stop PrEP become susceptible individuals again, at a rate $\theta$. Individuals under PrEP may suffer of natural death at a rate $\mu$. The model is given by the following system of ordinary differential equations:

[TABLE]

We consider the biologically feasible region

[TABLE]

3.1. Existence and stability of the disease-free equilibrium

Model (10) has a disease-free equilibrium, given by

[TABLE]

The linear stability of $\Sigma_{0}$ can be obtained using the next-generation method on system (10). Following [34], the basic reproduction number for model (10) is given by (3), that is,

[TABLE]

where $\xi_{1}=\alpha+\mu+d$, $\xi_{2}=\omega+\mu$, and $\xi_{3}=\rho+\phi+\mu$. Thus, from Theorem 2 of [34], the following result is established.

Lemma 3.1**.**

The disease free equilibrium $\Sigma_{0}$ of model (10), given by (11), is locally asymptotically stable if $R_{0}<1$, and unstable if $R_{0}>1$.

Biologically speaking, Lemma 3.1 implies that HIV infection can be eliminated from the community (when $R_{0}<1$) if the initial size of the population is in the basin of attraction of $\Sigma_{0}$. To ensure that elimination of HIV infection is independent of the initial size of the population, it is necessary to show that the disease free equilibrium is globally asymptotically stable [21]. This is obtained in what follows.

Theorem 3.2**.**

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

Proof.

Consider the following Lyapunov function:

[TABLE]

The time derivative of $V$ computed along the solutions of (1) is given by

[TABLE]

After some simplifications, we have

[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 the equations for $S$ and $E$ in system (10) shows, respectively, that $S\to S^{0}$ and $E\to E^{0}$ as $t\to\infty$. Thus, it follows from LaSalle’s Invariance Principle [19] that every solution of system (10) with initial conditions in $\Omega$ approaches the disease free equilibrium $\Sigma_{0}$ as $t\to\infty$ whenever $R_{0}<1$. ∎

The epidemiological significance of Theorem 3.2 is that HIV infection will be eliminated from the population if the threshold quantity, $R_{0}$, can be brought to a value less than unity.

3.2. Existence and stability of the endemic equilibrium

The unique endemic equilibrium of model (10) exists whenever $R_{0}>1$ and is given by

[TABLE]

where $\mathcal{F}=(\mathcal{N}-\rho d\xi_{2})\theta+(\mathcal{D}-\rho d\xi_{2})\psi-\mu(\rho d\xi_{2}-\mathcal{N})$ and $\xi_{4}=\theta+\mu$. We investigate the global stability of the endemic equilibrium of model (10) for the case when the associated AIDS-induced mortality is negligible ($d=0$) and there is a strict adherence to PrEP, that is, $\theta=0$. Adding the equations of model (10) with $d=0$ and $\theta=0$ gives $\dot{N}=\Lambda-\mu N$, so that $N\to\frac{\Lambda}{\mu}$ as $t\to\infty$. Thus, $\frac{\Lambda}{\mu}$ is an upper bound of $N(t)$, provided $N(0)\leq\frac{\Lambda}{\mu}$. Further, if $N(0)>\frac{\Lambda}{\mu}$, then $N(t)$ decreases to this level. Using $N=\frac{\Lambda}{\mu}$ in the force of infection $\lambda=\frac{\beta}{N}\left(I+\eta_{C}\,C+\eta_{A}A\right)$ gives a limiting (mass action) system. Then, the force of infection becomes

[TABLE]

Therefore, we consider the following model:

[TABLE]

For system (12), the basic reproduction number is given by $R_{0}=\frac{\Lambda\mathcal{N}_{1}}{(\mu+\psi)\mathcal{D}}$ with $\mathcal{N}_{1}=\beta_{1}\,\left(\xi_{2}\left(\xi_{1}+\rho\,\eta_{A}\right)+\eta_{C}\,\phi\,\xi_{1}\right)$. When $R_{0}>1$, system (12) has a unique endemic equilibrium $\tilde{\Sigma}_{+}=(\tilde{S},\tilde{I},\tilde{C},\tilde{A},\tilde{E})$ given by

[TABLE]

where $\mathcal{D}_{1}=\beta_{1}\mu(\rho\xi_{2}+\xi_{1}(\phi+\xi_{2}))(\xi_{1}(\xi_{2}+\eta_{C}\phi)+\eta_{A}\rho\xi_{2})$. Since $E$ does not appear in the first four equations of system (12), we only need to consider the reduced system

[TABLE]

Whenever $R_{0}>1$, system (14) has a unique endemic equilibrium $\Sigma_{1}=(\tilde{S},\tilde{I},\tilde{C},\tilde{A})$ with $\tilde{S}$, $\tilde{I}$, $\tilde{C}$ and $\tilde{A}$ defined by (13). Consider the Lyapunov function (7) given by

[TABLE]

Differentiating $V$ with respect to time gives

[TABLE]

Substituting the expressions for the derivatives in $\dot{V}$, it follows from (12) that

[TABLE]

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

[TABLE]

which can be simplified to

[TABLE]

Using the relations

[TABLE]

at the steady state, after some simplifications we have

[TABLE]

The terms between the larger brackets are less than or equal to zero because the geometric mean is less than or equal to the arithmetic mean. Equality $\frac{dV}{dt}=0$ holds if and only if $S=S^{*}$, $I=I^{*}$, $C=C^{*}$ and $A=A^{*}$. By LaSalle’s Invariance Principle [19], every solution to the equations in model (14) with initial conditions in $\left\{(S,I,C,A)\in\mathbb{R}_{+0}^{4}\,:\,S\leq\frac{\Lambda}{\psi+\mu},N\leq\frac{\Lambda}{\mu}\right\}\backslash\left\{S=I=C=A\right\}$ approaches the endemic equilibrium $\tilde{\Sigma}_{+}$, which implies that the endemic equilibrium $\tilde{\Sigma}_{+}$ of system (13) is globally asymptotically stable on $\Omega_{P}\backslash\Omega_{P0}$, where $\Omega_{P0}=\{\left(S,I,C,A\right)\in\Omega_{P}\,:\,I=C=A=0\}$. We have just proved the following result.

Theorem 3.3**.**

The endemic equilibrium $\tilde{\Sigma}_{+}$ of model (12) is globally asymptotically stable in $\Omega_{P}\backslash\Omega_{P0}$ whenever $R_{0}>1$.

A generalization of Theorem 3.3 to the case $d>0$ and $\theta>0$ remains an open question.

3.3. Numerical simulations

In this section, we investigate the impact of PrEP in the reduction of HIV transmission. We assume that the total population is constant, that is, $\Lambda=\mu N$ with $\mu=1/69.54$ and $d=0$. The initial conditions are given by

[TABLE]

We start our numerical simulations assuming that 10 per cent of the susceptible individuals take PrEP, that is, $\psi=0.1$ and the default rate takes the value $\theta=0.001$. Consider $\beta=0.582$, $\eta_{C}=0.04$, $\eta_{A}=1.35$, and that the parameters $\omega$, $\rho$, $\phi$ and $\alpha$ take the values of Table 2. We compare the case $(\psi,\theta)=(0.1,0.001)$ with $(\psi,\theta)=(0,0)$ for $t\in[0,t_{f}]$, $t_{f}=25$ years, where $(\psi,\theta)=(0,0)$ means that no susceptible individual was under PrEP.

From Figure 3, we observe that PrEP reduces the number of individuals with HIV infection. In fact, for $(\psi,\theta)=(0.1,0.001)$, we have $S(25)\simeq 1687$, $I(25)\simeq 67$, $C(25)\simeq 626$, $A(25)\simeq 20$ and $E(25)=7800$. It is important to note that $I(t)$ is a decreasing function for all $t\in[0,25]$ and that the maximum value for the number of HIV-infected individuals with AIDS symptoms is less than 22 for $t\in[0,25]$. On the other hand, at the end of the 25 years, the number of individuals that are taking PrEP is equal to 7800, which is highly expensive (the PrEP drug costs between $8,000 and$14,000 per year for each individual). Therefore, it is of most importance to establish what is the optimal proportion of susceptible individuals that should take PrEP, taking into consideration its costs. In Section 4, we formulate this problem mathematically, using the theory of optimal control, and we study it both analytically and numerically.

4. Optimal control problem with a mixed state control constraint

Substantial gaps remain in understanding the trade-offs between costs and benefits of choosing alternative HIV prevention strategies, such as the initiation of PrEP by high risk uninfected individuals [13]. Following WHO, making PrEP drugs available for safe, effective prevention outside the clinical trial setting is the current challenge. However, it is important to highlight and recall that PrEP is not for everyone: only people who are HIV-negative and at very high risk for HIV infection should take PrEP [40]. Moreover, PrEP is highly expensive and it is still not approved in many countries like, for example, by the European Medicines Agency (EMEA) [32]. Therefore, the number of individuals that should take PrEP is limited at each instant of time for a fixed interval of time $[0,t_{f}]$. In order to study this health public problem, from an optimal point of view, we formulate an optimal control problem with a mixed state control constraint. For the usefulness of such problems in epidemiology see, e.g., [5]. We consider the model with PrEP (10) and formulate an optimal control problem with the aim to determine the PrEP strategy $\psi$ that minimizes the number of individuals with pre-AIDS HIV-infection $I$ as well as the costs associated with PrEP. We assume that the fraction of individuals that takes PrEP, at each instant of time, is a control function, that is, $\psi\equiv u(t)$ with $t\in[0,t_{f}]$, and that the total population $N$ is constant: the recruitment rate is proportional to the natural death rate, $\Lambda=\mu N$, and there are no AIDS-induced deaths ($d=0$). Precisely, we consider the model with control $u(t)$ given by

[TABLE]

and formulate an optimal control problem with the aim to determine the PrEP strategy $u$ over a fixed interval of time $[0,t_{f}]$ that minimizes the cost functional

[TABLE]

where the constants $w_{1}$ and $w_{2}$ represent the weights associated with the number of HIV infected individuals $I$ and on the cost associated with the PrEP prevention treatment, respectively. It is assumed that the control function $u$ takes values between 0 and 1. When $u(t)=0$, no susceptible individual takes PrEP at time $t$; if $u(t)=1$, then all susceptible individuals are taking PrEP at time $t$. Let $\vartheta$ denote the total number of susceptible individuals under PrEP for a fixed time interval $[0,t_{f}]$. This constraint is represented by

[TABLE]

which should be satisfied at almost every instant of time during the whole PrEP program. Let

[TABLE]

The optimal control problem consists to find the optimal trajectory $\tilde{x}$, associated with the control $\tilde{u}$, satisfying the control system (17), the initial conditions (16),

[TABLE]

the constraint (19), and where the control $\tilde{u}\in\Omega$ minimizes the objective functional (18) with

[TABLE]

The control system can be rewritten in the following way:

[TABLE]

with

[TABLE]

and

[TABLE]

where $b=\left(-1\,0\,0\,0\,1\right)^{T}$ and $Z=0$ with [math] the $5\times 4$ null matrix and $f=\left(f_{1}\,f_{2}\,0\,0\,0\right)$ with

[TABLE]

and

[TABLE]

It follows from Theorem 23.11 in [7] that problem (17)–(20) has a solution (see also [5]). Let $(\tilde{x},\tilde{u})$ denote such solution. To determine it, we apply the Pontryagin Maximum Principle (see, e.g., Theorem 7.1 in [8]): there exist multipliers $\lambda_{0}\leq 0$, $\lambda\in AC([0,t_{f}];{\mathbb{R}}^{5})$, and $\nu\in L^{1}([0,t_{f}];{\mathbb{R}})$, such that

• •

$\min\{|\lambda(t)|\,:\,t\in[0,t_{f}]\}>\lambda_{0}$ (nontriviality condition);

• •

$\frac{d\lambda(t)}{dt}=-\frac{\partial\mathcal{H}}{\partial x}(\tilde{x}(t),\tilde{u}(t),\lambda_{0},\lambda(t),\nu(t))$ (adjoint system);

• •

$\lambda(t)B\tilde{x}(t)+\nu(t)\tilde{x}_{1}(t)+\lambda_{0}w_{2}\tilde{u}^{2}(t)\in\mathcal{N}_{[0,1]}(\tilde{u}(t))$ a.e. and

[TABLE]

(minimality condition);

• •

$\nu(t)(\tilde{x}_{1}(t)\tilde{u}(t)-\vartheta)=0$ and $\nu(t)\leq 0$ a.e.

• •

$\lambda(t_{f})=(0,\ldots,0)$ (transversality conditions);

where the Hamiltonian $\mathcal{H}$ for problem (17)–(20) is defined by

[TABLE]

and $\mathcal{N}_{[0,1]}(\tilde{u}(t))$ stands for the normal cone from convex analysis to $[0,1]$ at the optimal control $\tilde{u}(t)$ (see, e.g., [7]). The optimal solution $(\tilde{x},\tilde{u})$ is normal (see [5] for details), so we can choose $\lambda_{0}=1$. The unique optimal control $\tilde{u}$ is given by

[TABLE]

where the adjoint functions satisfy

[TABLE]

We solve the optimal control problem (17)–(20) numerically for concrete parameter values and initial conditions. The initial conditions are given by (16) and the HIV transmission parameters take the values $\beta=0.582$, $\eta_{C}=0.04$, and $\eta_{A}=1.35$. We assume that the default PrEP rate is equal to $\theta=0.001$ and the parameters $\omega$, $\rho$, $\phi$, $\alpha$ take the values given by Table 2. The weight constants take the values $w_{1}=w_{2}=1$. We start by solving the optimal control problem without the mixed state control constraint and compare the extremals with the case where the fraction of individuals under PrEP is constant for all $t\in[0,25]$ and equal to $\psi=0.1$ and $\psi=0.9$.

In Figure 4, we observe that the number of HIV-infected individuals associated with the optimal control solution $\tilde{I}$, $\tilde{C}$, $\tilde{A}$ is lower than the respective number associated with the constant values $\psi=0.1$ or $\psi=0.9$. On the other hand, we observe that the optimal control starts taking the maximum value $1$, which means that all susceptible individuals should be under PrEP (see Figure 5a) but this situation does not respect the mixed state control constraint (see Figure 5b) and implies higher implementation costs.

Now, consider the mixed state-control constraint

[TABLE]

We start by comparing the optimal control solutions $\tilde{I}$, $\tilde{C}$, $\tilde{A}$ and $\tilde{E}$ associated with the optimal control $\tilde{u}$ with the solution of model (10) with fixed values for the fraction of susceptible individuals under PrEP, $\psi=0.1,0.9$.

In Figure 6, we observe that optimal control solutions are not associated with the lowest values of the number of individuals with HIV infection. In fact, if we consider the case where 90 per cent of the susceptible population is under PrEP, the number of individuals with HIV infection is lower than the corresponding values associated with the optimal control solution. This is related to the mixed constraint (21). In Figure 7b, we observe that the mixed constraint (21) is satisfied for all $t\in[0,25]$. The optimal control $\tilde{u}$ starts with the value $0.2$ and is an increasing function for $t\in[0,4.05]$ years; for $t\in[4.05,6.88]$, the optimal control $\tilde{u}$ takes the maximum value one, which means that all susceptible individuals should be taking PrEP. For the rest of the simulation period, the optimal control $\tilde{u}$ is a decreasing function, taking values less than one. The mean value of the optimal control is approximately $0.61$ (see Figure 7a).

In Figure 8, we simulate the case where the fraction of susceptible individuals under PrEP is constantly equal to the mean value of the optimal control $\tilde{u}$, that is, $\psi=0.61$, and compare the solutions of model (10) with the solutions of the optimal control problem with mixed state control constraint. We observe that the number of HIV infected individuals is lower for $\psi=0.61$ than the ones associated with the optimal control solutions. However, the total number of individuals under PrEP at each instant of time is much bigger for $\psi=0.61$ (see Figure 7c). In fact, we have $\int_{0}^{25}\tilde{u}(t)\tilde{S}(t)dt=12553$ and $\int_{0}^{25}0.61\,S(t)dt=13201$, where $S(t)$ denotes the solution of model (10) for $\psi=0.61$. We conclude that the cost associated with the PrEP strategy $\psi=0.61$ is bigger than the one associated with the optimal control strategy. On the other hand, for the optimal control problem with no mixed state control constraint, we have $\int_{0}^{25}\tilde{u}(t)\tilde{S}(t)dt=12836$, which means that this solution has a bigger cost than the one that satisfies the mixed constraint.

Let us now suppose that we increase the value of $w_{2}$, that is, the weight associated with the cost of submitting susceptible individuals to PrEP. An increase of $w_{2}$ is associated with a decrease of the number of individuals that will be under PrEP. For example, for $w_{2}=10000$, the maximum value of the optimal control $\tilde{u}$ is approximately $0.17$ and $\int_{0}^{25}\tilde{u}(t)\tilde{S}(t)dt=6652.5$, that is, the cost associated with the PrEP strategy decreases significantly. However, the number of HIV-infected individuals increases, namely $\tilde{I}(25)\simeq 110$, $\tilde{C}(25)\simeq 813$ and $\tilde{A}(25)\simeq 29$.

5. Conclusion and future work

We have proved the global stability of the endemic equilibrium point of the SICA model for HIV/AIDS transmission proposed in [30], in the case where the AIDS-induced death rate is negligible. Two different values for the relative infectiousness of HIV-infected individuals under ART treatment were considered, based on the studies reported in [9, 11] ($\eta_{C}=0.04$ and $\eta_{C}=0.015$) and two values for the relative infectiousness of HIV-infected individuals with AIDS symptoms ($\eta_{A}=1.3$ and $\eta_{A}=1.35$), based in [37]. We have shown that taking the parameter values given in Table 2, the SICA model describes well the cumulative cases of infection by HIV and AIDS in Cape Verde for the period from 1987 to 2014 [25, 42]. Furthermore, we generalized the SICA model to a HIV/AIDS-PrEP model (SICAE model) including PrEP as an HIV prevention strategy. We proved the existence and uniqueness of disease-free and endemic equilibrium points, depending on the value of the basic reproduction number $R_{0}$. Through Lyapunov functions and LaSalle’s Invariance Principle, we proved the global stability of the disease-free equilibrium for $R_{0}<1$ and the global stability of the endemic equilibrium point when $R_{0}>1$ for negligible AIDS-induced death rate and strict adherence to PrEP. The number of HIV-infected individuals for the SICA and SICAE models was compared. We concluded that PrEP reduces the number of new HIV infections. However, only people who are HIV-negative and at very high risk for HIV infection should take PrEP. Therefore, the number of individuals that should take PrEP must be limited at each instant of time, for a fixed interval of time. In order to study this health public problem, we formulated an optimal control problem with a mixed state control constraint, where the objective is to determine the optimal PrEP strategy that satisfies the mixed constraint (the total number of individuals under PrEP at each instant of time is limited) and minimizes the number of pre-AIDS HIV-infected individuals as well as the cost associated with the implementation of PrEP. The optimal control problems with and without the mixed state control constraint were solved and compared with the optimal solutions in the case where the fraction of individuals under PrEP is constant. Optimal control theory gave us PrEP strategies that minimize the number of HIV-infected individuals, the cost associated with PrEP and satisfy the limitations on the number of total individuals that should be under PrEP at each instant of time.

It remains open the questions of how to prove the global stability of the endemic equilibrium point of the SICA model for positive AIDS-induced death rate and of the SICAE model for positive AIDS-induced death and PrEP default rates. As for the optimal control problem with mixed state control constraint, we believe it will be also interesting to consider a $L^{1}$ cost functional and a variable total population size. These and other questions are under investigation and will be addressed elsewhere.

Acknowledgments

This research was partially supported by Portuguese Foundation for Science and Technology (FCT) through the R&D unit CIDMA, reference UID/MAT/04106/2013, 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 the FCT post-doc fellowship SFRH/BPD/72061/2010. The authors are grateful to two anonymous referees for several pertinent questions and suggestions.