Forecasting the Effect of Pre-Exposure Prophylaxis (PrEP) on HIV Propagation with a System of Differential–Difference Equations with Delay

Abstract

The HIV/AIDS epidemic is still active worldwide with no existing definitive cure. Based on the WHO recommendations stated in 2014, a treatment, called Pre-Exposure Prophylaxis (PrEP), has been used in the world, and more particularly in France since 2016, to prevent HIV infections. In this paper, we propose a new compartmental epidemiological model with a limited protection time offered by this new treatment. We describe the PrEP compartment with an age-structure hyperbolic equation and introduce a differential equation on the parameter that governs the PrEP starting process. This leads us to a nonlinear differential–difference system with discrete delay. After a local stability analysis, we prove the global behavior of the system. Finally, we illustrate the solutions with numerical simulations based on the data of the French Men who have Sex with Men (MSM) population. We show that the choice of a logistic time dynamics combined with our Hill-function-like model leads to a perfect data fit. These results enable us to forecast the evolution of the HIV epidemics in France if the populations keep using PrEP.

1. Introduction

Since its onset in the early 1980s and its clear identification in 1983, Human Immunodeficiency Virus (HIV) and then Acquired Immune Deficiency Syndrome (AIDS) have still comprised one of the most deadly active worldwide epidemics. In the 2019 UNAIDS study, about 38 million persons lived with HIV, $1.7$ million became infected, and 690,000 died of AIDS-related diseases (UNAIDS, Global HIV & AIDS statistics—2020 fact sheet: https://www.unaids.org/en/resources/fact-sheet) (accessed on 13 September 2022), becoming one of the most serious public health challenges.

It is well known now that this infection evolves in three stages: first a short acute phase where flu-like symptoms appear, followed by a symptom-free chronic phase that lasts between 10 and 15 years. It eventually ends up with AIDS when the virus has killed enough TCD4, leading to the failure of the immune system (https://www.hiv.gov/hiv-basics/overview/about-hiv-and-aids/symptoms-of-hiv) (accessed on 13 September 2022).

Note here that the latency stage appears as one of the crucial problems. Indeed, the dormant period of the virus drastically delays HIV/AIDS diagnosis if not detected and plays a major role in the epidemic’s spread (see Figure 1 in [1]).

Despite extensive investigations, there is still no existing therapy that helps the organism to fully get rid of the virus.

However, since 1996, AntiRetroviral Therapy (ART), consisting of a combination of three (or more) antiretroviral agents taken daily, has appeared to reduce significantly the presence of the virus under a detectable threshold. Currently, not only efficient at avoiding the fatality of the disease, ART allows also HIV-infected individuals stop spreading the disease. In 2019, $25.4$ million of them successfully accessed antiretroviral therapy, (UNAIDS, Global HIV & AIDS statistics—2020 fact sheet: https://www.unaids.org/en/resources/fact-sheet (accessed on 13 September 2022); World Health Organization (WHO)—HIV/AIDS: https://www.who.int/news-room/fact-sheets/detail/hiv-aids (accessed on 13 September 2022).

In the late 2010s, Pre-Exposure Prophylaxis (PrEP) was introduced as a new way of preventing HIV transmission. This successful treatment addresses HIV-free individuals showing a high risk of becoming infected. PrEP consists of taking a combination of two antiretrovirals in two different protocols: either on a daily basis (continuous treatment) or on demand (discrete treatment), which is just at least 2 h before and 2 days after sexual intercourse. Every 3 months, continuous PrEP users need to have a global Sexually Transmitted Infections (STIs) screening in order to obtain a prescription for another 3 months (World Health Organization (WHO)—HIV/AIDS: https://www.who.int/news-room/fact-sheets/detail/hiv-aids; https://www.who.int/teams/global-hiv-hepatitis-and-stis-programmes/hiv/prevention/pre-exposure-prophylaxis) (accessed on 13 September 2022). Theoretically, they can stop anytime, but they are practically more likely to choose at the end of each quarter. We remind here that PrEP is preventive and does not cure HIV. However, it effectively reduces HIV transmission (World Health Organization (WHO)—HIV/AIDS: https://www.who.int/news-room/fact-sheets/detail/hiv-aids) (accessed on 13 September 2022).

Our objective here was to design a new model to mathematically forecast the influence of PrEP users among high-risk profiles on the epidemic. The application of our study focused only on the French population, where PrEP has been used since 2016, and is currently followed by an average of 20,000 individuals [2].

In the past few decades, many epidemiological models have been used to describe HIV’s dynamics. Some of them used the standard Susceptible (S)–Infected (I)–Infected under ART (C) with low viral charge–Infected with AIDS symptoms (A) ($SICA$) model (see [3] for instance), and for a good review, we suggest [4,5].

It has only been in the past few years that new models including a PrEP user compartment appeared and started to receive our attention. A good example is the preliminary work [6], where the authors suggested a vaccination approach, which could be seen as a precursor to PrEP. Their numerical simulations showed that the epidemic spread could be controlled thanks to vaccination, even if the basic reproduction number remains above 1. In [7,8], an explicit PrEP compartment was included in the $SICA$ model. The parameters in [8] were adjusted with clinical data and PrEP’s effectiveness mathematically proven. In [9], the authors divided the PrEP compartment according to the adherence of users to the treatment. They showed that, with at least $70\%$ of PrEP users in the male homosexual population, the HIV epidemic could be effectively controlled. Finally, in [10], PrEP was combined with screening and the result was compared to Portuguese data showing that these two processes appeared necessary to control the HIV epidemic.

In our paper, we introduce a new approach with a Susceptible–Infected–Protected ($SIP$) model. Inspired by [8,11], it includes an age structure on the PrEP (protected) compartment, corresponding to the time spent after the onset of the new 3-month treatment period. Once the term is reached, the user chooses either to keep taking PrEP or to stop it, becoming susceptible again. On the other hand, individuals of the S compartment may decide to start or resume this preventative process and, hence, reach the P population.

Our goal was to propose a more accurate model than that of [8] that still possible to handle analytically by considering a more complex susceptible interaction $F(t,S)$ depending both on time t (through a function $\psi$) and the total S population (through a function f), which can be seen as a political or economic decision made depending on the size of the HIV-exposed population(see Section 2 for details). After standard modifications, we prove in Section 2 that our new model can be written as a nonlinear differential–difference system with discrete delay. As presented below in Section 3 and Section 4, the stability analysis is then possible when the $\psi$ term is constant. However, since it may take a while to convince a population set to start PrEP, it appears natural to us that $\psi$ may follow a logistic equation. The model becomes then too challenging to be investigated in an analytical way, and thus, numerical simulations take over this work in Section 5. Thanks to this new assumption, we show eventually that it is possible to fit the data with realistic values. Before this, we proceed to its well-posedness and stability study and investigate the constant and non-constant joining PrEP rate subcases in Section 3 and Section 4. Then, we compare our theoretical work to official clinical French data in order to understand the role of PrEP’s dynamics on the HIV/AIDS epidemic in Section 5. With some numerical simulations, we observe that, by choosing a logistic equation for the nonlinear joining PrEP rate and a Hill function for the dynamics of $f\left(S\right)$, we perfectly fit our data.

2. Introducing Our Model

In order to keep the model as simple as possible, we considered only three compartments: Susceptible (S), individuals that may be infected by HIV; Infected (I) and Protected (P), individuals on PrEP treatment (see Figure 1). The total population is

$\begin{array}{c}\hfill N\left(t\right)=S\left(t\right)+I\left(t\right)+P\left(t\right).\end{array}$

To take into account the limited protection duration offered by PrEP, we structured the compartment of protected individuals by age. Here, age a corresponds to the time since the last test and the renewal of PrEP treatment. The maximum duration of the protection period is generally three months, $\tau =3$ months. Therefore, $a\in [0,\tau ]$. If an individual decides to renew the treatment, the age a is reset to zero. We denote by $p(t,a)$ the population on PrEP at time t who started their new treatment period a days ago. Thus, the total population of PrEP users at time t is

$\begin{array}{c}\hfill P\left(t\right)={\int }_0^{\tau }p(t,a)da.\end{array}$

The model was designed taking into account the input and output rates of each compartment with appropriate parameters (see Figure 1). At each time t, there is a constant number of individuals (source term $\sigma$) who become susceptible by reaching the age of sexual consent, becoming single again, or deciding to start new intercourse experiences. At the same time, some susceptible individuals become infected with HIV at a rate $\beta$. A part,

$F(t,S(t\left)\right)=\psi \left(t\right)f\left(S\right(t\left)\right),$

of the susceptible individuals decide to start PrEP. The choice of functions f and $\psi$ allows us to consider many scenarios. The function f is an increasing function of the total population of susceptible individuals with $f\left(0\right)=0$. Examples could be

$f\left(S\right)=S^n\mathrm{or}f\left(S\right)={\displaystyle \frac{S^n}{1+S^n}},n\ge 1.$

The first choice is well adapted to rich countries or countries with a low number of susceptible individuals (the proposed treatment, for example the case $f\left(S\right)=S$, $n=1$, is proportional to the total population of susceptible individuals). The second choice is better adapted to poor countries or countries with a large number of susceptible individuals, which, for budgetary reasons, cannot increase the treatment indefinitely. The fact that the function $\psi$ is time-dependent is dictated by the data we wish to use on PrEP, which began in 2016 in France and continues to increase. For cost reasons, this PrEP prescription rate will eventually reach a maximum threshold. A well-fitting example of function $\psi$ is

$\psi \left(t\right)={\psi }_0+{\displaystyle \frac{{\kappa }_1{t}^m}{t^m+{\kappa }_2}},m\ge 1,{\kappa }_1,{\kappa }_2>0.$

For the analytical study of the model, we assumed that this threshold has already been reached and, therefore, considered

$\psi \equiv K:={\psi }_0+{\kappa }_1$

to be a constant. We can also choose $t\mapsto \psi \left(t\right)$ as a solution of the logistic equation with K as the carrying capacity. Finally, we supposed that the disease and treatment do not affect mortality (death rate $\mu$ is consequently the same for all). The system satisfied by S and I is then given, for all $t>0$, by

(1)$\left\{\begin{array}{ccc}{S}^{\prime }\left(t\right)\hfill & =& {\displaystyle \sigma -\beta I\left(t\right)S\left(t\right)-\mu S\left(t\right)-Kf\left(S\right(t\left)\right)+(1-\theta )p(t,\tau ),}\hfill \\ I^{\prime }\left(t\right)\hfill & =& {\displaystyle \beta I\left(t\right)S\left(t\right)-\mu I\left(t\right).}\hfill \end{array}\right.$

The density p satisfies, for $t>0$, the following age-structured partial differential equation:

(2)$\left\{\begin{array}{cc}{\displaystyle \frac{\partial p}{\partial t}}(t,a)+{\displaystyle \frac{\partial p}{\partial a}}(t,a)=-\mu p(t,a),\hfill & 0<a<\tau ,\hfill \\ p(t,0)=Kf\left(S\right(t\left)\right)+\theta p(t,\tau ).\hfill \end{array}\right.$

We remind here that $a\in [0,\tau ]$ is the time elapsed from the new period of time under treatment. If an individual decides to renew the treatment, the age a is reset to zero. We denote by $p(t,a)$ the population under PrEP at time t who started their new treatment period a days ago. The initial condition is given by

$\begin{array}{c}\hfill S\left(0\right)=S_0,I\left(0\right)=I_0\mathrm{and}p(0,a)=p_0\left(a\right),0<a<\tau .\end{array}$

The boundary condition

$p(t,0)=Kf\left(S\right(t\left)\right)+\theta p(t,\tau ),$

represents individuals who renewed their PrEP, $\theta p(t,\tau )$, or those who started it, $Kf\left(S\right(t\left)\right)$.

Using the characteristics method and the boundary condition from (2) (see [12]), we obtain, for $t>0$ and $a\in [0,\tau ]$,

(3)$p(t,a)=\left\{\begin{array}{cc}{e}^{-\mu a}p(t-a,0),\hfill & t>a,\hfill \\ e^{-\mu t}p(0,a-t)=e^{-\mu t}{p}_0(a-t),\hfill & 0\le t\le a.\hfill \end{array}\right.$

To write Equation (3) in a convenient form, we supposed that the initial condition $p_0$ is continuous on $[0,\tau ]$ and satisfies the compatibility condition

$p_0\left(0\right)=Kf\left(S_0\right)+\theta p_0\left(\tau \right).$

We put

$\phi \left(t\right):=e^{-\mu t}{p}_0(-t),-\tau \le t\le 0,$

and we define the function

$u\left(t\right):=p(t,0),t\ge 0.$

We then obtain from (3) the expression

$p(t,\tau )=e^{-\mu \tau }\left\{\begin{array}{cc}u(t-\tau ),\hfill & t>\tau ,\hfill \\ \phi (t-\tau ),\hfill & 0\le t\le \tau .\hfill \end{array}\right.$

This means that we can extend the function u to the interval $[-\tau ,0]$, by putting $u\left(t\right)=\phi \left(t\right)$ for $t\in [-\tau ,0]$. Then, we have directly $p(t,\tau )=e^{-\mu \tau }u(t-\tau )$, for all $t\ge 0$. Consequently, using the boundary condition of the system (2), u satisfies the difference equation:

$u\left(t\right)=\left\{\begin{array}{cc}Kf\left(S\left(t\right)\right)+\theta e^{-\mu \tau }u(t-\tau ),\hfill & t>0,\hfill \\ \phi \left(t\right),\hfill & -\tau \le t\le 0.\hfill \end{array}\right.$

Furthermore, the total population $P\left(t\right)$ can be explicitly expressed, for $t\ge 0$, in terms of the function u:

$P\left(t\right)={\int }_0^{\tau }{e}^{-\mu a}u(t-a)da=e^{-\mu t}{\int }_{t-\tau }^t{e}^{\mu a}u\left(a\right)da.$

Thus, we obtain the complete model, for $t>0$,

(4)$\left\{\begin{array}{ccc}{S}^{\prime }\left(t\right)\hfill & =& {\displaystyle \sigma -\beta I\left(t\right)S\left(t\right)-\mu S\left(t\right)-Kf\left(S\left(t\right)\right)+(1-\theta )e^{-\mu \tau }u(t-\tau ),}\hfill \\ I^{\prime }\left(t\right)\hfill & =& {\displaystyle \beta I\left(t\right)S\left(t\right)-\mu I\left(t\right),}\hfill \\ u\left(t\right)\hfill & =& Kf\left(S\left(t\right)\right)+\theta e^{-\mu \tau }u(t-\tau ),\hfill \end{array}\right.$

with the initial condition:

(5)$S\left(0\right)=S_0,I\left(0\right)=I_0\mathrm{and}u\left(t\right)=\phi \left(t\right),t\in [-\tau ,0].$

Knowing that ${\displaystyle P\left(t\right)=e^{-\mu t}{\int }_{t-\tau }^t{e}^{\mu a}u\left(a\right)da}$, we obtain

$P^{\prime }\left(t\right)=-\mu P\left(t\right)+Kf\left(S\left(t\right)\right)-(1-\theta )e^{-\mu \tau }u(t-\tau ).$

We then have

$N^{\prime }\left(t\right)=S^{\prime }\left(t\right)+I^{\prime }\left(t\right)+P^{\prime }\left(t\right)=\sigma -\mu N\left(t\right).$

We easily obtain

(6)$N\left(t\right)=\left(N_0-\frac{\sigma }{\mu }\right)e^{-\mu t}+{\displaystyle \frac{\sigma }{\mu }}\to {\displaystyle \frac{\sigma }{\mu }},\mathrm{as}t\to +\infty .$

3. Mathematical Analysis

In this section, we discuss the well-posedness of the model, the steady-states, the basic reproduction number, and the local asymptotic stability.

3.1. Well-Posedness of the Model

The problem (4) and (5) is a coupled system of nonlinear differential and difference equations with discrete delay. We directly derive the following proposition about the well-posedness of the model.

Proposition 1.

For any nonnegative initial condition $(S_0,I_0,\phi )$, where φ is a nonnegative continuous function on $[-\tau ,0]$, there exists a unique nonnegative solution to Problem (4) and (5) defined on $[0,+\infty )$. Moreover, this solution is uniformly bounded on $[0,+\infty )$.

Proof. 

We first solve the system (4) and (5) on the interval $[0,\tau ]$. In this case, since u is completely defined, the Cauchy–Lipschitz theorem gives the existence and uniqueness of the solution on $[0,\tau ]$, and this solution is nonnegative. Then, we repeat this method on each interval of type $[k\tau ,(k+1\left)\tau \right]$ with $k\in \mathbb{N}$. Thus, we obtain a unique nonnegative solution on the interval $[0,+\infty )$. Furthermore, using (6) and the fact that $S\left(t\right)+I\left(t\right)+u\left(t\right)\le S\left(t\right)+I\left(t\right)+P\left(t\right)=N\left(t\right)$, we deduce that $(S,I,u)$ is uniformly bounded on $[0,+\infty )$. □

3.2. Steady-States and Basic Reproduction Number

In this subsection, we compute the steady-states of the system and study their local asymptotic stability. Let us consider a steady-state $(S^{*},I^{*},u^{*})$ of (4). We then have

(7)$\begin{array}{c}\hfill \left\{\begin{array}{ccc}0\hfill & =& {\displaystyle \sigma -\beta I^{*}{S}^{*}-\mu S^{*}-Kf\left(S^{*}\right)+(1-\theta )e^{-\mu \tau }{u}^{*},}\hfill \\ 0\hfill & =& (\beta S^{*}-\mu )I^{*},\hfill \\ u^{*}\hfill & =& Kf\left(S^{*}\right)+\theta e^{-\mu \tau }{u}^{*}.\hfill \end{array}\right.\end{array}$

The third equation of the system (7) leads to

$u^{*}={\displaystyle \frac{Kf\left(S^{*}\right)}{1-\theta e^{-\mu \tau }}}.$

We remark that, as $\theta \in (0,1)$, we always have $1-\theta e^{-\mu \tau }>0$.

Assume that $I^{*}=0$. Then, the first equation of (7) implies that $S^{*}$ satisfies the equation:

(8)$\mathsf{\Phi }\left(S^{*}\right):=S^{*}+{\displaystyle \frac{K(1-e^{-\mu \tau })}{\mu (1-\theta e^{-\mu \tau })}}f\left(S^{*}\right)={\displaystyle \frac{\sigma }{\mu }}.$

Since $\mathsf{\Phi }$ is increasing, tends to infinity, and satisfies $\mathsf{\Phi }\left(0\right)=0$, we obtain the existence of a unique $S^{*}$ such that $\mathsf{\Phi }\left(S^{*}\right)=\sigma /\mu$. We can therefore write $S^{*}={\mathsf{\Phi }}^{-1}(\sigma /\mu )$. Thus, we obtain the disease-free steady-state:

(9)$(S^{*},I^{*},u^{*})=\left({\mathsf{\Phi }}^{-1}\left(\sigma /\mu \right),0,{\displaystyle \frac{Kf\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right)}{1-\theta e^{-\tau \mu }}}\right).$

Now, suppose that $S^{*}=\mu /\beta$. We immediately obtain that

$u^{*}={\displaystyle \frac{Kf(\mu /\beta )}{1-\theta e^{-\mu \tau }}}.$

We inject these two last expressions into the first equation of (7) to obtain

$I^{*}={\displaystyle \frac{\sigma }{\mu }}-\mathsf{\Phi }\left({\displaystyle \frac{\mu }{\beta }}\right),$

where the function $\mathsf{\Phi }$ is given by the expression (8). Then, $I^{*}>0$ exists if and only if

${\mathsf{\Phi }}^{-1}\left({\displaystyle \frac{\sigma }{\mu }}\right)>{\displaystyle \frac{\mu }{\beta }}.$

Let us define the basic reproduction number ${\mathcal{R}}_0$ as a threshold for the existence of the endemic steady-state:

${\mathcal{R}}_0:={\displaystyle \frac{\beta }{\mu }}{\mathsf{\Phi }}^{-1}\left({\displaystyle \frac{\sigma }{\mu }}\right).$

Hence, we obtain the existence of a unique endemic steady-state:

(10)$(S^{*},I^{*},u^{*})=\left({\displaystyle \frac{\mu }{\beta }},{\displaystyle \frac{\sigma }{\mu }}-\mathsf{\Phi }\left({\displaystyle \frac{\mu }{\beta }}\right),{\displaystyle \frac{Kf(\mu /\beta )}{1-\theta e^{-\mu \tau }}}\right)$

if and only if

${\mathcal{R}}_0>1.$

We will need to know the variation of ${\mathcal{R}}_0$ as a function of $\tau \in [0,+\infty )$. We have the following lemma.

Lemma 1.

The basic reproduction number $\tau \in [0,+\infty )\mapsto {\mathcal{R}}_0\left(\tau \right)$ is a decreasing function such that

${\displaystyle {\mathcal{R}}_0\left(0\right)=\frac{\beta \sigma }{{\mu }^2}and{\mathcal{R}}_0(\infty )=\frac{\beta }{\mu }{\mathsf{\Phi }}_{\infty }^{-1}\left(\frac{\sigma }{\mu }\right),}$

with

${\mathsf{\Phi }}_{\infty }\left(S\right):=S+\frac{K}{\mu }f\left(S\right).$

Proof. 

Consider the dependence of the function $\mathsf{\Phi }$ on $\tau$ by putting

$\mathsf{\Phi }(S,\tau ):=S+{\displaystyle \frac{K(1-e^{-\mu \tau })}{\mu (1-\theta e^{-\mu \tau })}}f\left(S\right).$

We want to study the variation of ${\displaystyle \tau \mapsto {\mathcal{R}}_0\left(\tau \right)=\frac{\beta }{\mu }{\mathsf{\Phi }}^{-1}\left(\frac{\sigma }{\mu },\tau \right)}$. Then, we set $S\left(\tau \right):={\mathsf{\Phi }}^{-1}\left(\sigma /\mu ,\tau \right)$. In fact, we have $\mathsf{\Phi }\left(S\right(\tau ),\tau )=\sigma /\mu$. By differentiating, we deduce that

$\frac{\partial \mathsf{\Phi }}{\partial S}(S\left(\tau \right),\tau )\frac{d}{d\tau }S\left(\tau \right)+\frac{\partial \mathsf{\Phi }}{\partial \tau }(S\left(\tau \right),\tau )=0.$

Thus, we obtain

$\frac{d}{d\tau }{\mathcal{R}}_0\left(\tau \right)=\frac{\beta }{\mu }\frac{d}{d\tau }S\left(\tau \right)=-\frac{\beta }{\mu }{\displaystyle \frac{\frac{\partial \mathsf{\Phi }}{\partial \tau }(S\left(\tau \right),\tau )}{\frac{\partial \mathsf{\Phi }}{\partial S}(S\left(\tau \right),\tau )}}<0.$

Then, the function $\tau \mapsto {\mathcal{R}}_0\left(\tau \right)$ is decreasing on $[0,+\infty )$. Furthermore, by taking $\tau =0$ and $\tau \to +\infty$ in the equality:

$\mathsf{\Phi }(S,\tau ):=S\left(\tau \right)+{\displaystyle \frac{K(1-e^{-\mu \tau })}{\mu (1-\theta e^{-\mu \tau })}}f\left(S\left(\tau \right)\right)={\displaystyle \frac{\sigma }{\mu }},$

we obtain that

${\displaystyle {\mathcal{R}}_0\left(0\right)=\frac{\beta \sigma }{{\mu }^2}\mathrm{and}{\mathcal{R}}_0(\infty )=\frac{\beta }{\mu }{\mathsf{\Phi }}_{\infty }^{-1}\left(\frac{\sigma }{\mu }\right),}$

with

${\mathsf{\Phi }}_{\infty }\left(S\right):=S+\frac{K}{\mu }f\left(S\right).$

□

We conclude from Lemma 1 that ${\mathcal{R}}_0(\infty )<{\mathcal{R}}_0\left(\tau \right)<{\mathcal{R}}_0\left(0\right)$, for all $\tau \in (0,+\infty )$, and we immediately obtain the following result.

Proposition 2.

Suppose that $\tau \ge 0$. Then, we have three cases:

• If ${\mathcal{R}}_0\left(0\right)<1$, then we have ${\mathcal{R}}_0\left(\tau \right)<1$, for all $\tau \ge 0$.

• If ${\mathcal{R}}_0(\infty )>1$, then we have ${\mathcal{R}}_0\left(\tau \right)>1$, for all $\tau \ge 0$.

• If ${\mathcal{R}}_0\left(0\right)>1$ and ${\mathcal{R}}_0(\infty )<1$, then there exists a unique $\overline{\tau }>0$ such that ${\mathcal{R}}_0\left(\tau \right)>1$ for $0\le \tau <\overline{\tau }$ and ${\mathcal{R}}_0\left(\tau \right)<1$ for $\tau >\overline{\tau }$, with ${\mathcal{R}}_0\left(\overline{\tau }\right)=1$.

We linearize System (4) around any equilibrium $(S^{*},I^{*},u^{*})$, and we obtain

$\begin{array}{c}\hfill \left\{\begin{array}{ccc}{S}^{\prime }\left(t\right)\hfill & =& -\beta I^{*}S\left(t\right)-\beta S^{*}I\left(t\right)-\mu S\left(t\right)-K{f}^{\prime }\left(S^{*}\right)S\left(t\right)+(1-\theta )e^{-\mu \tau }u(t-\tau ),\hfill \\ I^{\prime }\left(t\right)\hfill & =& \beta I^{*}S\left(t\right)+\beta S^{*}I\left(t\right)-\mu I\left(t\right),\hfill \\ u\left(t\right)\hfill & =& K{f}^{\prime }\left(S^{*}\right)S\left(t\right)+\theta e^{-\mu \tau }u(t-\tau ).\hfill \end{array}\right.\end{array}$

We look for solutions of the form $S\left(t\right)=e^{-\lambda t}{S}_0$, $I\left(t\right)=e^{-\lambda t}{I}_0$ and $u\left(t\right)=e^{-\lambda t}{u}_0$, $t>0$, and $\lambda \in \mathbb{C}$. Then, we obtain the following characteristic equation

$\begin{array}{c}\hfill \mathsf{\Delta }(\lambda ,\tau ):=\left|\begin{array}{ccc}\lambda +\beta I^{*}+\mu +K{f}^{\prime }\left(S^{*}\right)& \beta S^{*}& -(1-\theta )e^{-\mu \tau }{e}^{-\lambda \tau }\\ -\beta I^{*}& \lambda -\beta S^{*}+\mu & 0\\ -K{f}^{\prime }\left(S^{*}\right)& 0& 1-\theta e^{-\mu \tau }{e}^{-\lambda \tau }\end{array}\right|=0.\end{array}$

We develop the determinant and obtain

$\begin{array}{ccc}{\displaystyle \mathsf{\Delta }(\tau ,\lambda )}\hfill & =\hfill & (\lambda -\beta S^{*}+\mu )\left[(\lambda +\mu )(1-\theta e^{-\mu \tau }{e}^{-\lambda \tau })+K{f}^{\prime }\left(S^{*}\right)(1-e^{-\mu \tau }{e}^{-\lambda \tau })\right]\hfill \\ \hfill & \hfill & +(\lambda +\mu )\beta I^{*}(1-\theta e^{-\mu \tau }{e}^{-\lambda \tau }),\hfill \\ \hfill & =\hfill & 0.\hfill \end{array}$

First, consider the disease-free equilibrium (9). Then, the characteristic equation becomes

$\begin{array}{c}[\lambda -\mu ({\mathcal{R}}_0-1)]\left[(\lambda +\mu )(1-\theta e^{-\mu \tau }{e}^{-\lambda \tau })+K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}\left(\sigma /\mu \right)\right)(1-e^{-\mu \tau }{e}^{-\lambda \tau })\right]=0.\hfill \end{array}$

We have an immediate real eigenvalue

${\lambda }_0=\mu ({\mathcal{R}}_0-1).$

The other eigenvalues are solutions of the transcendental equation:

(11)$\overline{\mathsf{\Delta }}\left(\lambda \right):=P\left(\lambda \right)-Q\left(\lambda \right)e^{-\lambda \tau }=0,$

with

$P\left(\lambda \right):=\lambda +\mu +K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right)\mathrm{and}Q\left(\lambda \right):=e^{-\mu \tau }[\theta \lambda +\theta \mu +K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right)].$

We have the following properties:

Property 1.

• (i) 

  The equation (11) has no purely imaginary root.

• (ii) 

  All roots ${\lambda }_n$, $n\in \mathbb{N}$, of (11) such that $|{\lambda }_n|\to +\infty$, have a negative real part for n large enough:

Proof. 

(i) First, we seek for pure imaginary roots $\lambda =\pm i\omega$, $\omega >0$, of (11). Splitting in (11) the real and the imaginary part, we obtain

$\left\{\begin{array}{ccc}{\displaystyle \gamma cos\left(\omega \tau \right)e^{-\mu \tau }+\theta \omega sin\left(\omega \tau \right)e^{-\mu \tau }}\hfill & =& \mu +K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right),\hfill \\ -\theta \omega cos\left(\omega \tau \right)e^{-\mu \tau }+\gamma sin\left(\omega \tau \right)e^{-\mu \tau }\hfill & =& -\omega ,\hfill \end{array}\right.$

where $\gamma :=\theta \mu +K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right)$. Solving this system, we find that

$\left\{\begin{array}{ccc}{\displaystyle cos\left(\omega \tau \right)}\hfill & =& {\displaystyle \frac{1}{({\theta }^2{\omega }^2+{\gamma }^2)e^{-\mu \tau }}}\left[\theta {\omega }^2+\gamma (\mu +K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right)\right],\hfill \\ {\displaystyle sin\left(\omega \tau \right)}\hfill & =& {\displaystyle \frac{1}{({\theta }^2{\omega }^2+{\gamma }^2)e^{-\mu \tau }}}\left[-\omega \gamma +\theta \omega (\mu +K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right)\right].\hfill \end{array}\right.$

Using the trigonometric identity ${cos}^2\left(\omega \tau \right)+{sin}^2\left(\omega \tau \right)=1$, we obtain that $\omega$ satisfies

$({\theta }^2{\omega }^2+{\gamma }^2){\left(e^{-\mu \tau }\right)}^2={\omega }^2+{(\mu +K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right))}^2.$

Then, we obtain

(12)${\omega }^2={\displaystyle \frac{(\gamma e^{-\mu \tau }-(\mu +K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right)))((\mu +K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right))+\gamma e^{-\mu \tau })}{1-{\theta }^2{e}^{-2\mu \tau }}}.$

Notice that $1-{\theta }^2{e}^{-2\mu \tau }>0$ and

$\gamma e^{-\mu \tau }-(\mu +K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right))=[\theta \mu +K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right)]e^{-\mu \tau }-[\mu +K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right)]<0.$

Thus, the right-hand side of (12) is negative, which is absurd. We conclude that the Equation (11) has no purely imaginary root.

(ii) As the function $\overline{\mathsf{\Delta }}$ is entire, the set of its roots cannot have accumulation points. Furthermore, $\overline{\mathsf{\Delta }}$ has at most countably many zeros. If $\overline{\mathsf{\Delta }}$ has infinitely many different zeros, ${\lambda }_n\in \mathbb{C}$, $n\in \mathbb{N}$, then

$\underset{n\to +\infty }{lim}\left|{\lambda }_n\right|=+\infty .$

The sequence ${\lambda }_n$, $n\in \mathbb{N}$, satisfies

$e^{{\lambda }_n\tau }\left({\displaystyle \frac{P\left({\lambda }_n\right)}{{\lambda }_n}}\right)-{\displaystyle \frac{Q\left({\lambda }_n\right)}{{\lambda }_n}}=0.$

It is clear that

$\underset{n\to \infty }{lim}\left|{\displaystyle \frac{P\left({\lambda }_n\right)}{{\lambda }_n}}-1\right|=0\mathrm{and}\underset{n\to \infty }{lim}\left|{\displaystyle \frac{Q\left({\lambda }_n\right)}{{\lambda }_n}}-\theta e^{-\mu \tau }\right|=0.$

Then, there is a sequence ${\left({\lambda }_n^{\prime }\right)}_{n\in \mathbb{N}}$ of roots of the equation

$e^{\lambda \tau }-\theta e^{-\mu \tau }=0$

such that

$\underset{n\to \infty }{lim}({\lambda }_n-{\lambda }_n^{\prime })=0.$

This means that the closed-forms of the roots ${\lambda }_n$ are the complex numbers

${\overline{\lambda }}_p={\displaystyle \frac{1}{\tau }}ln\left(\theta \right)-\mu +{\displaystyle \frac{2p\pi }{\tau }}i,p\in \mathbb{Z}.$

Since $\theta \in (0,1)$, we have that $\mathcal{R}e\left({\overline{\lambda }}_p\right)<0$. Thus, all roots such that $|{\lambda }_n|\to +\infty$ have a negative real part for n large enough. □

We state the stability of the disease-free equilibrium (9):

Corollary 1.

• (i) 

  Suppose that there exists ${\tau }_0\ge 0$ such that ${\mathcal{R}}_0\left({\tau }_0\right)>1$. Then, the disease-free equilibrium (9) is unstable for $\tau ={\tau }_0$.

• (ii) 

  Suppose that there exists ${\tau }_1\ge 0$ such that ${\mathcal{R}}_0\left({\tau }_1\right)<1$ and the disease-free equilibrium (9) is locally asymptotically stable for $\tau ={\tau }_1$. Then, it is locally asymptotically stable for all $\tau \ge {\tau }_1$.

• (iii) 

  Suppose that ${\mathcal{R}}_0\left(0\right)<1$. Then, the disease-free equilibrium (9) is locally asymptotically stable for all $\tau \ge 0$.

Proof. 

(i) The real eigenvalue

${\lambda }_0=\mu ({\mathcal{R}}_0\left({\tau }_0\right)-1)$

is positive because ${\mathcal{R}}_0\left({\tau }_0\right)>1$. Then, the disease-free equilibrium (9) is unstable for $\tau ={\tau }_0$.

(ii) Suppose that there exits ${\tau }_1\ge 0$ such that ${\mathcal{R}}_0\left({\tau }_1\right)<1$ and the disease-free equilibrium (9) is locally asymptotically stable for $\tau ={\tau }_1$. As $\tau \mapsto {\mathcal{R}}_0\left(\tau \right)$ is decreasing, we have ${\mathcal{R}}_0\left(\tau \right)<1$ for all $\tau \ge {\tau }_1$. Then, the real eigenvalue ${\lambda }_0$ is negative for all $\tau \ge {\tau }_1$. The other eigenvalues given by

$\overline{\mathsf{\Delta }}(\tau ,\lambda ):=P(\tau ,\lambda )-Q(\tau ,\lambda )e^{-\lambda \tau }=0,$

with

$P(\tau ,\lambda ):=\lambda +\mu +K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right)\mathrm{and}Q(\tau ,\lambda ):=e^{-\mu \tau }[\theta \lambda +\theta \mu +K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right)],$

