**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

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, *τ* = 3 months. Therefore, *a* ∈ [0, *τ*]. 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

$$P(t) = \int\_0^\tau p(t, a) da.$$

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 *σ*) 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 *β*. A part,

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

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

$$f(\mathcal{S}) = \mathcal{S}^n \quad \text{or} \quad f(\mathcal{S}) = \frac{\mathcal{S}^n}{1 + \mathcal{S}^n}, \quad 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*(*S*) = *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 *ψ* 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 *ψ* is

$$
\psi(t) = \psi\_0 + \frac{\kappa\_1 t^m}{t^m + \kappa\_2}, \quad 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 \mathbb{K} := \psi\_0 + \kappa\_1,
$$

to be a constant. We can also choose *t* → *ψ*(*t*) 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 *μ* is consequently the same for all). The system satisfied by *S* and *I* is then given, for all *t* > 0, by

$$\begin{cases} \begin{array}{rcl} S'(t) &=& \sigma - \beta I(t)S(t) - \mu S(t) - Kf(S(t)) + (1 - \theta)p(t, \tau), \\\ I'(t) &=& \beta I(t)S(t) - \mu I(t). \end{array} \end{cases} \tag{1}$$

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

$$\begin{cases} \frac{\partial p}{\partial t}(t,a) + \frac{\partial p}{\partial a}(t,a) = -\mu p(t,a), & 0 < a < \tau, \\\ p(t,0) = Kf(S(t)) + \theta p(t,\tau). \end{cases} \tag{2}$$

We remind here that *a* ∈ [0, *τ*] 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

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

The boundary condition

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

represents individuals who renewed their PrEP, *θp*(*t*, *τ*), or those who started it, *K f*(*S*(*t*)).

Using the characteristics method and the boundary condition from (2) (see [12]), we obtain, for *t* > 0 and *a* ∈ [0, *τ*],

$$p(t,a) = \begin{cases} \ e^{-\mu a} p(t-a,0), & t > a, \\\ e^{-\mu t} p(0,a-t) = e^{-\mu t} p\_0(a-t), & 0 \le t \le a. \end{cases} \tag{3}$$

To write Equation (3) in a convenient form, we supposed that the initial condition *p*<sup>0</sup> is continuous on [0, *τ*] and satisfies the compatibility condition

$$p\_0(0) = Kf(S\_0) + \theta p\_0(\tau).$$

We put

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

and we define the function

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

We then obtain from (3) the expression

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

This means that we can extend the function *u* to the interval [−*τ*, 0], by putting *u*(*t*) = *ϕ*(*t*) for *<sup>t</sup>* <sup>∈</sup> [−*τ*, 0]. Then, we have directly *<sup>p</sup>*(*t*, *<sup>τ</sup>*) = *<sup>e</sup>*−*μτu*(*<sup>t</sup>* <sup>−</sup> *<sup>τ</sup>*), for all *<sup>t</sup>* <sup>≥</sup> 0. Consequently, using the boundary condition of the system (2), *u* satisfies the difference equation:

$$u(t) = \begin{cases} Kf(S(t)) + \theta e^{-\mu \tau} u(t - \tau), & t > 0, \\ \varphi(t), & -\tau \le t \le 0. \end{cases}$$

Furthermore, the total population *P*(*t*) can be explicitly expressed, for *t* ≥ 0, in terms of the function *u*:

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

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

$$\begin{cases} S'(t) &= \sigma - \beta I(t)S(t) - \mu S(t) - Kf(S(t)) + (1 - \theta)e^{-\mu \tau}u(t - \tau), \\ \quad I'(t) &= \beta I(t)S(t) - \mu I(t), \\ u(t) &= \, Kf(S(t)) + \theta e^{-\mu \tau}u(t - \tau), \end{cases} \tag{4}$$

with the initial condition:

$$S(0) = S\_0, \quad I(0) = I\_0 \quad \text{and} \quad \mu(t) = \varphi(t), \quad t \in [-\tau, 0]. \tag{5}$$

Knowing that *P*(*t*) = *e*−*μ<sup>t</sup> t t*−*τ e μa u*(*a*)*da*, we obtain

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

We then have

$$N'(t) = S'(t) + I'(t) + P'(t) = \sigma - \mu N(t).$$

We easily obtain

$$N(t) = \left(N\_0 - \frac{\sigma}{\mu}\right)e^{-\mu t} + \frac{\sigma}{\mu} \to \frac{\sigma}{\mu'} \quad \text{as} \ t \to +\infty. \tag{6}$$
