# Stability Analysis of a New Differential-Difference Model Applied to the Pre-exposure Prophylaxis (PrEP) Effect on the Spread of HIV

Mostafa Adimy, Abdennasser Chekroun, Grégoire Ranson, Laurent Pujo-Menjouet

► To cite this version:

Mostafa Adimy, Abdennasser Chekroun, Grégoire Ranson, Laurent Pujo-Menjouet. Stability Analysis of a New Differential-Difference Model Applied to the Pre-exposure Prophylaxis (PrEP) Effect on the Spread of HIV. Qualitative Theory of Dynamical Systems, 2024, 23 (5 (238)), 10.1007/s12346-024-01093-x. hal-04777788

HAL Id: hal-04777788

https://inria.hal.science/hal-04777788v1

Submitted on 12 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Stability analysis of a new differential-difference
model applied to the Pre-Exposure Prophylaxis
(PrEP) effect on the spread of HIV

Mostafa Adimy<fnref n="1" />*†, Abdennasser Chekroun<fnref n="2" />†, Grégoire Ranson<fnref n="3" />,<fnref n="4" />†,
Laurent Pujo Menjouet<fnref n="3" />†

Abstract

In this research, an examination is conducted on a model derived from the inquiry
into the efficacy of HIV Pre-Exposure Prophylaxis (PrEP) within high-risk pop-
ulations. To achieve this objective, we employ an SI model coupled with an
age-structured equation to delineate the dynamics of individuals under PrEP pro-
tection, incorporating an infected-dependent rate of new users recruited from the
susceptible population. This nonlinear term is contingent upon a factor dedicated
to the political or economic context of a government. Local asymptotic stability
for both disease-free and endemic equilibria is established, and global asymp-
totic stability for the disease-free steady-state is demonstrated. To address the

<footnote n="1">¹*Inria, Universite Claude Bernard Lyon 1, CNRS, Ecole Centrale de Lyon, INSA Lyon, Universite Jean Monnet, ICJ UMR5208, Villeurbanne, F-69200, France.</footnote>

<footnote n="2">²Laboratoire d'Analyse Non linéaire et Mathématiques Appliquées,
University of Tlemcen, Tlemcen, 13000, Algeria.</footnote>

<footnote n="3">³Universite Claude Bernard Lyon 1, CNRS, Ecole Centrale de Lyon, INSA Lyon, Université Jean Monnet, ICJ UMR5208, Inria, Villeurbanne, F-69622, France.</footnote>

<footnote n="4">⁴Laboratory for Industrial and Applied Mathematics (LIAM),
Department of Mathematics and Statistics, York University, Toronto,
ON M3J 1P3, Canada.</footnote>

<footnote>*Corresponding author(s). E-mail(s): mostafa.adimy@inria.fr;</footnote>

<footnote>Contributing authors: chekroun@math.univ-lyon1.fr; granson@yorku.ca;
pujo@math.univ-lyon1.fr;</footnote>

<footnote>†These authors contributed equally to this work.</footnote>
system's behavior, a reduction of the partial differential equation is undertaken, presenting it as a coupled system of differential equations and a delayed difference equation. Lastly, persistence is substantiated when the endemic equilibrium is realized.

**Keywords:** Age-structured PDE, Delay differential-difference system, Local and global asymptotic stability, SI epidemic model.

**MSC Classification:** 34K60, 37N25, 92D30

# 1 Introduction

HIV (Human Immunodeficiency Virus) is a retrovirus that infects humans and leads to AIDS (Acquired Immunodeficiency Syndrome). AIDS marks the advanced stage of the infection, wherein the virus has depleted a substantial count of CD4 T-cells, leading to immune system failure. It is crucial to note that this stage can be fatal if left untreated. Typically, this stage manifests between 5 and 10 years after initial infection, with an average survival time of 9.2 months [1]. At present, there is no cure or vaccine available for HIV, and treatment options become limited once AIDS has developed. In 2019, UNAIDS reported that around 38 million people worldwide were living with HIV. Within that year, approximately 1.7 million individuals contracted the virus, while 690,000 deaths resulted from AIDS-related illnesses [2]. HIV can be transmitted through various means, including unprotected sex, mother-to-child transmission, and exposure to contaminated blood through activities like drug injections. While prevention methods may seem straightforward in theory, implementing them systematically in practice proves more challenging for several reasons. One key factor is that younger generations are less inclined to consistently use condoms during sexual activity [3] while older generations [4] tend to underestimate their risk of infection.

A number of initiatives are already in place to help combat the current HIV epidemic. Among these, two treatments appear to be effective in preventing the spread of the epidemic. Namely, PrEP (Pre-Exposure Prophylaxis) and TasP (Treatment as Prevention).

PrEP is a treatment taken to reduce the risk of HIV infection. It is used by people who are not HIV-positive but who are at high risk of exposure to HIV through sexual intercourse or injecting drug use. PrEP combines several antiretroviral drugs and is available in both oral and injectable forms. The oral form can be taken according to continuous (daily) or discrete (on-demand) protocols. The on-demand regimen involves taking the drug approximately 2 hours before and up to 2 days after each sexual encounter. This worldwide approach has proven highly effective, offering excellent protection against HIV transmission [5]. PrEP users have a significantly reduced risk of HIV infection for three months. However, as PrEP does not protect against other Sexually Transmitted Infections (STIs) [6], they are required to undergo STI screening after this period, to renew their prescription of PrEP. They can then choose to continue
treatment if they wish. When correctly implemented, PrEP has a proven impact on HIV transmission [7].

TasP, on the other hand, targets people who are already infected (HIV-positive). Individuals taking this drug are then unable to transmit the disease, thus protecting their sexual partners [8, 9]. TasP is able to significantly reduce the viral load to undetectable levels and maintain a nearly unchanged lifespan [10]. TasP, when implemented at the population level, can effectively reduce the rate of transmission and bring the epidemic under control.

PrEP and TasP have in common the need for total adherence to treatment to fully exploit their potential to prevent the spread of HIV. However, many factors can lead to treatment failure, such as forgetting to take the medication, stigma-related problems, and side effects [11–13].

Many mathematical models have been investigated with the intent to bring light to difficult questions arising from HIV/AIDS conditions [14, 15]. For a good review, we suggest [16–18]. Delay differential equations have been one of the main tools in HIV immunology modeling over the years [19, 20]. In [21], researchers examined the influence of spatial factors on the MSM (Men who have Sex with Men) population, segmented into immunity-based compartments. They integrated this approach with an optimal control case to investigate the effects of varying PrEP costs. Similarly, [22] explored a combined approach involving both ART (Anti-Retroviral Therapy) and PrEP. They proposed an age-structured model incorporating infection ages and conducted qualitative analysis alongside numerical simulations utilizing data from China.

In [23], the authors developed a model for HIV/AIDS to evaluate the impact of combining PrEP and ART drugs. They conducted mathematical and numerical analyses to explore the effects of these combined controls in the context of PrEP drug resistance. The findings indicated a significant reduction in new HIV infections when both PrEP and ART are applied at high levels. However, the study also showed that resistance to PrEP drugs could diminish or negate PrEP's effects, especially at lower efficacy levels.

In [24], a mathematical model was proposed to understand HIV/AIDS transmission dynamics within a worker population, divided into isolated and non-isolated groups. Isolated workers were HIV-negative individuals who disclosed their status and committed to maintaining it, while non-isolated workers included those either HIV-positive or not wishing to reveal their status. The model extended the work of [23] by considering PrEP use and effectiveness. The study aimed to assess the impact of isolated groups on HIV/AIDS dynamics, revealing that high individual awareness and PrEP efficacy significantly reduce the epidemic.

In [25], the authors designed and calibrated a compartmental model using Canadian data to observe changes in serosorting patterns due to PrEP. The model underscored the critical importance of maintaining high PrEP effectiveness.

Finally, in [26], the authors investigated the issue of reduced PrEP efficacy, primarily due to non-adherence, and its effect on HIV/AIDS infection dynamics. The study aimed to determine the critical threshold of PrEP efficacy needed to prevent backward bifurcation, which occurs when efficacy falls below 100%. Conversely, with
100% efficacy, indicating full adherence, the model exhibited a forward bifurcation. The numerical analysis highlighted the crucial role of adherence in enhancing PrEP efficacy.

Here, we base our model on the classic Kermack-McKendrick SI model [27], integrating an age-structured phase to depict the impact of PrEP only. We do not deal with other STIs or other treatments. These latter will be the object of further work. The age variable is exclusive to the PrEP user compartment, symbolizing the duration spent between prescription renewals. To the best of our knowledge, this is the first model that includes the time between two consecutive prescriptions. Individuals who remain protected during this period are described by an age-structured equation. The age structure allows for a more comprehensive system by adding the age variable alongside the time variable (in our case, the protection age). This system extends a previously published model outlined in [28]. As detailed in [28], it can be reformulated into a delay differential-difference system, where the delay corresponds to the duration of HIV protection provided by PrEP, approximately three months between two prescriptions. The primary alteration concerns the nonlinear force of infection [29] among PrEP users recruited from susceptible individuals. While previously considered as a constant parameter [30], we argue that it should depend not only on the infected population (increasing with the number of infected individuals) but also on factors associated with the political or economic context of a government. Additionally, we propose that only a fraction of the population can afford to utilize PrEP due to economic and political constraints, which we will now elucidate.

On one hand, a country might opt to limit the number of PrEP users, either due to the treatment's cost or because the prevalence of infection is deemed too low to justify widespread use. On the other hand, some governments may choose to fully subsidize the treatment and actively promote its utilization among high-risk individuals. Conversely, other governments might allow individuals to decide whether to pursue PrEP at their own expense or opt not to include any PrEP options in their healthcare plans [31]. Our study aims to examine the behavior of solutions with a rate of new PrEP users, recruited from the susceptible population, assumed to model a government's political or economic context. We are investigating whether, under specific conditions, the endemic equilibrium of this system may become unstable. Could a specific series of political and economic measures lead to a destabilization of this equilibrium, allowing periodic oscillations to occur? We meticulously analyze the stability of each steady state, achieving local asymptotic stability for both disease-free and endemic equilibria. Furthermore, we demonstrate global asymptotic stability for the disease-free steady state. Finally, we establish persistence when the endemic equilibrium exists. Although our study focuses on a simplified model, we believe the results hold validity for the PDE system. Given certain compatibility conditions on initial conditions, equivalence between the two systems can be achieved.