which cannot cross the imaginary axis when we increase $\tau$ starting from ${\tau }_1$. The disease-free equilibrium (9), thus, keeps the same stability as for $\tau ={\tau }_1$, i.e., it remains locally asymptotically stable for all $\tau \le {\tau }_1$.

(iii) Let us consider the case $\tau =0$. Then, Equation (11) becomes

$(\lambda +\mu )(1-\theta )=0.$

This means that $\lambda =-\mu <0$. Then, the disease-free equilibrium (9) is locally asymptotically stable for $\tau =0$. We use the previous property, (ii), to conclude that the disease-free equilibrium (9) is locally asymptotically stable for all $\tau \ge 0$.

One can prove this corollary directly by proving that Equation (11) has no root with a nonnegative real part, and therefore, the stability of the endemic equilibrium is determined by the sign of ${\mathcal{R}}_0-1$. Indeed, let us suppose by contradiction the existence of such a root $\lambda =\alpha +i\omega$, with $\alpha \ge 0$. Then, we obtain

${\displaystyle \frac{\left|P(\alpha +i\omega )\right|}{\left|Q(\alpha +i\omega )\right|}}=e^{-\alpha \tau }\le 1.$

This gives the following inequality:

$\left[P\left(\alpha \right)-Q\left(\alpha \right)\right]\left[P\left(\alpha \right)+Q\left(\alpha \right)\right]+\left(1-e^{-2\mu \tau }{\theta }^2\right){\omega }^2\le 0,$

with

$P\left(\alpha \right)-Q\left(\alpha \right)=(1-e^{-\mu \tau }\theta )(\alpha +\mu )+(1-e^{-\mu \tau })K{f}^{\prime }\left({\mathsf{\Phi }}^{-1}(\sigma /\mu )\right)>0.$

This leads to a contradiction. □

Now, we focus on the endemic equilibrium (10):

$(S^{*},I^{*},u^{*})=\left({\displaystyle \frac{\mu }{\beta }},{\displaystyle \frac{\sigma }{\mu }}-\mathsf{\Phi }\left({\displaystyle \frac{\mu }{\beta }}\right),{\displaystyle \frac{Kf(\mu /\beta )}{1-\theta e^{-\mu \tau }}}\right).$

Proposition 3.

Suppose that ${\mathcal{R}}_0\left(0\right)>1$. Then, for all $\tau \ge 0$ such that ${\mathcal{R}}_0\left(\tau \right)>1$, the endemic steady-state is locally asymptotically stable.

Proof. 

We suppose that ${\mathcal{R}}_0\left(\tau \right)>1$, and we consider the endemic equilibrium. The characteristic equation becomes

(13)$\begin{array}{c}\mathsf{\Delta }(\tau ,\lambda ):={\lambda }^2+A\lambda +\mu \beta I^{*}-e^{-\mu \tau }{e}^{-\lambda \tau }\left[\theta {\lambda }^2+(\theta A+K{f}^{\prime }\left(S^{*}\right))\lambda +\theta \mu \beta I^{*}\right]=0,\hfill \end{array}$