The organization of this paper is as follows. In the next section (Section 2), we formulate the age-structured model by introducing a protection phase by PrEP. We carefully detail the role of each term and parameter. Section 3 focuses on establishing the existence of two steady-states. A threshold basic reproduction number is also
defined. In Section 4, we provide stability of the disease and disease-free equilibria: local, global, and persistence. We conclude with a discussion of our results in Section 5.

## 2 Age-structured model with fixed maximum age

Let us first remind the model introduced in [28]. We consider a population with three categories of individuals: susceptible $S(t)$, infected $I(t)$, and PrEP users $p(t, a)$, where $a \in [0, \tau)$ represents the time since an individual started a new prescription (usually $\tau = 3$ months is the recommended duration between two prescription renewals). Then, the total protected population of all ages is given by

$$P(t) = \int_{0}^{\tau} p(t, a) da,$$

and the entire population studied by

$$N(t) := S(t) + I(t) + P(t).$$

In order to generalize a previous work, we define the recruitment rate of starting PrEP as a function depending on infected $I$ and we note $\varphi(I)$.

We assume that a large number of infected people in a population will push susceptible people to protect themselves. Then, $I \in \mathbb{R}^{+} \mapsto \psi(I)$ is an increasing function. Consequently, we assume $\psi(0) = 0$. If there is no infection there is no individual seeking PrEP.

For political and economic reasons, we consider that only a fraction of the population of susceptible is able to start taking PrEP and we note $f(S)$, where $f$ is an increasing function.

Then, the part of the population starting PrEP is represented by $\psi(I)f(S)$. This term allows us to consider different scenarios. We are giving interpretations in 2, 3 and 4.

At the end of the prescription period $\tau$, only a proportion $0 < \theta < 1$ of PrEP users decide to renew their treatment and remain protected. Thus, $(1 - \theta)p(t, \tau)$ represents the proportion of individuals who decide to stop renewing the prescription and lose their protection. They re-enter the susceptible phase. We denote $\mu$ as the constant death rate common to the three categories and $\Lambda$ be the susceptible recruitment rate. We neglect here vertical (mother/child) transmission. We denote by $\beta$ the rate at which the susceptibles become infected. Now and using the previous description (see Figure 1), we are ready to set up the whole system for $t > 0$ as follows,

$$ 
\left\{ 
\begin{aligned}
\frac{dS(t)}{dt} &= \Lambda - \mu S(t) - \beta I(t)S(t) - \psi(I(t))f(S(t)) + (1 - \theta)p(t, \tau), \\
\frac{dI(t)}{dt} &= \beta I(t)S(t) - \mu I(t), \\
\frac{\partial p}{\partial t}(t, a) + \frac{\partial p}{\partial a}(t, a) &= -\mu p(t, a), && 0 < a < \tau, \\
p(t, 0) &= \psi(I(t))f(S(t)) + \theta p(t, \tau),
\end{aligned}
\right. 
\quad (1) $$
with initial conditions

$$S(0) = S_0, \quad I(0) = I_0 \quad \text{and} \quad p(0, a) = p_0(a), \quad 0 < a < \tau.$$

<figure><img src="image_2.png" /><figcaption>Fig. 1 Schematic representation of the process. The population splits into susceptible $S(t)$, infected $I(t)$, and protected individuals under PrEP $p(t, a)$, where $a$ stands for the elapsed time (finite time: less than $\tau$, between to prescription renewals) in the PrEP compartment.</figcaption></figure>

We can reduce the model by using the method of characteristics (see [28, 30, 32] for the details) and introducing the function $\phi$ defined by

$$\phi(t) = e^{-\mu t} p_0(-t), \quad t \in [-\tau, 0],$$

and putting the new state

$$u(t) := p(t, 0), \quad t > 0,$$

it is possible to obtain

$$p(t, \tau) = e^{-\mu \tau} \begin{cases} u(t - \tau), & t > \tau, \\ \phi(t - \tau), & t \in [0, \tau]. \end{cases}$$
Then, we have the following coupled differential-delayed difference equations, for $t > 0$,

$$ egin{cases} 
\frac{dS(t)}{dt} = \Lambda - \mu S(t) - \beta S(t)I(t) - \psi(I(t))f(S(t)) + (1-\theta)e^{-\mu\tau}u(t-\tau), \\ 
\frac{dI(t)}{dt} = -\mu I(t) + \beta S(t)I(t), \\ 
u(t) = \psi(I(t))f(S(t)) + \theta e^{-\mu\tau}u(t-\tau), 
\end{cases} 
\quad (2) $$

which is completed by the initial conditions

$$ S(0) = S_0, \quad I(0) = I_0 \quad \text{and} \quad u(t) = \phi(t), \quad \text{for} \quad -\tau \le t \le 0. $$

The existence, uniqueness, and positivity of the solution of this system may be deduced in a classic way [30] or [32].

However, in the context of our objective to generalize the study of [28] where only the homogeneous steady-states were investigated, we decided to analyze the dynamics of the age-distributed equilibria from the original age-structured model (1). Moreover, both systems can be linked by the same states $S$, $I$ and considering, for $a \in (0, \tau)$,

$$ p(t, a) = e^{-\mu a} \begin{cases} \phi(t-a), & 0 < t \le a, \\ u(t-a), & t > a. \end{cases} $$

### 3 Existence of steady-states

Beforehand, observe first that the equation on $p(t, a)$ can be simplified by integrating over the age, for $t > 0$,

$$ \frac{dP(t)}{dt} = -\mu P(t) - p(t, \tau) + p(t, 0) = -\mu P(t) + \psi(I(t))f(S(t)) - (1-\theta)p(t, \tau). $$

Then, by adding it to the two first equations of the system (1), we obtain, for $t > 0$,

$$ \frac{dN(t)}{dt} = \frac{dS(t)}{dt} + \frac{dI(t)}{dt} + \frac{dP(t)}{dt} = \Lambda - \mu(S(t) + I(t) + P(t)) = \Lambda - \mu N(t). $$

This means that

$$ \lim_{t \to +\infty} N(t) = \frac{\Lambda}{\mu}. \quad (3) $$
Then, we denote the steady-states by $(S^*, I^*, p^*(a))$, $0 < a < \tau$. They satisfy the following system

$$ egin{cases} 
\Lambda - \mu S^* - \beta I^* S^* - \psi(I^*) f(S^*) + (1-\theta)p^*(\tau) = 0, \\
\beta I^* S^* - \mu I^* = 0, \\
\frac{d}{da} p^*(a) = -\mu p^*(a), \\
p^*(0) = \psi(I^*) f(S^*) + \theta p^*(\tau). 
\end{cases} 
\quad (4) $$

The equation of $p^*$ has an analytical expression given by

$$ p^*(a) = e^{-\mu a} p^*(0). $$

Two possible steady-states may arise then: the disease-free equilibrium (DFE) characterized by the fact that $I^* = 0$, and the endemic one for which $I^* > 0$.

Let us investigate the disease-free steady-state. Considering $I^* = 0$ yields to

$$ egin{cases} 
\Lambda - \mu S^* + (1-\theta)e^{-\mu\tau}p^*(0) = 0, \\
(1-\theta)e^{-\mu\tau})p^*(0) = 0. 
\end{cases} $$

Since $\theta < 1 < e^{\mu\tau}$, then $1 - \theta e^{-\mu\tau} > 0$ and the second equation of the system implies that

$$ p^*(0) = 0 \quad \text{and} \quad S^* = \frac{\Lambda}{\mu}. $$

It follows that there is a unique DFE and is given by

$$ \mathcal{E}_0 = \left( \frac{\Lambda}{\mu}, 0, 0 \right). $$

As previously mentioned, for an endemic steady-state, we consider $I^* > 0$. From the second equation of (4) we have

$$ S^* = \frac{\mu}{\beta}. $$

Substituting the expression of $S^*$ in (4) leads to

$$ egin{cases} 
\Lambda - \frac{\mu^2}{\beta} - \mu I^* - \psi(I^*) f\left(\frac{\mu}{\beta}\right) + (1-\theta)e^{-\mu\tau}p^*(0) = 0, \\
p^*(a) = e^{-\mu a}p^*(0), \\
p^*(0) = \frac{1}{1-\theta e^{-\mu\tau}}\psi(I^*)f\left(\frac{\mu}{\beta}\right). 
\end{cases} $$

Then, by replacing the expression of $p^*(0)$ into the first equation, we obtain

$$ \Lambda - \frac{\mu^2}{\beta} - \mu I^* - \left[ 1 - \frac{(1-\theta)e^{-\mu\tau}}{1-\theta e^{-\mu\tau}} \right] \psi(I^*) f\left(\frac{\mu}{\beta}\right) = 0. $$
To prove the existence of a solution to this equation, we consider the function $\mathcal{H}$ defined on $\mathbb{R}^+$ by

$$ \mathcal{H}(I^*) := \frac{\beta\Lambda}{\mu^2 + \beta\mu I^* + \left[ \frac{\beta(1 - e^{-\mu\tau})}{1 - \theta e^{-\mu\tau}} \right] \psi(I^*) f\left(\frac{\mu}{\beta}\right)}. $$

The existence of an endemic steady-state is then equivalent to solving the equation

$$ \mathcal{H}(I^*) = 1. $$

Since all the parameters of the system are positive, we easily check that $\mathcal{H}$ is a decreasing function on $\mathbb{R}^+$ with $\lim_{I^* \to +\infty} \mathcal{H}(I^*) = 0$. We are then able to define the basic reproduction number $\mathcal{R}_0$ as

$$ \mathcal{R}_0 := \mathcal{H}(0) = \frac{\beta\Lambda}{\mu^2}. $$

Therefore, we conclude that