with $A:=\mu +\beta I^{*}>0$. We follow the same steps as in the proof of the stability of the disease-free equilibrium. First, we consider the case $\tau =0$. Then, the endemic steady-state becomes

$(S^{*},I^{*},u^{*})=\left({\displaystyle \frac{\mu }{\beta }},{\displaystyle \frac{\sigma }{\mu }}-{\displaystyle \frac{\mu }{\beta }},{\displaystyle \frac{Kf(\mu /\beta )}{1-\theta }}\right),$

with ${\mathcal{R}}_0\left(0\right)=\beta \sigma /{\mu }^2$, and the characteristic equation is reduced to

$\begin{array}{c}\hfill \mathsf{\Delta }(\lambda ,\tau ):=(1-\theta )\left|\begin{array}{cc}\lambda +\beta I^{*}+\mu & \beta S^{*}\\ -\beta I^{*}& \lambda -\beta S^{*}+\mu \end{array}\right|=0.\end{array}$

We obtain two eigenvalues ${\lambda }_1=-\mu <0$ and ${\lambda }_2=-\mu -\beta (I^{*}-S^{*})=-\mu ({\mathcal{R}}_0\left(0\right)-1)<0$. Hence, the endemic steady-state is locally asymptotically stable for $\tau =0$. We then look for pure imaginary roots $\lambda =\pm i\omega$, $\omega >0$. Then, splitting the real and imaginary parts in (13), we have the following system:

$\left\{\begin{array}{ccc}\omega \mathsf{\Gamma }cos\left(\omega \tau \right)e^{-\mu \tau }+{\rho }_{\omega }sin\left(\omega \tau \right)e^{-\mu \tau }\hfill & =& \omega (A+\psi f^{\prime }(\mu /\beta )),\hfill \\ {\rho }_{\omega }cos\left(\omega \tau \right)e^{-\mu \tau }-\omega \mathsf{\Gamma }sin\left(\omega \tau \right)e^{-\mu \tau }\hfill & =& {\omega }^2-\mu \beta I^{*},\hfill \end{array}\right.$

with ${\rho }_{\omega }:=({\omega }^2-\mu \beta I^{*})\theta e^{-\mu \tau }$ and $\mathsf{\Gamma }=(\theta A+\psi f^{\prime }(\mu /\beta ))e^{-\mu \tau }$. Using the identity ${cos}^2\left(\omega \tau \right)+{sin}^2\left(\omega \tau \right)=1$, we obtain the polynomial equation in $\omega$:

$\begin{array}{ccc}\hfill [1-{\left(\theta e^{-\mu \tau }\right)}^2]{\omega }^4+[2\mu \beta I^{*}{\left(\theta e^{-\mu \tau }\right)}^2-2\mu \beta I^{*}+{(A+\psi f^{\prime }(\mu /\beta ))}^2-{\mathsf{\Gamma }}^2]{\omega }^2& & \\ \hfill +{\left(\mu \beta I^{*}\right)}^2(1-{\left(\theta e^{-\mu \tau }\right)}^2)& =& 0.\hfill \end{array}$

Let us set

$N:=\frac{1}{1-{\left(\theta e^{-\mu \tau }\right)}^2}[2\mu \beta I^{*}{\left(\theta e^{-\mu \tau }\right)}^2-2\mu \beta I^{*}+{(A+\psi f^{\prime }(\mu /\beta ))}^2-{\mathsf{\Gamma }}^2]\mathrm{and}X:={\omega }^2.$

Then, we obtain the following polynomial:

(14)$X^2+NX+{\left(\mu \beta I^{*}\right)}^2=0.$

First, if $N>0$, then the Routh–Hurwitz criterion implies that all roots of (14) have a negative real part. This leads to a contradiction.

Suppose now that $N\le 0$. We compute the discriminant of Equation (14) to obtain

$\mathsf{\Delta }=N^2-4{\left(\mu \beta I^{*}\right)}^2=(N+2\mu \beta I^{*})(N-2\mu \beta I^{*}).$

As we supposed that $N\le 0$, we have that $N-\mu \beta I^{*}<0$. On the other hand, we have that

$(1-{\left(\theta e^{-\mu \tau }\right)}^2)(N+2\mu \beta I^{*})={(A+K{f}^{\prime }(\mu /\beta ))}^2-{\mathsf{\Gamma }}^2.$

But, as $\theta <1$ and $e^{-\mu \tau }<1$, we obtain

${\mathsf{\Gamma }}^2={\left[(\theta A+K{f}^{\prime }(\mu /\beta ))e^{-\mu \tau }\right]}^2<{(A+K{f}^{\prime }(\mu /\beta ))}^2.$

We conclude that $\mathsf{\Delta }<0$. Therefore, there is no positive real root of (14) and $\omega >0$ cannot exist. Furthermore, as for the disease-free equilibrium, if $\mathsf{\Delta }$ has infinitely many different zeros, ${\lambda }_n\in \mathbb{C}$, $n\in \mathbb{N}$, then

$\underset{n\to +\infty }{lim}\left|{\lambda }_n\right|=+\infty .$

The sequence ${\lambda }_n$, $n\in \mathbb{N}$, satisfies

$e^{{\lambda }_n\tau }\left({\displaystyle \frac{P\left({\lambda }_n\right)}{{\lambda }_n^2}}\right)-{\displaystyle \frac{Q\left({\lambda }_n\right)}{{\lambda }_n^2}}=0,$

with

$P\left(\lambda \right):={\lambda }^2+A\lambda +\mu \beta I^{*}\mathrm{and}Q\left(\lambda \right):=e^{-\mu \tau }\left[\theta {\lambda }^2+(\theta A+K{f}^{\prime }\left(S^{*}\right))\lambda +\theta \mu \beta I^{*}\right].$

We have

$\underset{n\to \infty }{lim}\left|{\displaystyle \frac{P\left({\lambda }_n\right)}{{\lambda }_n^2}}-1\right|=0\mathrm{and}\underset{n\to \infty }{lim}\left|{\displaystyle \frac{Q\left({\lambda }_n\right)}{{\lambda }_n^2}}-\theta e^{-\mu \tau }\right|=0.$

Then, as for the case of the disease-free steady-state, the closed forms of the roots ${\lambda }_n$ are

${\overline{\lambda }}_p={\displaystyle \frac{1}{\tau }}ln\left(\theta \right)-\mu +{\displaystyle \frac{2p\pi }{\tau }}i,p\in \mathbb{Z}.$

Thus, all roots such that $|{\lambda }_n|\to +\infty$ have a negative real part for n large enough. We conclude that all roots of (13) have a negative real part, and then, the endemic equilibrium is locally asymptotically stable. □

Remark 1.

As for the proof of Corollary 1, we can also prove directly by contradiction that there is no root $\lambda =\alpha +i\omega$ of (13) with $\alpha \ge 0$.

4. Global Stability Analysis

We consider in this section the case:

$f\left(S\right)=S.$

The model becomes

(15)$\left\{\begin{array}{ccc}{S}^{\prime }\left(t\right)\hfill & =& {\displaystyle \sigma -\beta I\left(t\right)S\left(t\right)-\mu S\left(t\right)-KS\left(t\right)+(1-\theta )e^{-\mu \tau }u(t-\tau ),}\hfill \\ I^{\prime }\left(t\right)\hfill & =& {\displaystyle \beta I\left(t\right)S\left(t\right)-\mu I\left(t\right),}\hfill \\ u\left(t\right)\hfill & =& KS\left(t\right)+\theta e^{-\mu \tau }u(t-\tau ),\hfill \end{array}\right.$

with the initial condition:

$S\left(0\right)=S_0,I\left(0\right)=I_0\mathrm{and}u\left(t\right)=\phi \left(t\right),t\in [-\tau ,0].$

The function $\mathsf{\Phi }$ given by (8) becomes

$\mathsf{\Phi }\left(S^{*}\right)=\left[1+{\displaystyle \frac{K(1-e^{-\mu \tau })}{\mu (1-\theta e^{-\mu \tau })}}\right]S^{*}.$

The corresponding disease-free equilibrium is

(16)$(S^0,I^0,u^0):=\left(\frac{\sigma (1-\theta e^{-\mu \tau })}{\mu +K-(\mu \theta +K)e^{-\mu \tau }},0,\frac{K\sigma }{\mu +K-(\mu \theta +K)e^{-\mu \tau }}\right),$

and the endemic steady-state is

(17)$(S^{*},I^{*},u^{*})=\left({\displaystyle \frac{\mu }{\beta }},{\displaystyle \frac{\sigma }{\mu }}-\mathsf{\Phi }\left({\displaystyle \frac{\mu }{\beta }}\right),{\displaystyle \frac{K\mu }{\beta (1-\theta e^{-\mu \tau })}}\right).$

The basic reproduction number becomes

(18)${\mathcal{R}}_0:={\displaystyle \frac{\beta }{\mu }}{\mathsf{\Phi }}^{-1}\left({\displaystyle \frac{\sigma }{\mu }}\right)={\displaystyle \frac{\beta }{\mu }}{S}^0={\displaystyle \frac{\beta \sigma (1-\theta e^{-\mu \tau })}{{\mu }^2(1-\theta e^{-\mu \tau })+\mu K(1-e^{-\mu \tau })}}.$

A model similar to (15) was studied in [11] as a generalization of a Kermack–McKendrick SIR model with an age-structured protection phase. We established the global asymptotic stability of the two steady-states. We showed that if ${\mathcal{R}}_0<1$, then the disease-free steady-state is globally asymptotically stable, and if ${\mathcal{R}}_0>1$, then the endemic steady-state is globally asymptotically stable. We constructed quadratic and logarithmic Lyapunov functions to establish this global stability. In the following subsections, we also present results on the global asymptotic stability of the two steady-states. We used similar techniques as those used for the proofs of the global stability results in [11].

4.1. Global Asymptotic Stability of the Disease-Free Steady-State

In this subsection, we assumed that ${\mathcal{R}}_0<1$, and we show that the disease-free equilibrium $(S^0,I^0,u^0)$ is globally asymptotically stable. The proof is based on the use of the following auxiliary differential–difference system, for $t>0$:

(19)$\left\{\begin{array}{ccc}{\displaystyle \frac{d{S}^{+}}{dt}}\left(t\right)\hfill & =& \sigma -(\mu +K)S^{+}\left(t\right)+(1-\theta )e^{-\mu \tau }{u}^{+}(t-\tau ),\hfill \\ u^{+}\left(t\right)\hfill & =& K{S}^{+}\left(t\right)+\theta e^{-\mu \tau }{u}^{+}(t-\tau ),\hfill \end{array}\right.$

with the same initial condition as the main system:

$S^{+}\left(0\right)=S_0\mathrm{and}{u}^{+}\left(t\right)=\phi \left(t\right),t\in [-\tau ,0].$

By using a comparison principle, we obtain that $S\left(t\right)\le S^{+}\left(t\right)$ and $u\left(t\right)\le u^{+}\left(t\right)$ for all $t>0$. Furthermore, the system (19) has a unique equilibrium that corresponds to the disease-free equilibrium of the system (15). The basic reproduction number associated with (19) is also given by (18). The global asymptotic stability of the auxiliary system (19) was established in [11]. Indeed, by setting

$\widehat{S}\left(t\right):=S^{+}\left(t\right)-S^0,\widehat{u}\left(t\right):=u^{+}\left(t\right)-u^0$

and considering the Lyapunov function $V:\mathbb{R}\times \mathcal{C}([-\tau ,0],\mathbb{R})⟶{\mathbb{R}}^{+}$,

$V(S_0,\phi )=\frac{S_0^2}{2}+\xi {\int }_{-\tau }^0{\phi }^2\left(s\right)ds,\mathrm{with}\xi ={\displaystyle \frac{1}{2K^2}}\left[\mu (1-{\theta }^2{e}^{-2\mu \tau })+K\right]>0,$

we proved in [11] the following lemma.

Lemma 2.

Suppose that ${\mathcal{R}}_0<1$. Then, the unique steady-state $(S^0,u^0)$ of the system (19) is globally asymptotically stable.

Thanks to a comparison principle and the global attractivity of the set:

$\begin{array}{ccc}{\mathsf{\Omega }}_{\epsilon }\hfill & :=& \left\{(S,I,u)\in {\mathbb{R}}^{+}\times {\mathbb{R}}^{+}\times \mathcal{C}([-\tau ,0],{\mathbb{R}}^{+}):\right.\hfill \\ \hfill & & \left.0\le S\le S^0+\epsilon ,0\le u\left(s\right)\le u^0+\epsilon ,s\in [-\tau ,0]\right\},\hfill \end{array}$

for $\epsilon >0$ small enough, we are able to obtain the global asymptotic stability of the disease-free equilibrium of (15). More precisely, the proof of the global asymptotic stability of the disease-free steady-state (16) is based on the following lemma.

Lemma 3.

Suppose that ${\mathcal{R}}_0<1$. Then, for any sufficiently small $\epsilon >0$, the subset ${\mathsf{\Omega }}_{\epsilon }$ of ${\mathbb{R}}^{+}\times {\mathbb{R}}^{+}\times \mathcal{C}([-\tau ,0],{\mathbb{R}}^{+})$ is a global attractor for the system (15).

Proof. 

The solutions of (15) satisfy, for all $t>0$, the system:

$\left\{\begin{array}{ccc}{S}^{\prime }\left(t\right)\hfill & \le & \sigma -(\mu +K)S\left(t\right)+(1-\theta )e^{-\mu \tau }u(t-\tau ),\hfill \\ u\left(t\right)\hfill & =& KS\left(t\right)+\theta e^{-\mu \tau }u(t-\tau ).\hfill \end{array}\right.$

By the comparison principle and the positivity of the solutions, we have $0\le S\left(t\right)\le S^{+}\left(t\right)$ and $0\le u\left(t\right)\le u^{+}\left(t\right)$, for all $t>0$, where $(S^{+},u^{+})$ is the solution of the system (19). Lemma 2 implies that $(S^{+}\left(t\right),u^{+}\left(t\right))\to (S^0,u^0)$ as $t\to +\infty$. This convergence means that the subset ${\mathsf{\Omega }}_{\epsilon }$, with $\epsilon >0$ small enough, is a global attractor for the system (15). □

Lemma 3 allows us to restrict the analysis of the global stability of the disease-free steady-state (16) of (15) to the subset ${\mathsf{\Omega }}_{\epsilon }$ with $\epsilon >0$ small enough.

Theorem 1.

Suppose that ${\mathcal{R}}_0<1$. Then, the disease-free steady-state $(S^0,0,u^0)$ given by (16) of the model (15) is globally asymptotically stable.

Proof. 

We consider the solutions of System (15) in the subset ${\mathsf{\Omega }}_{\epsilon }$, with $\epsilon >0$ small enough. It is clear that the second equation of (15) implies that

${\displaystyle I^{\prime }\left(t\right)\le \beta (S^0+\epsilon )I\left(t\right)-\mu I\left(t\right)=-\mu \left(1-\frac{\beta (S^0+\epsilon )}{\mu }\right)I\left(t\right).}$

Since ${\mathcal{R}}_0={\displaystyle \frac{\beta }{\mu }}{S}^0<1$, we can choose $\epsilon >0$ such that ${\displaystyle \frac{\beta (S^0+\epsilon )}{\mu }}<1$. This implies that ${lim}_{t\to +\infty }I\left(t\right)=0$. Then, there exists $T_{\epsilon }>0$ such that $0\le I\left(t\right)\le \epsilon$, for all $t>T_{\epsilon }$. We then have, for $t>T_{\epsilon }$,

$\left\{\begin{array}{ccc}{S}^{\prime }\left(t\right)\hfill & \ge & \sigma -(\mu +K)S\left(t\right)-\epsilon \beta S\left(t\right)+(1-\theta )e^{-\mu \tau }u(t-\tau ),\hfill \\ u\left(t\right)\hfill & =& KS\left(t\right)+\theta e^{-\mu \tau }u(t-\tau ).\hfill \end{array}\right.$

We again use the comparison principle to obtain $S\left(t\right)\ge S_{\epsilon }\left(t\right)$ and $u\left(t\right)\ge u_{\epsilon }\left(t\right)$, for all $t>T_{\epsilon }$, where $(S_{\epsilon },u_{\epsilon })$ is the solution of the problem, for $t>0$,

(20)$\left\{\begin{array}{ccc}\hfill S_{\epsilon }^{\prime }\left(t\right)& =& \sigma -(\mu +K)S_{\epsilon }\left(t\right)-\epsilon \beta S_{\epsilon }\left(t\right)+(1-\theta )e^{-\mu \tau }{u}_{\epsilon }(t-\tau ),\hfill \\ \hfill u_{\epsilon }\left(t\right)& =& K{S}_{\epsilon }\left(t\right)+\theta e^{-\mu \tau }{u}_{\epsilon }(t-\tau ),\hfill \\ \hfill S_{\epsilon }\left(0\right)& =& S_0,u_{\epsilon }\left(s\right)=\phi \left(s\right),\mathrm{for}-\tau \le s\le 0.\hfill \end{array}\right.$

We use the same techniques as for the proof of the lemmas 2 and 3 to show that $(S_{\epsilon }\left(t\right),u_{\epsilon }\left(t\right))\to (S_{\epsilon }^0,u_{\epsilon }^0)$ as $t\to +\infty$, with $(S_{\epsilon }^0,u_{\epsilon }^0)$ the steady-state of System (20). In fact, $(S_{\epsilon }^0,u_{\epsilon }^0)$ is given by the expression:

$(S_{\epsilon }^0,u_{\epsilon }^0):=\left(\frac{\sigma (1-\theta e^{-\mu \tau })}{{\mu }_{\epsilon }+K-({\mu }_{\epsilon }\theta +K)e^{-\mu \tau }},\frac{K\sigma }{{\mu }_{\epsilon }+K-({\mu }_{\epsilon }\theta +K)e^{-\mu \tau }}\right),$

with ${\mu }_{\epsilon }:=\mu +\epsilon \beta$. Then, there exists ${\tilde{T}}_{\epsilon }>T_{\epsilon }>0$ such that, for $t>{\tilde{T}}_{\epsilon }$,

$S_{\epsilon }^0-\epsilon \le S\left(t\right)\le S^0+\epsilon \mathrm{and}{u}_{\epsilon }^0-\epsilon \le u\left(t\right)\le u^0+\epsilon .$

As $\epsilon \to 0$, we have $S_{\epsilon }^0\to S^0$ and $u_{\epsilon }^0\to u^0$. Then, we obtain

$\underset{t\to +\infty }{lim}S\left(t\right)=S^0\mathrm{and}\underset{t\to +\infty }{lim}u\left(t\right)=u^0.$

Recall that $(S^0,0,u^0)$ is locally asymptotically stable. Then, it is globally asymptotically stable. □

4.2. Global Asymptotic Stability of the Endemic Equilibrium

In this subsection, we assumed that

${\mathcal{R}}_0>1,$

where ${\mathcal{R}}_0$ is given by (18), and we studied the global asymptotic stability of the endemic steady-state (17). Let us define

$\tilde{S}\left(t\right):=S\left(t\right)-S^{*}\mathrm{and}\tilde{u}\left(t\right):=u\left(t\right)-u^{*}.$

Then, the system (15) becomes, with $\beta S^{*}=\mu$,

(21)$\left\{\begin{array}{ccc}{\tilde{S}}^{\prime }\left(t\right)\hfill & =& -(\mu +K)\tilde{S}\left(t\right)-\beta \tilde{S}\left(t\right)I\left(t\right)-\beta S^{*}I\left(t\right)+\beta S^{*}{I}^{*}+(1-\theta )e^{-\mu \tau }\tilde{u}(t-\tau ),\hfill \\ I^{\prime }\left(t\right)\hfill & =& \beta \tilde{S}\left(t\right)I\left(t\right),\hfill \\ \tilde{u}\left(t\right)\hfill & =& K\tilde{S}\left(t\right)+\theta e^{-\mu \tau }\tilde{u}(t-\tau ).\hfill \end{array}\right.$

The endemic steady-state (17) of the system (15) becomes $(0,I^{*},0)$ (as a steady-state of the system (21)). We then obtain the following result.

Theorem 2.

Let us suppose that ${\mathcal{R}}_0>1$. Then, the endemic steady-state $(S^{*},I^{*},u^{*})$ given by (17) of the model (15) is globally asymptotically stable.

Proof. 

The proof of this theorem is based on the use of the following Lyapunov function $V:\mathbb{R}\times {\mathbb{R}}^{+}\times \mathcal{C}([-\tau ,0],\mathbb{R})\to {\mathbb{R}}^{+}$ defined by

$V(S_0,I_0,\phi )=\frac{S_0^2}{2}+\xi {\int }_{-\tau }^0{\phi }^2\left(s\right)ds+S^{*}\left[I_0-I^{*}-I^{*}ln\left(\frac{I_0}{I^{*}}\right)\right],$

with $\xi ={\displaystyle \frac{1}{2K^2}}\left[K+\mu (1-\left({\theta }^2{e}^{-2\mu \tau }\right)\right]$.

First, we note that the function $G:(0+\infty )\to [0+\infty )$ defined by

$G\left(I_0\right)=I_0-I^{*}-I^{*}ln\left({\displaystyle \frac{I_0}{I^{*}}}\right),I_0>0,$

satisfies $G\left(I_0\right)\ge 0$ for all $I_0>0$ and $G\left(I_0\right)=0$ if and only if $I_0=I^{*}$. Then, we have $V(S_0,I_0,u_0)=0$ if and only if $(S_0,I_0,u_0)=(0,I^{*},0)$. We set

$\tilde{V}\left(t\right):=V(\tilde{S}\left(t\right),I\left(t\right),{\tilde{u}}_t),t\ge 0,$

where $(\tilde{S}\left(t\right),I\left(t\right),{\tilde{u}}_t)$ is the solution of (21). Then, we obtain

(22)$\begin{array}{ccc}{\tilde{V}}^{\prime }\left(t\right)\hfill & =\hfill & -a{\tilde{S}}^2\left(t\right)+b\tilde{S}\left(t\right)\tilde{u}(t-\tau )-c{\tilde{u}}^2(t-\tau )-\beta I\left(t\right){\tilde{S}}^2\left(t\right),\hfill \\ \hfill & \le \hfill & {\displaystyle \frac{b^2-4ac}{4c}}{\tilde{S}}^2\left(t\right),\hfill \\ \hfill & =\hfill & -\kappa {\tilde{S}}^2\left(t\right),\hfill \end{array}$

with

$\left\{\begin{array}{ccc}a\hfill & =& \mu +K-\xi K^2,\hfill \\ b\hfill & =& (1-\theta )e^{-\mu \tau }+2\xi K\theta e^{-\mu \tau }>0,\hfill \\ c\hfill & =& \xi (1-\theta e^{-\mu \tau })(1+\theta e^{-\mu \tau })>0,\hfill \\ \kappa \hfill & =& (4ac-b^2)/4c>0.\hfill \end{array}\right.$

We conclude that the function $t\mapsto V(\tilde{S}\left(t\right),I\left(t\right),{\tilde{u}}_t)$ is nonincreasing. Then, we obtain

$\underset{t\to +\infty }{lim}V(\tilde{S}\left(t\right),I\left(t\right),{\tilde{u}}_t)=\underset{s\ge 0}{inf}V(\tilde{S}\left(s\right),I\left(s\right),{\tilde{u}}_s)=:V^{*}\ge 0.$

By integrating (22) from 0 to $t>0$, we obtain

$\kappa {\int }_0^t{\tilde{S}}^2\left(s\right)ds\le V(\tilde{S}\left(0\right),I\left(0\right),{\tilde{u}}_0)-V(\tilde{S}\left(t\right),I\left(t\right),{\tilde{u}}_t).$

It is clear that the limits on both sides exist when $t\to +\infty$, and we have

(23)$\underset{t\to +\infty }{lim}{\int }_0^t{\tilde{S}}^2\left(s\right)ds\le {\displaystyle \frac{1}{\kappa }}\left[V(\tilde{S}\left(0\right),I\left(0\right),{\tilde{u}}_0)-V^{*}\right].$

Furthermore, the first equation of System (21) implies that $t\mapsto {\tilde{S}}^{\prime }\left(t\right)$ is uniformly bounded. Then, the function $t\mapsto \tilde{S}\left(t\right)$ is uniformly continuous. We conclude from (23) that ${lim}_{t\to +\infty }\tilde{S}\left(t\right)$ = 0, and the fluctuations lemma implies that there exists a sequence ${\left(t_k\right)}_{k\in \mathbb{N}}$, $t_k\to +\infty$ such that ${lim}_{k\to +\infty }{\tilde{S}}^{\prime }\left(t_k\right)=0$. On the other hand, using the difference equation satisfied by $\tilde{u}$, we also have ${lim}_{t\to +\infty }\tilde{u}\left(t\right)=0$. Then, the first equation of System (21) implies that ${lim}_{k\to +\infty }I\left(t_k\right)=I^{*}$. We also have from the expression of V that

$\underset{t\to +\infty }{lim}G\left(I\left(t\right)\right)={\displaystyle \frac{V^{*}}{S^{*}}}.$

Furthermore, the continuity of the function G implies that ${lim}_{k\to +\infty }G\left(I\left(t_k\right)\right)=G\left(I^{*}\right)=0$. Then, $V^{*}=0$. This means that ${lim}_{t\to +\infty }G\left(I\left(t\right)\right)=0$. The properties of the function G imply that ${lim}_{t\to +\infty }I\left(t\right)=I^{*}$. We conclude that the endemic steady-state $(S^{*},I^{*},u^{*})$ is globally asymptotically stable. □

We assumed, for the first time, that the rate of the start of PrEP treatment was constant and that $f\left(S\right)=S$. Thus, the model admits two equilibrium points, the disease-free equilibrium and the endemic equilibrium. By considering ${\mathcal{R}}_0$ as a function of the delay, we showed that, if ${\mathcal{R}}_0<1$, the disease-free equilibrium is globally asymptotically stable and that, if ${\mathcal{R}}_0>1$, then it is the endemic equilibrium that becomes globally asymptotically stable.

5. Applications to French Clinical Data

In order to apply our model to real data, we chose the French MSM population, and later, we focused more specifically on the high-risk individuals within this group. We define here the high-risk MSM population as multi-partner gay men, as well as sexual intercourse while using drugs, usually called chemsex, or new generations not inclined to use any protection anymore (mainly condoms). Of course, this reduces our study to a small part of the population, but, because of their having the highest rate of infectivity, they belong to the target of PrEP treatments at the infectious disease center in France.

Furthermore, for a more realistic approach, $\psi$ was considered time-dependent, and $\psi$ is defined by the logistic equation for a more accurate data fit (details and parameters in the first subsection below). The PrEP users were collected between January 2016 (starting year of PrEP in France) and June 2019 (before the SARS-CoV-2 outbreak) (see [2]).

Section 5 is structured into three subsection: First, we introduce the data and the methods to approximate the parameters. Then, we propose a set of simulations adapted to the general French MSM population. The last subsection is focused on the so-called high-risk individuals as described previously.

It is important at this stage to remind that HIV is not curable and that PrEP does not have a direct influence on the decreasing number of the infected compartment over time, which is mainly due to natural death (with a life expectancy for treated persons as high as non-infected ones) or, more rarely now, thus neglected here for the French group, untreated AIDS individuals. However, PrEP policy in the population may have an important impact on the infectivity rate. Since PrEP’s efficiency can be measured in a long-term range and using incidence, we simulated over 15 years (that is, 180 months) for a forecast of this prevention strategy over time.

5.1. Choice of the Parameters

In Table 1, we give the parameters used in Model (1), as well as their meaning and values for the general and high-risk MSM community in France estimated from [2], assuming that $60\%$ of the the population considered in the dataset of the study are MSM.

PrEP treatment in France started in 2016, and since the SARS-CoV-2 epidemic drastically modified sexual habits (through successive lockdowns), we considered the population data until 2019 only. Then, after, we simulated our model based on the estimated parameter values for the following 15 years. Note here that we did not take the SARS-CoV-2 nor the monkey pox epidemic into account for three reasons: first, because these periods are a bias for our model and need to be carefully adjusted with updated parameter values; second, to the best of our knowledge, the official French data related to these past two years have not been released yet; finally, we preferred keeping our model as simple as possible for this present work. A more complex approach will be the object of a future work.

In [2], information between 2016 and 2019 was collected every 6 months. We remind here that, at the end of a 3-month period, the patient decides whether to give up the treatment or to continue. This is the reason why our simulations below were computed on a monthly basis.

Initial conditions for functions $t\mapsto S\left(t\right),t\mapsto I\left(t\right)$ and $t\mapsto u\left(t\right)$ were chosen to be January, the 1st of 2016, according to [2].

Because of the delay equation for u, the function $u_{\mathrm{init}}:[-3,0]\to {\mathbb{R}}^{+},t\mapsto u_{\mathrm{init}}\left(t\right)$ was adjusted with a cubic spline interpolation of 60% of each total number of the last column in Table 2. Note that $u_{\mathrm{init}}$ is the same for both populations’ (general MSM and its high-risk subset) simulations because we assumed here that most of the French PrEP users belong to the high-risk MSM group.

Regarding the susceptible: based on an AVAC study (https://www.avac.org/sites/default/files/u3/MSM_in_Europe_Euro_Rave.pdf, page 5) (accessed on 13 September 2022), 4–10% of French males declare belonging to the MSM group. Assuming that France has about 65,000,000 inhabitants, it seems then reasonable to assume that $S\left(0\right)=$ 2,600,000.

Furthermore, on the official website of UNAIDS (https://www.unaids.org/en/regionscountries/countries/france) (accessed on 13 September 2022), we obtained that 170,000 individuals were infected with HIV in France in 2016. We assumed that a majority belong to the MSM group, but not the entire number (indeed, a large portion of new cases was reported to be foreign heterosexuals, among other subcases). This is why we took $I\left(0\right)=$ 90,000. It is important to remind here that this value may not be the exact one, since the number of HIV infections detected is always below the real number of HIV infections (HIV testing is not compulsory, and thus, several individuals may carry the virus for years without knowing it).

Finally, when simulating the MSM high-risk population, the estimated susceptible number would be between 3000 and 50,000 depending on the literature. Thus, we took $S\left(0\right)=$ 40,000 according to [2]. According to French data (https://vih.org/20190328/stabilite-des-chiffres-du-vih-en-france/) (accessed on 13 September 2022), French infections have been stable since 2010. On average, there are 6500 contaminations per year, and infections within the MSM population represent 43%. Then, we chose $I\left(0\right)=9000$.

According to the official French data (https://www.insee.fr/fr/statistiques/2383440) (accessed on 13 September 2022), $\mu =0.000758333$ individuals.${\mathrm{months}}^{-1}$ at the country scale. We chose this value for the French MSM removal rate (death, as well as removed from sexual life). For the high-risk MSM, $\mu$ is considered much larger because a high-risk sexual life style does not last longer than the average MSM one. Thus, $\mu =0.0076$ individuals.${\mathrm{months}}^{-1}$ in this subgroup.

The 2016–2019 semester official French datasets are given in Table 3.

Thanks to Table 2 and Table 3, we are able to assess the values of functions $\theta$ and $\psi$ for the years 2016–2019. Indeed, for each semester, we obtained a different $\theta$ by using the following formula coming from its definition:

$\theta \left(\mathrm{semester}\right)=\frac{\mathrm{Number}\mathrm{of}\mathrm{renewal}\mathrm{treatment}\mathrm{of}\mathrm{the}\mathrm{current}\mathrm{semester}}{\mathrm{Total}\mathrm{of}\mathrm{PrEP}\mathrm{users}\mathrm{of}\mathrm{the}\mathrm{previous}\mathrm{semester}}.$

We assumed that this value is the same for every month of the semester. Then, we chose $\psi$, per month, as follows:

$\psi \left(\mathrm{months}\right)=\frac{\mathrm{Number}\mathrm{of}\mathrm{individuals}\mathrm{who}\mathrm{begin}\mathrm{the}\mathrm{treatment}}{S\left(0\right)\ast 6},$

where $S\left(0\right)$ is the initial condition for the susceptible compartment S.

These values are given in Table 4.

We selected $\theta =0.83$ as the average value of all the previous $\theta$ from Table 4.

As mentioned in Section 2, we assumed that $\psi$ is the solution of a logistic equation given by:

${\displaystyle \frac{d\psi }{dt}}=r\psi \left(1-\frac{\psi }{K}\right),$

where K is the carrying capacity and r exponential growth. We remind that the explicit form of $\psi$ is then written by the following expression:

$\psi \left(t\right)={\displaystyle \frac{K}{1+\left({\displaystyle \frac{K}{{\psi }_0}}-1\right)exp(-rt)}},$

where ${\psi }_0$ is the initial condition depending on the population type. Using the data for $\psi$ given in Table 4, we summarize our results in Table 5.

We assumed a Hill function for the dynamics of the function f:

$f\left(S\right)=S_{sat}\frac{S^n}{{\gamma }^n+S^n},$

with $S_{sat}$ the saturation of the Hill function, $\gamma$ the abscissa of the inflection point, and n the intensity of the slope. We summarize the values of these parameters in Table 6.

One of our goals was to estimate the HIV epidemic ${\mathcal{R}}_0$ for French MSM. We used the package in the R^® language untitled R0 (https://www.rdocumentation.org/packages/R0/versions/1.2-6) (accessed on 13 September 2022) and, precisely, the function est.R0.SB, which estimates ${\mathcal{R}}_0$ using a Bayesian approach following the idea developed in [13]. Thanks to our data of new HIV infections ([14,15]), we obtained ${\mathcal{R}}_0=0.93$ for the normal MSM French population.

The two remaining parameters, $\beta$ and $\sigma$, were estimated to keep ${\mathcal{R}}_0$ equal to $0.93$. We decided to choose $\sigma =3000$; in other words, each month, 3000 individuals might join the susceptible compartment (by reaching the sexual consent age or deciding to join the MSM group). Thus, $\beta =1.821\times {10}^{-10}{\left(\mathrm{individuals}.\mathrm{months}\right)}^{-1}$ was estimated for the French general MSM group.

On the other hand, according to the official French data (https://solidarites-sante.gouv.fr/IMG/pdf/argumentaire_depistage_vih_HAS_2009-2.pdf) (accessed on 13 September 2022), the basic reproduction number related to the high-risk MSM ranges between $1.27$ to $1.44$. We took ${\mathcal{R}}_0=1.31$ and estimated $\beta =2.85\times {10}^{-7}{\left(\mathrm{individuals}.\mathrm{months}\right)}^{-1}$, with $\mu =0.0076$ ${\mathrm{months}}^{-1}$ and $\sigma =562$ given by Table 1.

5.2. Numerical Simulations for the General French MSM

In Figure 2, protected (P) in green (graph on the right) keeps track of the exact data (crosses) for the first semesters, then reaches a threshold due to the French policy of regulating PrEP users. Indeed, the treatment is fully taken care of by the health insurance in France. Thus, a threshold needs to be set up (estimated at 60,000 MSM individuals here). With this threshold, it seems already that, within 15 years, the number of infected drops and keeps on decreasing. The number of susceptibles (blue curve on the left of the graph) keeps on growing with a lower slope between 50 and 100 months due to the increase of the protected and then the threshold. However, even with the increase of the susceptibles, the infected remain low. The protected compartment plays the role of a reservoir to prevent the HIV epidemic from growing back.

Just as a comparison point, we simulated our model with a function f defined as identity (and not a Hill function anymore) (see Figure 3). We easily observed that the graph of the protected population does not fit the data in the first months with the same realistic parameter sets. Then, it seems that f needs to be more complex than identity for the case of French MSM, and the Hill function seems to be validated here.

In Figure 4, we plot the incidence with or without PrEP treatment. It seems obvious that without the PrEP reservoir (in blue), the incidence keeps on growing with no chance for the HIV epidemic to be controlled. On the other hand, with the PrEP compartment (in red), it is possible to keep the incidence at a low level (and even to decrease it significantly for the first years). The slight increase at the end (after Month 125) is due to the threshold of the number of PrEP users imposed in our model and likely used by the French health insurance policy.

In Figure 5, we depict the effect of the variation of ${\mathcal{R}}_0$ as a function of $\psi$. It appears clearly here that, as $\psi$ increases (that is, the number of PrEP users rises), ${\mathcal{R}}_0$ declines drastically. This is, however, not linear, and a plateau may be reached after a certain threshold (above $0.2$), which indicates that the flow of new PrEP users does not need to expand drastically. Indeed, even at $\psi$ larger than $0.125$, we obtain ${\mathcal{R}}_0$ already below $0.3$, which is quite satisfactory. Augmenting $\psi$ would decrease ${\mathcal{R}}_0$, but not as fast as the first values of $\psi$.

5.3. Numerical Simulations for High-Risk Population

In this subsection, we focus our attention to the French MSM high-risk population. We remind that this MSM subgroup consists of individuals with multi-partner intercourse, sex relationships while using drugs (chemsex), etc. The HIV risk of infection is thus much higher than the rest of the MSM population. This is one of the reasons why this subgroup is the PrEP treatment target for the French health insurance. The protected compartment in this particular case plays a major role as a reservoir. This was clearly confirmed in our model simulations, as shown in Figure 6, where the number of susceptibles clearly drops as the number of protected individuals increases. Note that the infected population rises, but for the first time, reaches a threshold. This means that barely any new infected cases appear. We remind here that the model is an $SI$ and not an $SIR$, that is no recovered individuals can appear since $HIV$ is a disease with a lifetime treatment.

In Figure 7, the incidence evolution is plotted with or without PrEP treatment. The reservoir effect is even more explicit there in the sense that, with PrEP, the incidence declines to reach a low plateau, while without, it keeps on growing indefinitely.

Finally, as mentioned before, for the MSM high-risk population, the ${\mathcal{R}}_0$ in France is higher than 1. Figure 8 reveals that, as the number $\psi$ of new PrEP users increases, ${\mathcal{R}}_0$ not only decreases below one, but also keeps on falling off. Here, the critical value of $\psi$, which keeps ${\mathcal{R}}_0$ below one, is ${\psi }_{critical}=0.113\mathrm{person}.{\mathrm{months}}^{-1}$. According to our data, $\psi$ = 0.071 $\mathrm{person}.{\mathrm{months}}^{-1}$, as we assume the maximum threshold has been reached. The decay appears, however with a smaller slope, which means, as in the previous subsection, that, after a certain threshold ($0.2$ here), very large values of $\psi$ do not influence the reduction of ${\mathcal{R}}_0$ much.

This reservoir effect is so important that, if for any reason, there is the decision to stop PrEP treatment for this MSM group, the rise in new infected individuals would be inevitable, as shown in Figure 9.

6. Conclusions and Discussion

The aim of this work was to propose a brand-new model for prescribing PrEP for 3 months, with the choice to continue treatment or to stop it at any time after this period. This choice changes many things in the new modeling approach by introducing difference equations with delay and, thus, the possibility of describing the effect of inertia and the reservoir of protected individuals, not only analytically, but also numerically. This approach could also be generalized to other epidemiological models describing the effect of vaccine or booster doses for improved efficacy of protection.

Furthermore, the double insights of the MSM general and high-risk populations seemed very important to us, especially to depict the effect of PrEP treatment on the disease incidence of the disease and the French health insurance policy to target the high-risk subgroup. With our model and simulations, it appeared that these strategies are much more effective than random targeting of any MSM individual.

Our future work, in production, will take into account that the exponential growth of susceptibles is not realistic for a longer period of time, so it would be better to model it as a logistic growth and then numerically simulated on a longer period of time.

The data we worked on were from 2016 to 2019. It would be interesting to analyze the effect of the SARS-CoV-2 lockdown effect on the incidence and PrEP users. It would also be informative to analyze the booster effect of the non lockdown period (mid-2021 and after). For this purpose, it is necessary to wait a few more years to collect a sufficient amount of data.

A future model should also consider $\theta$ time dependent. Indeed, depending on certain periods (a lockdown, for example, or a seasonal effect), the number of PrEP users willing to take the treatment at the end of their 3-month prescription may vary over time. Some users may also take it for a short period of time (e.g., 2 years on average) and discontinue treatment after choosing a long-term closed relationship.

The addition of a compartment of the infected population under tritherapy should also be more realistic. Indeed, with tritherapy, HIV individuals can no longer be infectious and can be removed from the I compartment. The estimate of $\beta$ would then be more accurate and reflect a more realistic situation.

Finally, it would be interesting to add the cost of treatment and the cost-efficient strategy with the threshold of PrEP users in the compartment for optimal disease control and, eventually, its disappearance. For this purpose, our next goal is to compare data and policies across multiple countries (some with full financial support for treatment, some with partial support, and some with no support).

Author Contributions

Supervision, M.A., L.P.-M. and J.W.; Writing—original draft, J.M.; Writing—review & editing, G.R. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Schematic representation of the compartmental model (1). Continuous arrows represent movements between compartments. Dashed ones represents the transmission of the infection.

Figure 2. Plot of the evolution of the French MSM population (4) along the time (over 15 years). The crosses in the last plot represent the real values of PrEP users taken from Table 3. Function [Forumla omitted. See PDF.] verifies the logistic equation, and f is a Hill function.

Figure 3. Plot of the evolution of the different compartments of the model (4) along time (over 15 years). The crosses in the last plot represent the real values of PrEP users taken from Table 3. [Forumla omitted. See PDF.] verifies the logistic equation, and f is identity.

Figure 4. Plot of the evolution of the incidence with PrEP (in red) and without (in blue) from the system (4) over 15 years among the general French MSM.

Figure 5. Plot of the [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.].

Figure 6. Plot of the evolution of the different compartments of the high-risk french population along time (over 15 years). The crosses in the last plot represent the real values of the number of PrEP users obtained from Table 4. [Forumla omitted. See PDF.] verifies the logistic equation, and f is a Hill function.

Figure 7. Plot of the evolution of the incidence with and without PrEP (4) over 15 years among French high-risk MSM.

Figure 8. Plot of the [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.].

Figure 9. Plot of the susceptibles and infected if PrEP were stopped in 15 years. We considered that all the protected individuals become susceptibles again.

Table 1

Parameters used in Model (1) with values applied to the French population estimated from [2].

Symbols	Signification	Value for General MSM	Value for High-Risk MSM	Unity
$\theta$	Probability to keep the PrEP treatment	$0.83$	$0.83$
$\sigma$	Recruitment	$3000$	562	indiv.months${}^{-1}$
$\mu$	Removal rate from the compartments	$0.000758333$	$0.0076$	indiv.months${}^{-1}$
$\beta$	HIV transmission rate per infected individual	$1.821\times {10}^{-10}$	$2.85\times {10}^{-7}$	(indiv.months)${}^{-1}$
$\tau$	Duration of the period of PrEP taking	3	3	months

Table 2

Summary of initial conditions for susceptibles and infected among general MSM and MSM high-risk populations.

Initial Conditions	Value for MSM	Value for High-Risk MSM
$S\left(0\right)$	2,600,000	40,000
$I\left(0\right)$	90,000	9000

Table 3

Total number of PrEP users in France since 2016, given by semester (see Table 3 in [2] ).

Semester	Initiation of PrEP	Renewal of the Treatment	Total PrEP Users
S1—2016	1166	###	1166
S2—2016	1826	911	2737
S1—2017	2193	2273	4466
S2—2017	2564	3807	6371
S1—2018	3138	5413	8551
S2—2018	4488	7647	12,135
S1—2019	5103	10,398	15,501

Table 4

Values of parameters $\psi$ and $\theta$ per month according to each semester, computed according to Table 2 and Table 3.

Semester	Values of $\mathit{\psi }$ for General MSM	Values of $\mathit{\psi }$ for High-Risk MSMS	Values of $\mathit{\theta }$
S1—2016	$0.0000747$	$0.0048583$	###
S2—2016	$0.000117$	$0.007608$	$0.7813$
S1—2017	$0.00014$	$0.00913$	$0.8305$
S2—2017	$0.000164$	$0.0106$	$0.8159$
S1—2018	$0.000201$	$0.0130$	$0.8496$
S2—2018	$0.00029$	$0.0187$	$0.8943$
S1—2019	$0.00033$	$0.0216$	$0.8569$

Table 5

Values of parameters K and r depending on the population.

Parameters	Values for General MSM	Values for High-Risk MSM
K	$0.0007$	$0.072$
r	$0.0000222$	$0.000026$

Table 6

Values of parameters $S_{sat}$, $\gamma$, and n depending on the general French MSM and the high-risk French MSM.

Parameters	Values for MSM	Values for High-Risk MSM
$S_{sat}$	$5\times {10}^6$	230,000
$\gamma$	120	50
n	2	$1.56$

References

1. Van Manen, D.; Delaneau, O.; Koostra, A.; Boeser-Nunnink, B.; Limou, S.; Bol, S.; Burger, J.; Zwinderman, A.; Moerland, P.; van’t Slot, R. et al. Genome-Wide Association Scan in HIV-1-Infected Individuals Identifying Variants Influencing Disease Course. PLoS ONE; 2011; 6, e22208. [DOI: https://dx.doi.org/10.1371/journal.pone.0022208] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21811574]

2. Billioti de Gage, S.; Le Tri, T.; Dray-Spira, R. Suivi de l’utilisation de Truvada® ou Génériques Pour une Prophylaxie Pré-exposition (PrEP) au VIH à Partir des Données du Système National des Données de Santé (SNDS)—Actualisation des Données Jusqu’au 30 Juin 2019 (27/11/2019). 2019; Available online: https://www.epi-phare.fr/app/uploads/2020/04/Billioti-de-Gage-et-al.-2019-Suivi-de-l-%E2%80%99utilisation-de-Truvada%C2%AE-ou-ge%CC%81ne%CC%81riques-p.pdf (accessed on 13 September 2022).

3. Silva, C.J.; Torres, D.F.M. A SICA compartmental model in epidemiology with application to HIV/AIDS in Cape Verde. Ecol. Complex.; 2017; 30, pp. 70-75. [DOI: https://dx.doi.org/10.1016/j.ecocom.2016.12.001]

4. Akpa, O.; Oyejola, B. Modeling the transmission dynamics of HIV/AIDS epidemics: An introduction and a review. J. Infect. Dev. Ctries.; 2010; 4, pp. 597-608. [DOI: https://dx.doi.org/10.3855/jidc.542] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21045351]

5. Rivadeneira, P.; Moog, C.; Stan, G.B.; Brunet, C.; Raffi, F.; Ferre, V.; Costanza, V.; Mhawej, M.J.; Biafore, F.; Ouattara, D.A. et al. Mathematical Modeling of HIV Dynamics After Antiretroviral Therapy Initiation: A Review. BioRes. Open Access; 2014; 3, pp. 233-241. [DOI: https://dx.doi.org/10.1089/biores.2014.0024] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25371860]

6. Liu, D.; Wang, B. A novel time delayed HIV/AIDS model with vaccination & antiretroviral therapy and its stability analysis. Appl. Math. Model.; 2013; 37, pp. 4608-4625. [DOI: https://dx.doi.org/10.1016/j.apm.2012.09.065]

7. Nsuami, M.U.; Witbooi, P.J. A model of HIV/AIDS population dynamics including ARV treatment and pre-exposure prophylaxis. Adv. Differ. Equ.; 2018; 2018, pp. 1-12. [DOI: https://dx.doi.org/10.1186/s13662-017-1458-x]

8. Silva, C.J.; Torres, D.F.M. Modeling and optimal control of HIV/AIDS prevention through PrEP. Discret. Contin. Dyn. Syst.; 2018; 11, pp. 119-141. [DOI: https://dx.doi.org/10.3934/dcdss.2018008]

9. Simpson, L.; Gumel, A.B. Mathematical assessment of the role of pre-exposure prophylaxis on HIV transmission dynamics. Appl. Math. Comput.; 2017; 293, pp. 168-193. [DOI: https://dx.doi.org/10.1016/j.amc.2016.07.043]

10. Pinto, C.M.A.; Carvalho, A.R.M. The impact of pre-exposure prophylaxis (PrEP) and screening on the dynamics of HIV. J. Comput. Appl. Math.; 2018; 339, pp. 231-244. [DOI: https://dx.doi.org/10.1016/j.cam.2017.10.019]

11. Adimy, M.; Chekroun, A.; Ferreira, C. Global dynamics of a differential–difference system: A case of Kermack-McKendrick SIR model with age-structured protection phase. Math. Biosci. Eng.; 2020; 17, pp. 1329-1354. [DOI: https://dx.doi.org/10.3934/mbe.2020067] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32233581]

12. Webb, G.F. Theory of Nonlinear Age-Dependent Population Dynamics; M. Dekker New York: CRC Press: Boca Raton, FL, USA, 1985; 294.

13. Bettencourt, L.; Ribeiro, R. Real Time Bayesian Estimation of the Epidemic Potential of Emerging Infectious Diseases. PLoS ONE; 2008; 3, e2185. [DOI: https://dx.doi.org/10.1371/journal.pone.0002185] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18478118]

14. Cazein, F.; Pillonel, J.; Le Strat, Y.; Pinget, R.; Le Vu, S.; Brunet, S.; Thierry, D.; Brand, D.; Leclerc, M.; Benyelles, L. et al. Découvertes de séropositivité VIH et sida—France, 2003–2013. Bull. Epidemiol. Hebd.; 2015; 9–10, pp. 152-161.

15. Cazein, F.; Pillonel, J.; Sommen, C.; Bruyand, M.; Lydié, N.; Che, D.; Coignard, B.; Lot, F. Bulletin de Santé Publique. Découvertes de séropositivité VIH et diagnostique de SIDA—France, 2018. Bull. Santé Publique Vih/Sida; 2019.