* If $\mathcal{R}_0 > 1$, then there is a unique endemic steady-state.

* If $\mathcal{R}_0 < 1$, then the maximal value of the function $\mathcal{H}$ is $\mathcal{R}_0$ and there is no endemic steady-state.

We summarize the existence of equilibria in the following theorem.

**Theorem 1.** If $\mathcal{R}_0 > 1$, then the system (1) has two distinct steady-states: a disease-free steady-state $\mathcal{E}_0$ and an endemic steady-state $\mathcal{E}^* := (S^*, I^*, p^*(a))$. If $\mathcal{R}_0 < 1$, then $\mathcal{E}_0$ is the only steady-state.

Based on the equilibrium $\mathcal{E}^* := (S^*, I^*, p^*(a))$ of System (1) where

$$ \left\{ \begin{array}{l} S^* = \frac{\mu}{\beta}, \\ p^*(a) = e^{-\mu a} p^*(0), \\ p^*(0) = \frac{1}{1 - \theta e^{-\mu \tau}} \psi(I^*) f\left(\frac{\mu}{\beta}\right), \end{array} \right. $$

we notice that the positive equilibrium of System (2) is given by

$$ (S^*, I^*, u^*) = (S^*, I^*, p^*(0)). $$
## 4 Stability of the two steady-states and persistence analysis

To study the stability of the two steady-states, we linearize the system (1) about each steady-state $(S^*, I^*, p^*(a))$. We put, for $t > 0$,

$$ egin{cases} S(t) = S^* + s(t), \\ I(t) = I^* + i(t), \\ p(t, a) = p^*(a) + q(t, a). \end{cases} $$

Then, by deriving these expressions, the system (1) yields to

$$ egin{cases} \frac{dS(t)}{dt} = \Lambda - \mu S(t) - \beta I(t)S(t) - \psi(I(t))f(S(t)) + (1 - \theta)p(t, \tau), \\ \frac{dI(t)}{dt} = \beta I(t)S(t) - \mu I(t), \\ \frac{\partial q}{\partial t}(t, a) + \frac{\partial q}{\partial a}(t, a) + \frac{d}{da}p^*(a) = -\mu p(t, a), \quad 0 < a < \tau, \\ q(t, 0) + p^*(0) = \psi(I(t))f(S(t)) + \theta p(t, \tau). \end{cases} $$

We replace the functions $S$, $I$ and $p$ by their expressions

$$ egin{cases} \frac{dS(t)}{dt} = -(\mu + \beta I^* + \psi(I^*) \frac{df(S^*)}{dS})s(t) - (\beta S^* + \frac{d\psi(I^*)}{dI}f(S^*))i(t) \\ + (1 - \theta)q(t, \tau), \\ \frac{dI(t)}{dt} = \beta I^* s(t) + (\beta S^* - \mu) i(t), \\ \frac{\partial q}{\partial t}(t, a) + \frac{\partial q}{\partial a}(t, a) = -\mu q(t, a), \quad 0 < a < \tau, \\ q(t, 0) = \psi(I^*) \frac{df(S^*)}{dS} s(t) + \frac{d\psi(I^*)}{dI} f(S^*) i(t) + \theta q(t, \tau). \end{cases} $$

We look for solutions with exponential growth of the form

$$ s(t) = e^{\lambda t} s_0, \quad i(t) = e^{\lambda t} i_0, \quad \text{and} \quad q(t, a) = e^{\lambda t} q_0(a), \quad \lambda \in \mathbb{C}. $$

We substitute in the previous system and we simplify by $e^{\lambda t}$, to obtain

$$ egin{cases} \lambda s_0 = - \left( \mu + \beta I^* + \psi(I^*) \frac{df(S^*)}{dS} \right) s_0 - \left( \beta S^* + \frac{d\psi(I^*)}{dI} f(S^*) \right) i_0 \\ + (1 - \theta) q_0(\tau), \\ \lambda i_0 = \beta I^* s_0 + (\beta S^* - \mu) i_0, \\ \lambda q_0(a) + q'_0(a) = -\mu q_0(a), \\ q_0(0) = \psi(I^*) \frac{df(S^*)}{dS} s_0 + \frac{d\psi(I^*)}{dI} f(S^*) i_0 + \theta q_0(\tau). \end{cases} $$
The equation on $q_0$ can be solved analytically and we have the exact solution

$$q_0(a) = e^{-(\lambda+\mu)a}q_0(0).$$

This leads to the following linear system

$$\left\{
\begin{array}{l}
\left(\lambda + \mu + \beta I^* + \psi(I^*) \frac{df(S^*)}{dS}\right) s_0 + \left(\beta S^* + \frac{d\psi(I^*)}{dI} f(S^*)\right) i_0 \\
\qquad - (1-\theta)e^{-(\lambda+\mu)\tau} q_0(0) = 0, \\[1.5ex]
-\beta I^* s_0 + (\lambda - \beta S^* + \mu) i_0 = 0, \\[1.5ex]
-\psi(I^*) \frac{df(S^*)}{dS} s_0 - \frac{d\psi(I^*)}{dI} f(S^*) i_0 + (1-\theta)e^{-(\lambda+\mu)\tau} q_0(0) = 0.
\end{array}
\right.$$

This system can be written with a matrix of the form $A(s_0, i_0, q_0(0))^T = 0$, where

$$A = \begin{pmatrix}
\lambda + \mu + \beta I^* + \psi(I^*) \frac{df(S^*)}{dS} & \beta S^* + \frac{d\psi(I^*)}{dI} f(S^*) & -(1-\theta)e^{-(\lambda+\mu)\tau} \\
-\beta I^* & \lambda - \beta S^* + \mu & 0 \\
-\psi(I^*) \frac{df(S^*)}{dS} & -\frac{d\psi(I^*)}{dI} f(S^*) & 1 - \theta e^{-(\lambda+\mu)\tau}
\end{pmatrix}.$$

Then, the characteristic equation is

$$\det(A) = 0,$$

where $\det(A)$ is the determinant of the matrix $A$. We treat the two steady-states separately.

## 4.1 Disease-free equilibrium: $\mathcal{E}_0 := (\Lambda/\mu, 0, 0)$

For this case, we have to write and study the characteristic equation for $\mathcal{E}_0 := (\Lambda/\mu, 0, 0)$.

**Theorem 2.** Assume that $\mathcal{R}_0 < 1$. Then, $\mathcal{E}_0 := (\Lambda/\mu, 0, 0)$ is the only steady-state and is locally asymptotically stable. If $\mathcal{R}_0 > 1$, then the disease-free steady-state $\mathcal{E}_0 := (\Lambda/\mu, 0, 0)$ is unstable.

*Proof.* The characteristic equation is given, for $\mathcal{E}_0$, by

$$\left| \begin{array}{l}
\lambda + \mu \beta S^* + \frac{d\psi(0)}{dI} f\left(\frac{\Lambda}{\mu}\right) - (1-\theta)e^{-(\lambda+\mu)\tau} \\
\qquad 0 \qquad \lambda - \frac{\beta\Lambda-\mu^2}{\mu} \qquad 0 \\
\qquad 0 \qquad \frac{d\psi(0)}{dI} f\left(\frac{\Lambda}{\mu}\right) \qquad -1 + \theta e^{-(\lambda+\mu)\tau}
\end{array} \right| = 0.$$
This is equivalent to

$$ (\lambda + \mu) \left( \lambda - \frac{\beta\Lambda - \mu^2}{\mu} \right) (-1 + \theta e^{-(\lambda+\mu)\tau}) = 0. $$

Then, we have three eigenvalues

$$ -\mu < 0, \quad \mu(\mathcal{R}_0 - 1), \quad \frac{1}{\tau}\ln(\theta) - \mu < 0. $$

We easily conclude on the local asymptotic stability of the disease-free equilibrium.

• If $\mathcal{R}_0 < 1$, then the disease-free equilibrium $\mathcal{E}_0$ is locally asymptotically stable.

• If $\mathcal{R}_0 > 1$, then the disease-free equilibrium $\mathcal{E}_0$ is unstable.

□

## 4.2 Global asymptotic stability of $\mathcal{E}_0 := (\Lambda/\mu, 0, 0)$

Global attractivity of $\mathcal{E}_0$ may be seen with the following arguments. Note

$$ \left( \liminf_{t \to +\infty} S(t), \limsup_{t \to +\infty} S(t) \right) = (S_{\infty}, S^{\infty}), \quad S_{\infty} \le S^{\infty} < +\infty, $$

and

$$ \limsup_{t \to +\infty} I(t) = I^{\infty} < +\infty, \quad \limsup_{t \to +\infty} p(t, 0) = p^{\infty} < +\infty. $$

By (3), we get $S^{\infty} < \Lambda/\mu + \epsilon$, where $\epsilon > 0$ small enough. Then, for a large time, the function $I(t)$ obeys the following inequality

$$ \frac{dI(t)}{dt} \le \mu I(t) \left( \mathcal{R}_0 - 1 + \frac{\epsilon\beta}{\mu} \right). $$

Therefore, for $\mathcal{R}_0 < 1$, $I(t)$ is decreasing for $t$ large enough and tends to zero as $t$ goes to $+\infty$, keeping in mind that $\epsilon$ is very small. Moreover, from the age structured equation of $p$ (third equation of System (1)), when time is large, we have $p(t, \tau) = e^{-\mu\tau}p(t-\tau, 0)$. Then, using the fact that $\psi(0) = 0$, we obtain

$$ p^{\infty} = \frac{\psi(I^{\infty})f(S^{\infty})}{1 - \theta e^{-\mu\tau}} = 0. $$

We directly use Lemma 3.5 of [33] to prove that $p(t, 0)$ tends to zero as $t$ goes to $+\infty$. Using the fluctuation lemma (see Lemma A.14 of [34]), there exist two sequences $(t_n)_{n \in \mathbb{N}}$ and $(s_n)_{n \in \mathbb{N}}$ tending to infinity as $n$ goes to infinity such that $S(t_n) \to S^{\infty}$, $S(s_n) \to S_{\infty}$, $\frac{dS(t_n)}{dt} \to 0$ and $\frac{dS(s_n)}{dt} \to 0$. We conclude from the equation of $S$ (first equation of System (1)) that $S^{\infty} = S_{\infty} = \Lambda/\mu$. Consequently, we have the following theorem.
**Theorem 3.** Assume that $\mathcal{R}_0 < 1$. Then, the disease-free steady-state $\mathcal{E}_0 := (\Lambda/\mu, 0, 0)$ is globally asymptotically stable.

## 4.3 Persistence analysis

In this section, we examine the persistence of the solution of the system (1) (i.e., when the solution remains positive) in the presence of an unstable disease-free steady-state. For technical reasons, we consider the reduced system (2) instead of the original system (1).

We prove the following lemma on the weak persistence of the disease.

**Lemma 4.** Assume that $\mathcal{R}_0 > 1$. Then, there exists $\epsilon > 0$ such that

$$ \limsup_{t \to +\infty} I(t) > \epsilon, $$

for any initial condition $(S_0, I_0, \phi) \in \mathbb{R}^+ \times \mathbb{R}^{+*} \times C([-\tau, 0], \mathbb{R}^+)$ of the system (2).

*Proof.* Let us prove by contradiction that $\limsup_{t \to +\infty} I(t) > \epsilon$.

We assume that $\limsup_{t \to +\infty} I(t) \le \epsilon$. Then, there exists a sufficiently large $t_\epsilon > 0$ such that $I(t) \le \epsilon$, for all $t \ge t_\epsilon$. Therefore, we have for all $t \ge t_\epsilon$

$$ \frac{dS(t)}{dt} \ge \Lambda - \mu S(t) - \beta \epsilon S(t) - \psi(\epsilon) f(S(t)). $$

If we set $\liminf_{t \to +\infty} S(t) = S_\infty$, then there exists a sequence $(t_k)_{k \in \mathbb{N}}$ tending to infinity as $k$ tends to infinity such that $S(t_k) \to S_\infty$ and $S'(t_k) \to 0$. Necessarily, we have

$$ S_{\infty} \ge F_{\epsilon}^{-1}(\Lambda), \quad \text{with} \quad F_{\epsilon}(x) = \mu x + \beta \epsilon x + \psi(\epsilon) f(x), \quad x \ge 0. $$

Remember that for $\epsilon = 0$, we have $F_0^{-1}(\Lambda) = S^0 = \Lambda/\mu$ and $\mathcal{R}_0 = \beta S^0/\mu$. Hence, by using $\mathcal{R}_0 > 1$, we can choose $\epsilon > 0$ and $\nu > 0$ sufficiently smalls such that

$$ \mathcal{R}_0^\epsilon := \frac{\beta S_\epsilon^0}{\nu + \mu} > 1, \quad \text{where} \quad S_\epsilon^0 := F_\epsilon^{-1}(\Lambda). $$

Multiplying the equation of $I(t)$ by $e^{-\nu t}$ and integrating from $t_\epsilon$ to $+\infty$, we obtain

$$ 0 > -e^{-\nu t_\epsilon} I(t_\epsilon) > \left[-(\nu + \mu) + \beta S_\epsilon^0\right] \int_{t_\epsilon}^{+\infty} e^{-\nu t} I(t) dt > 0. $$

This is a contradiction with the positivity of the solution. Hence, the conclusion of the lemma holds true. □

We now give the following result on the strong persistence of the disease, see [34]. We use the same idea as in the paper [35, the proof of Theorem 1].
**Theorem 5.** Suppose that $\mathbb{R}_0 > 1$. There exists a positive constant $\epsilon_1 > 0$ such that $\liminf_{t \to +\infty} I(t) > \epsilon_1$, for any initial condition $(S_0, I_0, \phi) \in \mathbb{R}^+ \times \mathbb{R}^{+*} \times C([-\tau, 0], \mathbb{R}^+)$ of the system (2).

*Proof.* By using the weak persistence proved in Lemma 4, we have

$$ \limsup_{t \to +\infty} I(t) > \epsilon $$

for a given $\epsilon > 0$. This yields to the existence of an increasing sequence $(\eta_k)_{k \in \mathbb{N}}$, which is positive, with $\eta_k \to +\infty$ such that $I(\eta_k) > \epsilon$. Suppose by contradiction that for all $\tilde{\epsilon} \in (0, \epsilon]$ there exists an initial condition $(S_0, I_0, \phi) \in \mathbb{R}^+ \times \mathbb{R}^{+*} \times C([-\tau, 0], \mathbb{R}^+)$, such that

$$ \liminf_{t \to +\infty} I(t) \le \tilde{\epsilon}. $$

There exists a positive increasing sequence $(t_k)_{k \in \mathbb{N}}$ and a positive decreasing sequence $(t'_k)_{k \in \mathbb{N}}$ with $t_k > \eta_k$, $\lim_{k \to +\infty} t'_k = 0$ and

$$ I(t_k) < t'_k < \tilde{\epsilon} \le \epsilon. \quad (5) $$

Therefore, there exists a sequence $(\nu_k)_{k \in \mathbb{N}}$ where $\nu_k \in (\eta_k, t_k)$ such that

$$ I(\nu_k) = \epsilon \quad \text{and} \quad I(t) < \epsilon, \quad \text{for all } t \in (\nu_k, t_k). \quad (6) $$

Let $(I_k)_{k \in \mathbb{N}}$ and $(S_k)_{k \in \mathbb{N}}$ be the sequences such that $I_k := I(\nu_k) = \epsilon$ and $S_k := S(\nu_k) \in \mathbb{R}^+$. We point out that the sequence $(I_k)_{k \in \mathbb{N}}$ is constant and the sequence $(S_k)_{k \in \mathbb{N}}$ is uniformly bounded. Accordingly to the boundedness, there exist convergent subsequences of $(I_k)_{k \in \mathbb{N}}$ and $(S_k)_{k \in \mathbb{N}}$ (we keep the same notation for these two subsequences). Then, we have $I_k = \epsilon$ and $\lim_{k \to +\infty} S_k = l$, with $l \in \mathbb{R}^+$. Consider the following problem

$$ v(t) = \begin{cases} \psi(\epsilon)f(l) + \theta e^{-\mu\tau}v(t - \tau), & t > 0, \\ \phi(t), & t \in [-\tau, 0]. \end{cases} $$

It admits a unique continuous solution $v$ on $(0, +\infty)$, for any initial condition $\phi \in C([-\tau, 0], \mathbb{R}^+)$. Let $(u_k)_{k \in \mathbb{N}}$ be the functional sequence in $C([-\tau, 0], \mathbb{R}^+)$ defined by $u_k(\theta) := v(\nu_k + \theta)$, $\theta \in [-\tau, 0]$ with $\nu_k > \tau$ for $k$ large enough. Without loss of generality, we can consider $\nu_k > \tau$, for all $k \in \mathbb{N}$. Thus,

$$ u_k(\theta) = \psi(\epsilon)f(l) + \theta e^{-\mu\tau}u_k(\theta - \tau), \quad \text{for all } \theta \in [-\tau, 0]. $$

Clearly, we have the uniform boundedness of the sequence $(u_k)_{k \in \mathbb{N}}$. By taking $\alpha, \alpha' \in [-\tau, 0]$, we have

$$ \begin{aligned} |u_k(\alpha) - u_k(\alpha')| &= \theta e^{-\mu\tau} |u_k(\alpha - \tau) - u_k(\alpha' - \tau)|, \\ &\le (\theta e^{-\mu\tau})^{n_k+1} |\phi(\nu_k + \alpha - (n_k + 1)\tau) - \phi(\nu_k + \alpha' - (n_k + 1)\tau)|, \end{aligned} $$
$$ \leq |\phi(\nu_k + \alpha - (n_k + 1)\tau) - \phi(\nu_k + \alpha' - (n_k + 1)\tau)|, $$

where $n_k := \lfloor \nu_k / \tau \rfloor$. By the fact that $\phi$ is uniformly continuous on $[-\tau, 0]$, then the last inequality implies that the sequence $(u_k)_{k \in \mathbb{N}}$ is equicontinuous. It follows from the Ascoli-Arzela theorem that there exists a convergent subsequence of $(u_k)_{k \in \mathbb{N}}$ (we denote it again $(u_k)_{k \in \mathbb{N}}$). Then, $\lim_{k \to +\infty} u_k = u^*$ and $u^* \in C([-\tau, 0], \mathbb{R}^+)$. We consider now the solution, denoted $(S^\infty, I^\infty, u^\infty)$, of System (2), with the initial condition $S_0 = l$, $I_0 = \epsilon$ and $\phi = u^* \in C([-\tau, 0], \mathbb{R}^+)$. From Lemma 4, there exists $\varrho > 0$ and we can find $0 < m < \epsilon$, such that

$$ I^{\infty}(\varrho) > \epsilon \quad \text{and} \quad I^{\infty}(t) > m, \quad \text{for all } t \in (0, \varrho). \tag{7} $$

Now, we are going to prove that a contradiction occurs. For each $k \in \mathbb{N}$, we put $\tilde{I}_k(t) := I(\nu_k + t)$, $t > 0$. From (7), the continuity and the fact that

$$ \tilde{I}_k(0) = I_k = \epsilon, \quad \lim_{k \to +\infty} S_k = l \quad \text{and} \quad \lim_{k \to +\infty} u_k = u^*, $$

we have, for $k$ large enough and by $\lim_{k \to +\infty} t'_k = 0$,

$$ \tilde{I}_k(\varrho) > \epsilon \quad \text{and} \quad \tilde{I}_k(t) > m > t'_k, \quad \text{for all } t \in (0, \varrho). \tag{8} $$

In a parallel way, by putting $\tilde{t}_k := t_k - \nu_k$, we have from (5) and (6) that

$$ \tilde{I}_k(\tilde{t}_k) = I(t_k) < t'_k < \epsilon \quad \text{and} \quad \tilde{I}_k(t) = I(\nu_k + t) < \epsilon, \quad \text{for all } t \in (0, \tilde{t}_k). \tag{9} $$

We distinguish the following cases:

* If $\varrho < \tilde{t}_k$, then the second inequality of (9) which gives $\tilde{I}_k(\varrho) < \epsilon$ and the first inequality of (8) are in contradiction.

* If $\varrho = \tilde{t}_k$, then the first inequality of (8) and the first inequality of (9) are in contradiction.

* If $\tilde{t}_k < \varrho$, then the second inequality of (8) which gives $\tilde{I}_k(\tilde{t}_k) > t'_k$ and the first inequality of (9) are in contradiction.

Consequently, there exists $\tilde{\epsilon} \in (0, \epsilon]$ such that for any initial condition $(S_0, I_0, \phi) \in \mathbb{R}^+ \times \mathbb{R}^{+*} \times C([-\tau, 0], \mathbb{R}^+)$, we have

$$ \liminf_{t \to +\infty} I(t) > \tilde{\epsilon}. $$

□

This result suggests that given a set of parameters where the basic reproduction number $R_0 > 1$, the number of infected will not die out by itself, and the disease will become endemic. The number of infected will always be greater than a strictly positive constant $\epsilon$. Thus, any outbreak related to the disease will last in time because there is always infection. PrEP is then effective by reducing the reproduction number and has an impact on the persistence of the disease.
## 4.4 The dynamical behavior of the endemic steady state

Let us focus now on the stability of the endemic steady-state. Suppose that $\tau > 0$. We have

$$S^* = \frac{\mu}{\beta},$$

$$\mathcal{H}(I^*) := \frac{\beta\Lambda}{\mu^2 + \beta\mu I^* + \left[ \frac{\beta(1 - e^{-\mu\tau})}{1 - \theta e^{-\mu\tau}} \right] \psi(I^*) f\left(\frac{\mu}{\beta}\right)} = 1,$$

$$p^*(a) = \frac{e^{-\mu a}}{1 - \theta e^{-\mu \tau}} \psi(I^*) f\left(\frac{\mu}{\beta}\right).$$

Then, the characteristic equation is given by $\det(A) = 0$, where the determinant of the matrix $A$ is

$$\det(A) = \begin{vmatrix} \lambda + \mu + \beta I^* + \psi(I^*) \frac{df(S^*)}{dS} & \beta S^* + \frac{d\psi(I^*)}{dI} f(S^*) & -(1-\theta)e^{-(\lambda+\mu)\tau} \\ -\beta I^* & \lambda - \beta S^* + \mu & 0 \\ -\psi(I^*) \frac{df(S^*)}{dS} & -\frac{d\psi(I^*)}{dI} f(S^*) & 1 - \theta e^{-(\lambda+\mu)\tau} \end{vmatrix}.$$

To simplify the expression of the determinant, we add the three lines to the first line and we obtain

$$\begin{vmatrix} \lambda + \mu & \lambda + \mu & 1 - e^{-(\lambda+\mu)\tau} \\ -\beta I^* & \lambda - \beta S^* + \mu & 0 \\ -\psi(I^*) \frac{df(S^*)}{dS} & -\frac{d\psi(I^*)}{dI} f(S^*) & 1 - \theta e^{-(\lambda+\mu)\tau} \end{vmatrix} = 0.$$

This implies the following factorization

$$(\lambda + \mu) \begin{vmatrix} 1 & 1 & \frac{1 - e^{-(\lambda + \mu)\tau}}{\lambda + \mu} \\ -\beta I^* & \lambda - \beta S^* + \mu & 0 \\ -\psi(I^*) \frac{df(S^*)}{dS} & -\frac{d\psi(I^*)}{dI} f(S^*) & 1 - \theta e^{-(\lambda + \mu)\tau} \end{vmatrix} = 0.$$

Remark that the following limit is finite

$$\lim_{\lambda \to -\mu} \frac{1 - e^{-(\lambda + \mu)\tau}}{\lambda + \mu} = \tau \neq 0.$$
Therefore, $\lambda = -\mu < 0$ is a root of the characteristic equation. Other roots verify

$$ \left| \begin{array}{ccc} 1 & 0 & \frac{1 - e^{-(\lambda+\mu)\tau}}{\lambda+\mu} \\ -\beta I^* & \lambda - \beta(S^* - I^*) + \mu & 0 \\ -\psi(I^*) \frac{df(S^*)}{dS} & \psi(I^*) \frac{df(S^*)}{dS} - \frac{d\psi(I^*)}{dI} f(S^*) & 1 - \theta e^{-(\lambda+\mu)\tau} \end{array} \right| = 0. $$

By observing that $\beta S^* - \mu = 0$, we develop the determinant following the first line

$$ (\lambda + \beta I^*) (1 - \theta e^{-(\lambda+\mu)\tau}) + \frac{1 - e^{-(\lambda+\mu)\tau}}{\lambda + \mu} \left| \begin{array}{ccc} -\beta I^* & \lambda + \beta I^* & 0 \\ -\psi(I^*) \frac{df(S^*)}{dS} & \psi(I^*) \frac{df(S^*)}{dS} & -\frac{d\psi(I^*)}{dI} f(S^*) \end{array} \right| = 0. $$

This implies that

$$ (\lambda + \beta I^*) (1 - \theta e^{-(\lambda+\mu)\tau}) + \frac{1 - e^{-(\lambda+\mu)\tau}}{\lambda + \mu} \left( -\beta I^* (\psi(I^*) \frac{df(S^*)}{dS} - \frac{d\psi(I^*)}{dI} f(S^*)) + \frac{1 - e^{-(\lambda+\mu)\tau}}{\lambda + \mu} \left( \psi(I^*) \frac{df(S^*)}{dS} (\lambda + \beta I^*) \right) \right) = 0. $$

To simplify the notations, we note

$$ b := \beta I^* \frac{d\psi(I^*)}{dI} f(S^*) > 0, $$

$$ c := \psi(I^*) \frac{df(S^*)}{dS} > 0. $$

With these notations, we define the function $\mathcal{F}$ by

$$ \mathcal{F}(\lambda) := (\lambda + \mu)(\lambda + \beta I^*) (1 - \theta e^{-(\lambda+\mu)\tau}) + (c\lambda + b) (1 - e^{-(\lambda+\mu)\tau}) = 0. \quad (10) $$

We easily prove that $\mathcal{F}$ is strictly increasing on $[0, +\infty)$ due to the positive parameters. Moreover, we have

$$ \lim_{\lambda \to +\infty} \mathcal{F}(\lambda) = +\infty. $$

Furthermore,

$$ \mathcal{F}(0) = \mu\beta I^* (1 - \theta e^{-\mu\tau}) + (1 - e^{-\mu\tau})b > 0. $$
Then, the function $\mathcal{F}$ has no positive real root. We also remark

$$ \mathcal{F}(\lambda) := \left(1 - \theta e^{-(\lambda+\mu)\tau}\right) (c\lambda + b) \left[ \frac{(\lambda+\mu)(\lambda+\beta I^*)}{c\lambda + b} + \frac{1 - e^{-(\lambda+\mu)\tau}}{1 - \theta e^{-(\lambda+\mu)\tau}} \right], $$

with

$$ \lim_{\lambda \to -\infty} \frac{1 - e^{-(\lambda+\mu)\tau}}{1 - \theta e^{-(\lambda+\mu)\tau}} = \frac{1}{\theta}. $$

Then,

$$ \lim_{\lambda \to -\infty} \mathcal{F}(\lambda) = -\infty. $$

This means there is a negative real root and all real roots are negative.

- Case $\tau = 0$: the characteristic equation is simplified to

$$ \mathcal{F}(\lambda) := (\lambda + \mu)(\lambda + \beta I^*) (1 - \theta) = 0. $$

It has two roots $\lambda = -\mu < 0$ and $\lambda = -\beta I^* < 0$ which are negative. This implies that the endemic steady-state is locally asymptotically stable.

- Case $\tau > 0$: we need to go further in our analysis and investigate the behavior of complex roots. Before that, we check that there is no infinite sequence of complex roots $\{\lambda_n\}_n$ with strict positive real parts. We suppose by contradiction that there exists a sequence of roots $\{\lambda_n\}_n$ with

$$ \lim_{n \to +\infty} |\lambda_n| = +\infty. $$

In this case, every $\lambda_n$ satisfies

$$ \left(1 + \frac{\beta I^*}{\lambda_n}\right) \left(1 - \theta e^{-(\lambda_n + \mu)\tau}\right) + \frac{1 - e^{-(\lambda_n + \mu)\tau}}{\lambda_n + \mu} \left(c + \frac{b}{\lambda_n}\right) = 0. $$

By taking the limit for $n \to \infty$, we obtain

$$ \lim_{n \to \infty} \left(1 - \theta e^{-(\lambda_n + \mu)\tau}\right) = 0. $$

This implies that

$$ \lim_{n \to \infty} e^{\Re(\lambda_n)\tau} = \theta e^{-\mu\tau}. $$

Therefore

$$ \lim_{n \to \infty} \Re(\lambda_n) = \frac{1}{\tau} \ln(\theta e^{-\mu\tau}) = \frac{1}{\tau} \ln(\theta) - \mu < 0. $$

Then, the roots have negative real parts. At this step, we conclude with the following proposition.

**Proposition 6.** There exists some $\tau_0^* > 0$ such that for $\tau \in [0, \tau_0^*)$ all characteristic roots have negative real parts.
Next, we prove the existence or non-existence of complex roots with positive real parts.

- Case $\tau = 0$: we first notice again the non-existence of such roots.

- Case $\tau > 0$: we prove the non-existence or the existence of pure complex roots of the form $iw$, with $w > 0$. We rewrite the equation (10) as, for $\lambda \in \mathbb{C}$,

$$ \mathcal{F}(\lambda) := \lambda^2 + (\beta I^* + \mu + c)\lambda + \mu\beta I^* + b - e^{-(\lambda+\mu)\tau} [\theta\lambda^2 + (\theta\beta I^* + \theta\mu + c)\lambda + \theta\mu\beta I^* + b] = 0. \quad (11) $$

For $\lambda = iw$ and by separating the real and the imaginary part in (11), we get

$$ \begin{cases} A \cos(w\tau) - Bw \sin(w\tau) = w^2 - \mu\beta I^* - b, \\ Bw \cos(w\tau) + A \sin(w\tau) = w(\beta I^* + \mu + c), \end{cases} $$

where

$$ A = e^{-\mu\tau}[\theta w^2 - \theta\mu\beta I^* - b] \quad \text{and} \quad B = e^{-\mu\tau}[\theta\beta I^* + \theta\mu + c]. $$

It follows that

$$ \left\{ \begin{aligned} \cos(w\tau) &= \frac{A(w^2 - \mu\beta I^* - b) + w^2 B(\beta I^* + \mu + c)}{A^2 + (Bw)^2}, \\ \sin(w\tau) &= \frac{wA(\beta I^* + \mu + c) - Bw(w^2 - \mu\beta I^* - b)}{A^2 + (Bw)^2}. \end{aligned} \right. \quad (12) $$

Then, we get

$$ A^2 + w^2 B^2 = (w^2 - \mu\beta I^* - b)^2 + w^2(\beta I^* + \mu + c)^2. $$

Hence, it implies that $w$ must be a root of the following equation

$$ (1 - (\theta e^{-\mu\tau})^2)w^4 + (-2(\mu\beta I^* + b) + 2\theta e^{-2\mu\tau}(\theta\mu\beta I^* + b) + \vartheta(\tau))w^2 + (\mu\beta I^* + b)^2 - e^{-2\mu\tau}(\theta\mu\beta I^* + b)^2 = 0, $$

where $\vartheta(\tau) := (\beta I^* + \mu + c)^2 - e^{-2\mu\tau}(\theta\beta I^* + \mu\theta + c)^2 \ge 0$. We note $X = w^2$ and analyze a second-degree polynomial function for $X \ge 0$. If $(-2(\mu\beta I^* + b) + 2\theta e^{-2\mu\tau}(\theta\mu\beta I^* + b) + \vartheta(\tau))$ is nonnegative for $\tau \in [0, \tau_1^*)$ with $\tau_1^* > 0$, then the Routh-Hurwitz criterion implies that no positive real roots exist. So, we have the following property.

**Proposition 7.** Assume that $\mathcal{R}_0 > 1$ and $\vartheta(\tau) \ge \vartheta^*(\tau) := 2(\mu\beta I^* + b) - 2\theta e^{-2\mu\tau}(\theta\mu\beta I^* + b)$ for all $\tau \in [0, \tau_1^*)$ with $\tau_1^* > 0$. Then, the endemic steady-state $\mathcal{E}^* := (S^*, I^*, p^*(a))$ is locally asymptotically stable for $\tau \in [0, \tau_1^*)$.

To illustrate the previous proposition, we give the following two examples associated with the bilinear case and a saturated case with a function as in [29] which is more realistic.
**Example 1.** Let $\psi$ and $f$ given by

$$ \psi(I) = l_1 I, \quad f(S) = l_2 S, \quad \text{with } l_1, l_2 > 0. $$

This set-up would describe a population of small size or living in a rich country where no competition effect would occur when individuals are seeking treatment. The choice of the function $\psi$ means that recruitment is motivated by the number of infected with no restriction from political or economic decisions.

We note $l := l_1 l_2$. The condition $\vartheta(\tau) \ge 2(\mu\beta I^* + b) - 2\theta e^{-2\mu\tau}(\theta\mu\beta I^* + b)$ is always true for any $\tau \ge 0$. Then, the endemic equilibrium is always asymptotically stable (see Figure 2). In fact, in this case, we have

$$ c := \psi(I^*) \frac{df(S^*)}{dS} = lI^* \quad \text{and} \quad b := \beta I^* \frac{d\psi(I^*)}{dI} f(S^*) = \mu l I^* = \mu c. $$

Therefore, we get

$$ \begin{aligned} \vartheta(\tau) - \vartheta^*(\tau) &= (\beta I^* + \mu + c)^2 - e^{-2\mu\tau}(\theta\beta I^* + \mu\theta + c)^2 \\ & \quad - 2(\mu\beta I^* + \mu c) - 2\theta e^{-2\mu\tau}(\theta\mu\beta I^* + \mu c), \\ &= [(\beta + l)^2 - e^{-2\mu\tau}(\theta\beta + l)^2](I^*)^2 + \mu^2(1 - \theta^2 e^{-2\mu\tau}), \\ & \ge 0. \end{aligned} $$

Thus, the endemic steady-state $\mathcal{E}^* := (S^*, I^*, p^*(a))$ is locally asymptotically stable for $\tau \ge 0$.

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Left: we observe that the condition $\vartheta(\tau) \ge \vartheta^*(\tau)$ of Proposition 7 is satisfied for the selected parameters. Right: a simulation is carried out for the stability of the endemic equilibrium. Parameters are $\Lambda = 50$, $\mu = 0.1$, $\beta = 0.001$, $\theta = 0.1$, $l_1 = 0.01$, $l_2 = 0.3$ and $\tau = 3.5$. In this case $\mathcal{R}_0 = 5 > 1$.</figcaption></figure>
**Example 2.** Let $\psi$ and $f$ given by

$$ \psi(I) = \frac{a_I I}{b_I + I}, \quad f(S) = a_S S, \quad \text{with } a_I, a_S, b_I > 0. $$

The choice of function would describe a saturation effect in recruitment due to a high number of infected. It would describe an outbreak in a very specific population where the ones targeted by PrEP and the very persons supposed to begin treatment are unable as they are infected.

We check that the condition $\vartheta(\tau) \ge 2(\mu\beta I^* + b) - 2\theta e^{-2\mu\tau}(\theta\mu\beta I^* + b)$ is always true for any $\tau \ge 0$ (see Figure 3). In this case, we have

$$ c := \psi(I^*) \frac{df(S^*)}{dS} = \frac{a_S a_I}{b_I + I^*} I^* =: l(I^*) I^* $$

and

$$ b := \mu L(I^*) c, \quad \text{where } L(I^*) := \frac{b_I}{b_I + I^*} < 1. $$

After an easy computation, we get

$$ \begin{aligned} \vartheta(\tau) - \vartheta^*(\tau) &= (\beta I^* + \mu + c)^2 - e^{-2\mu\tau}(\theta\beta I^* + \mu\theta + c)^2 \\ &\quad - 2(\mu\beta I^* + \mu L(I^*)c) - 2\theta e^{-2\mu\tau}(\theta\mu\beta I^* + \mu L(I^*)c), \\ &= [(\beta + l(I^*))^2 - e^{-2\mu\tau}(\theta\beta + l(I^*))^2](I^*)^2 + \mu^2(1 - \theta^2 e^{-2\mu\tau}) \\ &\quad + 2\mu l(I^*) I^* (1 - L(I^*)) (1 - \theta e^{-2\mu\tau}), \\ &\ge 0. \end{aligned} $$

Consequently, the endemic steady-state $\mathcal{E}^* := (S^*, I^*, p^*(a))$ is locally asymptotically stable for $\tau \ge 0$.

Next, we suppose that $\vartheta(\tau) < \vartheta^*(\tau) = 2(\mu\beta I^* + b) - 2\theta e^{-2\mu\tau}(\theta\mu\beta I^* + b)$ for some $\tau > 0$. After tedious calculations, the discriminant is given by

$$ \Delta(\vartheta(\tau)) = \vartheta^2(\tau) - 4[(\mu\beta I^* + b) - \theta e^{-2\mu\tau}(\theta\mu\beta I^* + b)]\vartheta(\tau) \\ + [2e^{-\mu\tau}\theta(\mu\beta I^* + b) - 2e^{-\mu\tau}(\theta\mu\beta I^* + b)]^2. $$

Furthermore, the discriminant of the last polynomial of the second degree is

$$ \begin{aligned} \bar{\Delta} &= 16[(1 + \theta e^{-\mu\tau})(\mu\beta I^* + b) - (\theta e^{-2\mu\tau} + e^{-\mu\tau})(\theta\mu\beta I^* + b)] \\ &\quad \times [(1 - \theta e^{-\mu\tau})(\mu\beta I^* + b) - (\theta e^{-2\mu\tau} - e^{-\mu\tau})(\theta\mu\beta I^* + b)], \\ &\ge 0. \end{aligned} $$

Define

$$ \vartheta_-(\tau) := 2[(\mu\beta I^* + b) - \theta e^{-2\mu\tau}(\theta\mu\beta I^* + b)] - \frac{\sqrt{\bar{\Delta}}}{2} \ge 0. $$
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Left: we observe that the condition of Proposition 7 is satisfied for the selected parameters. Right: a simulation is carried out for the stability of the endemic equilibrium. Parameters are $\Lambda = 50$, $\mu = 0.15$, $\beta = 0.001$, $\theta = 0.1$, $a_S = 0.3$, $a_I = 0.5$, $b_I = 1$ and $\tau = 3.5$. In this case $\mathcal{R}_0 = 2.22 > 1$.</figcaption></figure>

Clearly, the discriminant $\Delta(\vartheta(\tau))$ is either negative for $\vartheta_-(\tau) < \vartheta(\tau) < \vartheta^*(\tau)$ or nonnegative for $0 \le \vartheta(\tau) \le \vartheta_-(\tau)$.

**Proposition 8.** Assume that $\mathcal{R}_0 > 1$ and $\vartheta_-(\tau) < \vartheta(\tau) < \vartheta^*(\tau)$ for all $\tau \in [0, \tau_2^*)$ with $\tau_2^* > 0$. Then, the endemic steady-state $\mathcal{E}^* := (S^*, I^*, p^*(a))$ is locally asymptotically stable for $\tau \in [0, \tau_2^*)$.

Now, we focus on the case where $0 \le \vartheta(\tau) \le \vartheta_-(\tau)$ for some $\tau > 0$ (see example given in Figure 4). Then, there are two possible roots $w_{\pm}(\tau)$ given by

$$w_{-}(\tau) = \left[ \frac{-(-2(\mu\beta I^{*} + b) + 2\theta e^{-2\mu\tau}(\theta\mu\beta I^{*} + b) + \vartheta(\tau)) - \sqrt{\Delta(\vartheta(\tau))}}{2(1 - (\theta e^{-\mu\tau})^2)} \right]^{\frac{1}{2}},$$

and

$$w_{+}(\tau) = \left[ \frac{-(-2(\mu\beta I^{*} + b) + 2\theta e^{-2\mu\tau}(\theta\mu\beta I^{*} + b) + \vartheta(\tau)) + \sqrt{\Delta(\vartheta(\tau))}}{2(1 - (\theta e^{-\mu\tau})^2)} \right]^{\frac{1}{2}}.$$

From (12), we define the angles $\Theta_{\pm}(\tau) \in [0, 2\pi]$ solution of

$$\begin{cases} \cos(\Theta_{\pm}(\tau)) = \frac{A(w_{\pm}^2(\tau) - \mu\beta I^* - b) + w_{\pm}^2(\tau)B(\beta I^* + \mu + c)}{A^2 + (Bw_{\pm}(\tau))^2}, \\ \sin(\Theta_{\pm}(\tau)) = \frac{w_{\pm}(\tau)A(\beta I^* + \mu + c) - Bw_{\pm}(\tau)(w_{\pm}^2(\tau) - \mu\beta I^* - b)}{A^2 + (Bw_{\pm}(\tau))^2}. \end{cases}$$
We get $w_{\pm}(\tau)\tau = \Theta_{\pm}(\tau) + 2k\pi$ with $k \in \mathbb{N}$. A deep analysis of the existence of such angles can be referred to [36]. We also get

$$ \Theta_{\pm}(\tau) = \arctan(H_1(\tau)/H_2(\tau)), \quad \text{if } \sin(\Theta_{\pm}(\tau)) > 0, \cos(\Theta_{\pm}(\tau)) > 0, $$

$$ \Theta_{\pm}(\tau) = \pi/2, \quad \text{if } \sin(\Theta_{\pm}(\tau)) = 1, \cos(\Theta_{\pm}(\tau)) = 0, $$

$$ \Theta_{\pm}(\tau) = \pi + \arctan(H_1(\tau)/H_2(\tau)), \quad \text{if } \cos(\Theta_{\pm}(\tau)) < 0, $$

$$ \Theta_{\pm}(\tau) = 3\pi/2, \quad \text{if } \sin(\Theta_{\pm}(\tau)) = -1, \cos(\Theta_{\pm}(\tau)) = 0, $$

$$ \Theta_{\pm}(\tau) = 2\pi + \arctan(H_1(\tau)/H_2(\tau)), \quad \text{if } \sin(\Theta_{\pm}(\tau)) < 0, \cos(\Theta_{\pm}(\tau)) > 0, $$

with

$$ H_1(\tau) := w_{\pm}(\tau)A(\beta I^* + \mu + c) - Bw_{\pm}(\tau)(w_{\pm}^2(\tau) - \mu\beta I^* - b), $$

$$ H_2(\tau) := A(w_{\pm}^2(\tau) - \mu\beta I^* - b) + w_{\pm}^2(\tau)B(\beta I^* + \mu + c). $$

We keep in mind that the endemic equilibrium exists for $\mathbb{R}_0 > 1$ and $\tau$ ranges on $\mathbb{R}^+$. To obtain a critical value of $\tau$ for the stability change, we define the following functions, for $k \in \mathbb{N}$,

$$ T_k^{\pm} : \mathbb{R}^{+} \to \mathbb{R}, \\ \tau \mapsto T_k^{\pm}(\tau) = \tau - \frac{\Theta_{\pm}(\tau) + 2k\pi}{w_{\pm}(\tau)}. $$

Note that a pair of simple conjugate pure imaginary roots $\lambda = \pm iw(\tau^*)$ exist at $\tau^*$ if and only if $\tau^*$ is a root of $T_k^+$ or $T_k^-$, $k \in \mathbb{N}$. Then, the stability may switch at any root of $T_k^\pm$. We are now ready to conclude with the following theorem.

**Theorem 9.** Let $\mathbb{R}_0 > 1$. If $T_k^+$ has no root on $\mathbb{R}^+$, for any $k \in \mathbb{N}$, then the endemic equilibrium is locally asymptotically stable for all $\tau \in \mathbb{R}^+$. If $T_k^+$ or $T_k^-$ has at least one root in $\mathbb{R}^+$ and if we denote by $\tau^*$ the smallest root among all, then the endemic equilibrium is locally asymptotically stable for $\tau \in [0, \tau^*)$.

**Example 3.** Let $\psi$ and $f$ given by

$$ \psi(I) = a_I I, \quad f(S) = \frac{a_S S}{b_S + S}, \quad \text{with } a_I, a_S, b_S > 0. $$

For this last example, the situation described would take place in a poor country or in a very large population where we would observe competition between individuals. For a given parameter, we check that the condition $0 \le \vartheta(\tau) \le \vartheta_-(\tau)$ is true for $\tau \in [0, 6]$ in Figure 4.

In this example, we remark that the functions $T_k^+$ plotted do not admit any roots and this is valid for all the values of the parameters. So, we conclude that the characteristic equation cannot have a purely imaginary root and therefore no destabilization of the endemic equilibrium is possible.
<figure><img src="image_5.png" /><figcaption>**Fig. 4** We observe that the condition $0 \le \vartheta(\tau) \le \vartheta_-(\tau)$ is satisfied for a selected parameters. Parameters are $\Lambda = 50, \mu = 0.1, \beta = 0.001, \theta = 0.1, a_S = 1, a_I = 2, b_S = 0.2$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Functions $T_k^{\pm}$ are drawn for $\tau \in [0, 6]$. We observe that these functions exhibit no root in this interval where a possible purely imaginary root could exist. In this case, no possible Hopf bifurcation occurs. Parameters are $\Lambda = 50, \mu = 0.1, \beta = 0.001, \theta = 0.1, a_S = 1, a_I = 2, b_S = 0.2$.</figcaption></figure>

# 5 Discussion

In this study, we concentrated on an epidemiological SI (Susceptible-Infected) model that characterizes HIV epidemics, with an additional external compartment denoted as $P$ to represent Pre-Exposure Prophylaxis (PrEP) treatment. We assumed the $P$ compartment to be age-structured, where age corresponds to the time elapsed since the last PrEP prescription. Our goal is to introduce a novel framework to assess the effects of PrEP on HIV transmission dynamics. Our approach builds upon prior models [28, 30, 37, 38], in which the initiation rate ($\psi$) of the protection phase remained constant. However, unlike these works, we now consider $\psi$ to be dependent on the infected population $I$, emphasizing the influence of the HIV epidemic on the protection phase.
Employing the method of characteristics, we derived a coupled system of differential-difference equations with delay. We pinpointed two steady states – the disease-free equilibrium and the endemic equilibrium – and scrutinized their stability and existence using the classic criterion associated with the basic reproduction number $R_0$. We confirmed the global asymptotic stability of the disease-free equilibrium and investigated the system's persistence when $R_0 > 1$.

Remarkably, the basic reproduction number of the system reflects that of any standard SI model due to $\psi$ being an $I$-dependent function, with the assumption $\psi(0) = 0$. This differs from a previous study [28], where $\psi$ was assumed to be constant, resulting in a more complex expression for $R_0$ involving delay dependence. Since $\psi$ now varies with $I$, the disease-free equilibrium relies solely on the susceptible component, unlike in [28]. However, the equilibrium and its stability are affected by the delay. We demonstrated the weak persistence of the disease, leading to the proof of strong persistence. Subsequently, we analyzed the stability of the endemic state, uncovering various stability regions. While the case with zero delay is straightforward, the non-trivial case presents challenges, and we could not fully determine the occurrence of Hopf bifurcations.

We maintain that the occurrence of a Hopf bifurcation in this system is unlikely, supported by heuristic reasoning from the theory of monotone systems. However, further investigation is necessary to confirm this. It is important to note that the behavior of the recruitment function was based on pre-COVID-19 data, and the pandemic likely affected data collection and access to treatment services, including PrEP. Nevertheless, it appears that the recruitment rate has returned to pre-pandemic levels. Future research will focus on the impact of various policies, economic factors, and cultural influences from different countries on the evolution of the epidemic, as well as on the development of STIs. We have included numerical simulations to visually demonstrate the established stability conditions.

As far as we know, this study presents a more realistic depiction of PrEP's influence on new HIV cases and introduces fresh theoretical perspectives on stability. In essence, our approach offers a more efficient method of modeling the HIV epidemic, considering political and economic decisions regarding PrEP treatment for high-risk populations.

**Acknowledgments.** M.A. thanks the Inria international associate team DIFAIRE. G.R. was supported by the Natural Sciences and Engineering Research Council and the York Research Chairs program. This work was performed within the framework of the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). A.C. thanks grant, PRFU: C00L03UN29012022002, from DGRSDT of Algeria. This work was partially supported by a grant from the IMU-CDC and Simons Foundation.

## Declarations

The authors declare there is no conflict of interest.
References

[1] Stoddart, C.A., Reyes, R.A.: Models of HIV-1 disease: A review of current status. Drug Discovery Today: Disease Models **3**(1), 113–119 (2006). doi.org/10.1016/j.ddmod.2006.03.016

[2] The path that ends AIDS: UNAIDS Global AIDS Update 2023. Geneva: Joint United Nations Programme on HIV/AIDS; 2023. Licence: CC BY-NC-SA 3.0 IGO. thepath.unaids.org/wp-content/themes/unaids2023/assets/files/2023_report.pdf

[3] Mashaphu, S., Wyatt, G.E., Zhang, M., Liu, H.: Condom use consistency among South African HIV serodiscordant couples following an HIV risk-reduction intervention. International journal of STD & AIDS **33**(5), 479–484 (2022). doi.org/10.1177/09564624221076617

[4] Pilowsky, D.J., Wu, L.T.: Sexual risk behaviors and HIV risk among americans aged 50 years or older: a review. Substance Abuse and Rehabilitation **6**, 51–60 (2015). doi.org/10.2147/sar.s78808

[5] Jourdain, H., Billioti de Gage, S., Desplas, D., Dray-Spira, R.: Real-world effectiveness of Pre-Exposure Prophylaxis in men at high risk of HIV infection in france: a nested case-control study. The Lancet **7**, 529–536 (2022). doi.org/10.1016/S2468-2667(22)00106-2

[6] Centers for Disease Control and Prevention. (2018). Pre-Exposure Prophylaxis for the prevention of HIV infection in the United States–2017 update. Atlanta, GA: CDC. www.cdc.gov/hiv/pdf/risk/prep/cdc-hiv-prep-guidelines-2021.pdf

[7] Fonner, V.A., Dalglish, S.L., Kennedy, C.E., Baggaley, R., O'Reilly, K., Koechlin, F., Rodolph, M., Hodges-Mameletzis, I., Grant, R.M.: Effectiveness and safety of oral HIV Pre-Exposure Prophylaxis for all populations. AIDS **30**(12), 1973–1983 (2016). doi.org/10.1097/qad.0000000000001145

[8] World Health Organization. (2012). Programmatic update: antiretroviral treatment as prevention (TasP) of HIV and TB: executive summary (No. WHO/HIV/2012.12). World Health Organization. iris.who.int/bitstream/handle/10665/70904/WHO_HIV_2012.12_eng.pdf

[9] Paschen-Wolff, M.M., Campbell, A.N.C., Tross, S. et al.: HIV Treatment Knowledge in the Context of “Treatment as Prevention” (TasP). AIDS Behav. **24**, 2984–2994 (2020). doi.org/10.1007/s10461-020-02849-7

[10] Hu, Q., Meyers, K., Xu, J., Chu, Z., Zhang, J., Ding, H., Han, X., Jiang, Y., Geng, W.: Efficacy and cost-effectiveness of early antiretroviral therapy and partners’ Pre-Exposure Prophylaxis among men who have sex with men in Shenyang,
China: a prospective cohort and costing study. BMC Infectious Diseases **19**(1) (2019). doi.org/10.1186/s12879-019-4275-x

[11] Zablotska, I., Grulich, A.E., Phanuphak, N., Anand, T., Janyam, S., Poonkaset-wattana, M., Baggaley, R., van Griensven, F. and Lo, Y.R.: PrEP implementation in the Asia-Pacific region: opportunities, implementation and barriers. Journal of the International AIDS Society **19**: 21119 (2016). doi.org/10.7448/IAS.19.7.21119

[12] Calabrese, S.K.: Understanding, Contextualizing, and Addressing PrEP Stigma to Enhance PrEP Implementation. Curr. HIV/AIDS Rep. **17**, 579–588 (2020). doi.org/10.1007/s11904-020-00533-y

[13] Wood, S., Gross, R., Shea, J.A. et al.: Barriers and Facilitators of PrEP Adherence for Young Men and Transgender Women of Color. AIDS Behav. **23**, 2719–2729 (2019). doi.org/10.1007/s10461-019-02502-y

[14] Anderson, R.M., Medley, G.F., May, R.M., Johnson, A.M.: A preliminary study of the transmission dynamics of the Human Immunodeficiency Virus (HIV), the causative agent of AIDS. Mathematical medicine and biology **3**(4), 229–263 (1986). doi.org/10.1093/imammb/3.4.229

[15] Perelson, A.S., Kirschner, D.E., De Boer, R.: Dynamics of hiv infection of CD4+ T-cells. Mathematical biosciences **114**(1), 81–125 (1993). doi.org/10.1016/0025-5564(93)90043-a

[16] Buratto, A., Cesaretto, R., Zamarchi, R.: Mathematical models for HIV treatment: A schematic review. Journal of Interdisciplinary Mathematics **23**(3), 707–725 (2020). doi.org/10.1080/09720502.2020.1732049

[17] Cassels. S, Clark, S.J., Morris, M.: Mathematical Models for HIV Transmission Dynamics: Tools for Social and Behavioral Science Research. JAIDS Journal of Acquired Immune Deficiency Syndromes **47**, 34–39 (2008). DOI:10.1097/QAI.0b013e3181605da3

[18] Isham, V.: Mathematical Modelling of the Transmission Dynamics of HIV Infection and AIDS: A Review. Journal of the Royal Statistical Society. Series A (Statistics in Society) **151**(1), 5–49 (1988). doi.org/10.2307/2982179

[19] Culshaw, R.V., Ruan, S.: A delay-differential equation model of HIV infection of CD4+ T-cells. Mathematical biosciences **165**(1), 27–39 (2000). doi.org/10.1016/s0025-5564(00)00006-7

[20] Nelson, P.B., Perelson, A.S.: Mathematical analysis of delay differential equation models of HIV-1 infection. Mathematical biosciences **179**(1), 73–94 (2002). doi.org/10.1016/s0025-5564(02)00099-8
[21] Wang, L., Din, A., Wu, P.: Dynamics and optimal control of a spatial diffusion HIV/AIDS model with antiretroviral therapy and Pre-Exposure Prophylaxis treatments. Math. Meth. Appl. Sci. 45(16), 10136-10161 (2022). doi:10.1002/mma.8359

[22] Wu, P., Ahmed, S., Wang, X., Wang, H.: PreEP Intervention in the Mitigation of HIV/AIDS Epidemics in China via a Data-Validated Age-Structured Model. Bull. Math. Biol. 85(5), 41 (2023). doi.org/10.1007/s11538-023-01145-4

[23] Afassinou, K., Chirove, F., Govinder, K.S.: Pre-Exposure Prophylaxis and antiretroviral treatment interventions with drug resistance. Mathematical biosciences 285, 92-101 (2017). doi.org/10.1016/j.mbs.2017.01.005

[24] Djomegni, P.M.T., Tekle, A., Dawed, M.Y.: Pre-Exposure Prophylaxis HIV/AIDS mathematical model with non classical isolation. Japan J. Indust. Appl. Math. 37, 781-801 (2020). doi.org/10.1007/s13160-020-00422-2

[25] Wang, L., Moqueet, N., Simkin, A., Knight, J., Ma, H., Lachowsky, N. J., Armstrong, H. L., Tan, D. H. S., Burchell, A. N., Hart, T. A., Moore, D. M., Adam, B. D., Macfadden, D. R., Baral, S., Mishra, S.: Mathematical modelling of the influence of serosorting on the population-level HIV transmission impact of Pre-Exposure Prophylaxis. AIDS (London, England), 35(7), 1113-1125 (2021). doi.org/10.1097/QAD.00000000000002826

[26] Chazuka, Z., Mudimu, E., Mathebula, D.: Stability and bifurcation analysis of an HIV model with Pre-Exposure Prophylaxis and treatment interventions. Scientific African 23, e01979, ISSN 2468-2276 (2024). doi.org/10.1016/j.sciaf.2023.e01979

[27] Kermack, W.O., McKendrick, A.G.: A contribution to the mathematical theory of epidemics. Proc. R. Soc. Lond. A 115(772), 700-721 (1927). doi.org/10.1098/rspa.1927.0118

[28] Adimy, M., Molina, J., Pujo-Menjouet, L., Ranson, G., Wu, J.: Forecasting the effect of Pre-Exposure Prophylaxis (PrEP) on HIV propagation with a system of differential-difference equations with delay. Mathematics 10(21), 4093 (2022). doi.org/10.3390/math10214093

[29] Capasso, V., Serio, G.: A generalization of the Kermack-McKendrick deterministic epidemic model. Mathematical Biosciences 42(1), 43-61 (1978). doi.org/10.1016/0025-5564(78)90006-8

[30] Adimy, M., Chekroun, A., Ferreira, C.P.: Global dynamics of a differential-difference system: a case of Kermack-McKendrick SIR model with age-structured protection phase. Mathematical Biosciences and Engineering 17(2), 1329-1354 (2020). doi:10.3934/mbe.2020067
[31] Dangerous inequalities: World AIDS Day report 2022. Geneva: Joint United Nations Programme on HIV/AIDS; 2022. Licence: CC BY-NC-SA 3.0 IGO www.unaids.org/sites/default/files/media_asset/dangerous-inequalities_en.pdf

[32] Adimy, M., Chekroun, A., Touaoula, T.M.: Age-structured and delay differential-difference model of hematopoietic stem cell dynamics. Discrete and Continuous Dynamical Systems - Series B 20(9), 2765–2791 (2015). doi:10.3934/dcdsb.2015.20.2765

[33] Cruz, M.A., Hale, J.K.: Stability of functional differential equations of neutral type. Journal of Differential Equations 7(2), 334–355 (1970). doi.org/10.1016/0022-0396(70)90114-2

[34] Smith, H.L., Thieme, H.R.: Dynamical Systems and Population Persistence. Graduate Studies in Mathematics, 118, American Mathematical Society, Providence R.I. (2011). doi.org/10.1090/gsm/118

[35] Freedman, H.I., Moson, P.: Persistence definitions and their connections. Proceedings of the American Mathematical Society 109, 1025–1033 (1990). doi:10.2307/2048133

[36] Fan, G., Maung Min, O.O., Wolkowicz Gail, S.K.: Hopf bifurcation of delay differential equations with delay dependent parameters. Canadian Appl. Math. Quart. 17, 37–60 (2009). math.mcmaster.ca/minoo/mypapers/FMW.pdf

[37] Adimy, M., Chekroun, A., Kuniya, T.: Global asymptotic stability for a distributed delay differential-difference system of a Kermack–McKendrick SIR model. Applicable Analysis 102(12), 3463–3475 (2023). doi.org/10.1080/00036811.2022.2075352

[38] Adimy, M., Chekroun, A., Dugourd-Camus, C., Meghelli, H.: Global Asymptotic Stability of a Hybrid Differential–Difference System Describing SIR and SIS Epidemic Models with a Protection Phase and a Nonlinear Force of Infection. Qual. Theory Dyn. Syst. 23, 34 (2024). doi.org/10.1007/s12346-023-00891-z