# Bifurcation analysis and optimal control of an infection age-structured epidemic model with vaccination and treatment

Yves Fotso Fotso, Josué Tchouanti, Janvier Ntahomvukiye, Samuel Bowong

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

## To cite this version:

Yves Fotso Fotso, Josué Tchouanti, Janvier Ntahomvukiye, Samuel Bowong. Bifurcation analysis and optimal control of an infection age-structured epidemic model with vaccination and treatment. International Journal of Biomathematics, 2025, 10.1142/S1793524525500937. hal-05305076

HAL Id: hal-05305076

https://hal.science/hal-05305076v1

Submitted on 18 Feb 2026

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Bifurcation analysis and optimal control of an infection
age-structured epidemic model with vaccination and treatment

Yves Fotso Fotso<sup>1,2,6</sup>, Josué Tchouanti <sup>3</sup>, Janvier Ntahomvukiye<sup>4</sup>, Samuel Bowong<sup>5,7</sup>

¹ Institut Agro, Univ. Rennes, INRAE, IGEPP, Rennes, France

² Univ. Rennes, Inria, CNRS, IRISA, Rennes, France

³ LIFEWARE, INRIA Saclay, France

⁴ Department of Mathematics, University of Burundi, Burundi

⁵ Department of Mathematics and Computer Science, University of Douala, Cameroon

⁶ Department of Mathematics and Computer Science, University of Dschang, Cameroon

⁷ IRD, Sorbonne University, UMMISCO, Bondy, France

October 31, 2025

Abstract

This paper addresses the problem of bifurcation analysis and optimal control of an age-structured epidemic model. We first develop and analyze a deterministic epidemiological model for the transmission of an infectious disease described by a susceptible – vaccinated – infected (SVI). The model takes into account the age structure of the infection, imperfect vaccination, and treatment of infected people who become susceptible again. The model without intervention strategies is completely analyzed. We compute the disease-free equilibrium and derive the basic reproduction number $R_0$. The analysis of the model reveals the existence of the phenomenon of backward bifurcation, where a stable disease-free equilibrium coexists with one or more stable endemic equilibria when the associated basic reproduction number is less than unity. In this case, controlling this infectious disease remains extremely difficult. Therefore, based on this continuous model, a control is formulated and solved as an optimal control theory problem, indicating how control terms on vaccination and treatment should be introduced in the population to minimize the number of infected individuals as well as the cost of applying controls. We establish the first-order necessary conditions for optimality and characterize optimal controls. The theory is supported by numerical simulations, which further suggest that the combination of intervention measures can significantly reduce the transmission dynamics of infectious diseases by keeping the population of infected individuals relatively low.

**Keywords: Age-structured model, Infectious disease, Backward bifurcation theory, Optimal control theory, Numerical simulations.**

# 1 Introduction

Infectious diseases are one of the greatest enemies of humans in the world. Although many efforts have been made in the fight against these diseases, some of which have been eliminated, many continue to develop and are responsible for the loss of many human lives. These diseases are caused by numerous pathogens, including bacteria, viruses, parasites, and fungi, which can evolve and harm patients if not detected and treated early at the population level. A combination of rapid and appropriate control
strategies can greatly help decision-makers limit the damage caused by these diseases. Vaccination and treatment are among the most important control strategies for reducing the spread of many infectious diseases. However, due to strain variation, vaccination may lose its effect. Given the economic and human damage caused by these diseases, it is necessary to develop frameworks for simulating their spread in order to better understand and predict their behavior, with a view to controlling them efficiently and at lower cost. In this case, mathematical modeling is an important tool for providing valuable information in the public health decision-making process by generating several scenarios through numerical simulations, taking into account different assumptions. It also enables the evaluation of various intervention measures before deploying them in the population. It has been established that the infectivity of certain infectious diseases depends on the time elapsed since infection. This time is generally referred to as the infection age. This age is widely regarded as a potentially useful variable for understanding the transmission mechanisms of certain infectious diseases, such as HIV/AIDS [26, 22, 33], Influenza disease [37], Tuberculosis (TB) [1, 19, 41] and see also [40, 31, 34, 1, 25] and references therein for general class of epidemic models. Thus, the study of age-structured epidemic models could provide a better understanding of the transmission mechanism of these diseases.

Age-structured epidemic models have been extensively studied recently. The main research subjects focus on the well-posedness of solutions which include the nonnegativity, boundedness and asymptotic smoothness of solutions, the existence of global compact attractors, the calculation of basic reproduction number, the existence of equilibria and their local and global stability, the permanence and persistence of the disease, and the existence of Hopf bifurcation, etc. [30, 9, 39, 36, 21, 42, 40]. Recently Wang et al. [40] studied the global dynamics of an Susceptible–Vaccinated–Infectious–Recovered (SVIR) epidemiological model with continuous infection age structure and imperfect vaccination. They have shown that the disease-free equilibrium is globally asymptotically stable when the basic reproduction rate is less than or equal to one, while when this threshold is greater than one, there is a unique (endemic) equilibrium point that is globally asymptotically stable in a feasible region. However, they do not take into account the temporary induced immunity of vaccination and recovery.

The study of optimal control in age-structured epidemic models has a significant meaning in mathe-
matical biology, and it has also attracted many authors. The existence and uniqueness of the control
strategy in some age-structured epidemic models with optimal control have been proved using the
maximum principle. We refer the reader to Refs; [5, 24, 4] for the general optimal control theory of such
systems, and also [34, 1, 2, 14, 27, 13, 35, 20, 15, 17] for applications of this theory in many branches of
science such as competing species, harvesting control, birth control, epidemic disease and plant-pest
interactions. Since that age-structured model is formulated as a partial differential equation, it is not
easy to study the related optimal control problem.

In this work, we intend to fill some of the gaps mentioned above to better take into account the
fact that the vaccination does not confer total immunity and that the infected individuals become
susceptible again after receiving treatment and can be vaccinated again. From the modeling perspective,
the model proposed in this work is a substantial extension of the model presented and studied in [40] by
considering (i) the vaccine-induced immunity acquired by vaccinated individuals is supposed to wane
and does not confer a permanent immunity and (ii) the recovery does not confer a permanent immunity
to infected individuals so that the recovered class is no longer necessary. The well-posedness of the
model, including the positivity and the boundedness of solutions, and the existence and uniqueness of a
solution, is studied. We compute the disease-free equilibrium and derive the explicit expression of the
basic reproduction number using the next generation method [10, 36, 21] . The long-term dynamics
of the model are studied as a function of this threshold. In general, reducing the basic reproduction
number below unity is a necessary and sufficient condition for eradicating infectious diseases. However,
several authors (see for instance [42, 25, 18, 7] and references therein) have demonstrated that it is
possible that in certain epidemiological models, when the basic reproduction number is less than unity,
the disease-free stable equilibrium point coexists simultaneously with an endemic stable equilibrium
point. This phenomenon, known as backward bifurcation, reveals that to eradicate a disease, it's not sufficient to reduce the basic reproduction number below unity, but below a threshold value. Determining the causes of this bifurcation in epidemiological models has attracted the attention of many researchers [42, 25, 18, 7, 30]. Using the Lyapunov-Schmidt theory introduced recently by Martcheva and Inaba in [30], we show that the model can exhibit the phenomenon of backward bifurcation, and we determine its cause. After, we formulate the optimal control problem with specific control strategies and provide a proof of the existence of an optimal control function. The aim is to minimize the number of people infected and the cost of implementing vaccination and treatment. We derive an optimality system that characterizes the optimal control. Numerical results of the continuous optimal vaccination and therapy by solving the optimality system.

The remainder of this manuscript is organized as follows. In Section 2, we present the model and the theoretical results. Section 3 is devoted to formulating the optimal control problem and establishing the existence and characterization of optimal control. Numerical simulations are provided in Section 4 to illustrate the theoretical results, and we end the paper with a conclusion.

## 2 Model formulation and mathematical analysis

This section is devoted to the main results of this manuscript about the qualitative analysis of the formulated model. These include the well-posedness of the model, as well as the existence and stability of equilibria.

### 2.1 The model

We here consider an infectious disease in a closed age-structured host population and assume that the vaccination of susceptible individuals does not confer permanent immunity to recipients, and the recovery of infected individuals confers temporary immunity on the recovered individuals. The age-structured SVIS model is formulated by the following system of equations:

$$ egin{cases} \frac{dS}{dt} = \Lambda - \mathcal{T}(\beta i(t, .))S + \mathcal{T}(\alpha i(t, .)) - (\xi + \mu)S, \\ \frac{dV}{dt} = \xi S - m \mathcal{T}(\beta i(t, .))V - \mu V, \\ \partial_t i(t, \tau) + \partial_\tau i(t, \tau) = -[\mu + \delta(\tau) + \alpha(\tau)]i(t, \tau), \\ i(t, 0) = [S + mV]\mathcal{T}(\beta i(t, .)), \end{cases} \quad (2.1) $$

where the state variables $S(t)$ and $V(t)$ denote the size of susceptible and vaccinated individuals at time $t > 0$ respectively. The state variable $i(t, \tau)$ represents the density of infectious individuals with the infection age $\tau \in \mathbb{R}_+$ at time $t > 0$. All recruitment is into the susceptible class and occurs at a constant rate $\Lambda$. The susceptible individuals are vaccinated at a constant rate $\xi$ and move to the class of vaccinated at a flow $\xi S$ with a vaccine efficiency $1 - m$ with the parameter $m \in (0, 1)$. The transmission of the disease to susceptible and vaccinated individuals occurs due to their contact with infected individuals. Then, the corresponding force of infection are given by $\mathcal{T}(\beta i)$ and $m\mathcal{T}(\beta i)$ respectively, where $\mathcal{T}(.)$ is the integral operator defined for some integrable function $f$ on $\mathbb{R}_+$ by $\mathcal{T}(f) = \int_{\mathbb{R}_+} f(\tau)d\tau$ and $\beta(\tau)$ is the age-structured disease transmission rate. The natural and disease-induced death of individuals are given by the rates $\mu$ and $\delta(\tau)$ respectively. The model considers that the treatment do not confer a permanent immunity and the treat recipients return to the susceptible class $S$ at flow $\mathcal{T}(\alpha i)$ where $\alpha(\tau)$ represents the treatment rate of infectious individuals of infection age $x > 0$.

Before stating our results, we make the following realistic assumptions on the positivity of the parameters and functions involved in the system (2.1). We make the following assumption.

**Assumption 2.1.** The parameters $\Lambda, \xi, \mu, m$ are all positives and the initial conditions $S(0)$ and $V(0)$ are nonnegative. Moreover, the functions $\alpha(\cdot), \beta(\cdot)$ and $\delta(\cdot)$ belong in $L_+^\infty(0, \infty)$, and the boundary condition $i(0, \cdot) = i_0(\cdot) \in L_+^1(0, \infty)$. Here $L_+^p(0, \infty)$ with $1 \le p \le \infty$ denotes the space of nonnegative functions in $L^p(0, \infty)$.
## 2.2 Well-posedness

Since the system (2.1) has a nonlinear initial condition on the variable $i(t, \tau)$, we use the integrated semigroup theory introduced recently by Thieme [39] in the context of age-structured models to establish the well-posedness of the model. It consists of removing the nonlinearity from the domain and incorporating it into the Lipschitz continuous perturbation function. Let us introduce the Banach space $X = \mathbb{R} \times \mathbb{R} \times L^1(0, \infty) \times \mathbb{R}$ endowed with the usual product norm. The positive cone of space $X$ is defined by $X_+ = \mathbb{R}_+ \times \mathbb{R}_+ \times L^1(0, \infty) \times \mathbb{R}_+$. We set $X_0 = \mathbb{R} \times \mathbb{R} \times L^1(0, \infty) \times \{0\}$ and denote $X_{0+} = X_0 \cap X_+$. The system (2.1) can be rewritten in the following abstract Cauchy problem

$$ \frac{dz(t)}{dt} = \mathcal{A}z(t) + H(z(t)), \quad \forall t \ge 0 \quad \text{and} \quad z(0) = z_0 \in X_{0+}, \qquad (2.2) $$

where $z(t) = (S(t), V(t), i(t, .), 0)^{\top}$ with $^{\top}$ denoting the transposition symbol. The linear differential operator $\mathcal{A}: D(\mathcal{A}) \subset X \to X$ is defined as follows:

$$ \mathcal{A} \begin{pmatrix} S \\ V \\ i \\ 0 \end{pmatrix} = \begin{pmatrix} -\mu S \\ -(\mu + \xi)V \\ -i' - pi \\ -i(0) \end{pmatrix}, \qquad (2.3) $$

with $p(\tau) = \mu + \delta(\tau) + \alpha(\tau)$ and the domain $D(\mathcal{A}) = \mathbb{R} \times \mathbb{R} \times W^{1,1}(0, \infty) \times \{0\}$ where $W^{1,1}(0, \infty) = \{f \in L^1(0, \infty) : D^k f \in L^1(0, \infty), \forall |k| \le 1\}$ is a Sobolev space. Note that the domain $D(\mathcal{A})$ is not dense in $X$ since it can be seen that $X_0 = \overline{D(\mathcal{A})} \ne X$. Finally the nonlinear map $H: X_0 \subset X \to X$ is defined by

$$ H \begin{pmatrix} S \\ V \\ i \\ 0 \end{pmatrix} = \begin{pmatrix} \Lambda - \mathcal{T}(\beta i)S + \mathcal{T}(\alpha i) \\ \xi S - m \mathcal{T}(\beta i)V \\ 0_{L^1(0, \infty)} \\ [S + mV]\mathcal{T}(\beta i) \end{pmatrix}. \qquad (2.4) $$

Now, we are ready to establish the well-posedness of the system (2.1), which is given by the following result:

**Theorem 2.1.** The system (2.1) represented by the abstract Cauchy problem (2.2) has a unique strongly continuous semiflow $\{\Phi(t, .)\}_{t \ge 0}$ on $X_{0+}$ such that for each $z_0 \in X_{0+}$, the map $z(\cdot) \in C([0, \infty), X_{0+})$ defined by $z(\cdot) : t \to z(t) = \Phi(t, z_0)$ is an integrated (or mild) solution of system (2.2) satisfying $\int_0^t z(s)ds \in \mathcal{D}$ for all $t \ge 0$ and $z(t) = z_0 + \int_0^t z(s)ds + \int_0^t F(z(s))ds$. In addition, the non-empty domain

$$ \Omega = \{(S, V, i, 0) : S(t) + V(t) + \mathcal{T}(i(t, .)) \le \Lambda/\mu\}, \qquad (2.5) $$

is positively invariant and attracts all nonnegative solutions. Moreover, the semiflow $\{\Phi(t, .)\}_{t \ge 0}$ is bounded dissipative, that is, there exists a bounded set $B \subset X_0$ such that for any bounded set $U \subset X_0$, we can find $\kappa = \kappa(U, B)$ such that $\Phi(t, U) \subset B$ for all $t \ge \kappa$.

*Proof.* The existence, uniqueness, and positivity of the integrated solution of system (2.2) can be obtained by applying a similar approach as in [28, 31, 32]. Let $z_0 \in X_{0+}$, then $\|\Phi(t, z_0)\|_X = S(t) + V(t) + \mathcal{T}(i(t, .))$. By integrating the third equation of system (2.1) on $\mathbb{R}_+$ with respect to $x$ and combining with the first–second equation of (2.1), one can easily get

$$ \frac{d\|\Phi(t, z_0)\|_X}{dt} \le \Lambda - \mu\|\Phi(t, z_0)\|_X. \qquad (2.6) $$

Then, using the Grönwall–Bellman inequality, we have

$$ \|\Phi(t, z_0)\|_X \le \Lambda/\mu - (\Lambda/\mu - \|z_0\|_X) e^{-\mu t}, \qquad (2.7) $$
which shows that $\Phi(t, z_0) \in \Omega$ holds for every solution of (2.2) satisfying $z_0 \in \Omega$. Hence, the set $\Omega$ is positively invariant. Furthermore, we have the bound $\limsup_{t \to +\infty} \|\Phi(t, z_0)\|_X \le \Lambda/\mu$ which implies that the semiflow $\{\Phi(t, z_0)\}_{t \ge 0}$ is bounded dissipative and $\Omega$ attracts all the points of space $X_{0+}$. $\square$

## 2.3 Basic reproduction number and equilibria

System (2.1) has always a disease-free equilibrium point $E_0 = (S^0, V^0, 0, 0) \in X_{0+}$ where $S^0 = \Lambda/(\mu + \xi)$ and $V^0 = \xi S^0/\mu$, corresponding to the equilibrium point without disease. To study the long-term behavior of the system (2.1), we need to calculate the basic reproduction number denoted $\mathcal{R}_0$. It is obtained by the so-called next-generation operator, which gives the distribution of secondary infections as a function of the distribution of the primary infected individuals. To calculate the expression of $\mathcal{R}_0$, we use the methodology developed in references [10, 21] where $\mathcal{R}_0$ corresponds to the spectral radius of the next-generation operator. Specifically, we linearize system (2.1) around the disease-free equilibrium point $E_0$ and obtain the following equations for the dynamics of the infected population:

$$ \begin{cases} \partial_t i(t, \tau) + \partial_\tau i(t, \tau) = -p(\tau)i(t, \tau), \\ i(t, 0) = [S^0 + mV^0]\mathcal{T}(\beta i). \end{cases} \tag{2.8} $$

Using the characteristics method, the solution of system (2.1) can be expressed as

$$ i(t, \tau) = i_0(\tau - t) \frac{\pi(\tau)}{\pi(\tau - t)} \mathbb{1}_{\{\tau > t\}} + i(t - \tau, 0) \pi(\tau) \mathbb{1}_{\{\tau < t\}}, \tag{2.9} $$

where the function $\pi(\tau) = e^{-\int_0^\tau p(\sigma)d\sigma}$ is the survival rate at which an individual stays in the infectious compartment after some time $\tau \ge 0$. Let $\bar{h}(t) = i(t, 0)$ be the number of newly infectious individuals at time $t > 0$. Inserting expression (2.9) into the boundary condition of (2.8), we get the renewal equation:

$$ \bar{h}(t) = \varphi(t) + \int_0^t \Psi(\tau)\bar{h}(t - \tau)d\tau, \tag{2.10} $$

where $\varphi(t) := [S^0 + mV^0] \int_t^\infty \beta(\tau) \frac{\pi(\tau)}{\pi(\tau - t)} i_0(\tau - t) d\tau$ and $\Psi(\tau) := [S^0 + mV^0]\beta(\tau)\pi(\tau)$. According to [10, 21], the basic reproduction number is calculated as the spectral radius of the next-generation operator $\int_0^\infty \Psi(\tau) d\tau$. Hence, the explicit expression of $\mathcal{R}_0$ is given by

$$ \mathcal{R}_0 = [S^0 + mV^0]\mathcal{T}(\beta\pi). \tag{2.11} $$

This threshold depends on the epidemiological parameters of the model, which ensures or not the outbreak of an epidemic process and measures the expected number of secondary cases produced by a typical infected individual during its entire period of infectiousness in a susceptible population [8, 10, 36]. However, this threshold is used to measure the transmission potential of an infectious disease. In general, it plays an important role in the determination and stability of equilibria. We have the following result:

**Theorem 2.2.** Let Assumption 2.1 be satisfied. If $\mathcal{R}_0 < 1$, the disease-free equilibrium point $E_0$ is locally asymptotically stable and it is unstable when $\mathcal{R}_0 > 1$.

*Proof.* Linearizing the system (2.1) at equilibrium point $E_0$ with $y(t), z(t)$ and $\omega(t, \tau)$ being the small perturbations, that is $y(t) = S(t) - S^0, z(t) = V(t) - V^0$ and $\omega(t, \tau) = i(t, \tau)$. We obtain the following abstract Cauchy problem:

$$ \frac{du(t)}{dt} = \mathcal{A}u(t) + DH[E_0](u(t)), \quad \forall t \ge 0 \quad \text{and} \quad u(0) \in X_{0+}, \tag{2.12} $$
where $u(t) = (y(t), z(t), \omega(t, .), 0)^T$ and the linear operator $DH[E_0] : X_0 \subset X \to X$ is defined for $u \in X_0$ by

$$DH[E_0](u) = \begin{pmatrix} -\mathcal{T}(\beta\omega)S^0 + \mathcal{T}(\alpha\omega) \\ \xi y - m\mathcal{T}(\beta\omega)V^0 \\ 0_{L^1(0,\infty)} \\ [S^0 + mV^0]\mathcal{T}(\beta\omega) \end{pmatrix}. \qquad (2.13)$$

Let us denote by $A_0$ the part of the linear operator $\mathcal{A}$ in $X_0$, i.e. $A_0: \psi \in X_0 \to A_0\psi = \mathcal{A}\psi \in X_0$ and let $\{T_{A_0}(t)\}_{t \ge 0}$ the semigroup generated by the linear operator $A_0$. It is easy to check by adapting the proof of Proposition 1 of [16] that $\|T_{A_0}(t)\| \le e^{-\mu t}$ for all $t \ge 0$. It follows that $\omega_{ess}(A_0)$, the essential growth rate of $\{T_{A_0}(t)\}_{t \ge 0}$ is less than or equal to $-\mu$. Let $\{T_{(A+DHE_0)_0}(t)\}_{t \ge 0}$ be the semigroup generated by $(A+DHE_0)_0$, the part of the linear operator $A+DHE_0$ in $X_0$. Since $DH[E_0]$ is a compact operator, it follows from the result in [11, Theorem 1.2] that $\omega_{ess}((A+DHE_0)_0) \le \omega_{ess}(A_0) \le -\mu < 0$. According to the result obtained in [39, Corollary 4.3], the equilibrium $E_0$ is locally asymptotically stable if all eigenvalues of the linear operator $(A+DH[E_0])_0$ have negative real part. In this case, the trajectories that start sufficiently close to $E_0$ remain close and converge to the equilibrium when time tends towards infinity. However, if at least one eigenvalue of $(A+DH[E_0])_0$ has a strictly positive real part, then $E_0$ is an unstable equilibrium. Now, we consider the exponential solutions of the linearized system (2.12) at the disease-free equilibrium $E_0$ by $y(t) = \bar{y}e^{\lambda t}$, $z(t) = \bar{z}e^{\lambda t}$ and $\omega(t) = \bar{\omega}(\tau)e^{\lambda t}$ with $(\bar{y}, \bar{z}, \bar{\omega}(\cdot)) \in \mathbb{R} \times \mathbb{R} \times W(0, \infty) \setminus \{0\}$ and $\lambda \in \mathbb{C}$ with $\mathfrak{Re}\lambda > -\mu$ (here the symbol $\mathfrak{Re}$ denotes the real part) to derive the characteristic equation. We get the following linear eigenvalue problem:

$$\left\{ \begin{array}{l} (\lambda + \mu + \xi)\bar{y} = -S^0\mathcal{T}(\beta\bar{\omega}) + \mathcal{T}(\alpha\bar{\omega}), \\ (\lambda + \mu)\bar{z} = \xi\bar{y} - mV^0\mathcal{T}(\beta\bar{\omega}), \\ \bar{\omega}'(\tau) = -[\lambda + p(\tau)]\bar{\omega}(\tau), \\ \bar{\omega}(0) = [S^0 + mV^0]\mathcal{T}(\beta\bar{\omega}). \end{array} \right. \qquad (2.14)$$

Since $\bar{y}$ and $\bar{z}$ do not interact on $\bar{\omega}$-equation, we can determine $\lambda$ follows the three-four equations of (2.14). From the third equation of (2.14), we get $\bar{\omega}(\tau) = \bar{\omega}(0)\pi(\tau)e^{-\lambda\tau}$. Putting this expression in the last equation of (2.14) and canceling $\bar{\omega}(0)$, we obtain the following characteristic equation:

$$G(\lambda) := [S^0 + mV^0] \int_0^\infty \beta(\tau)\pi(\tau)e^{-\lambda\tau}d\tau = 1. \qquad (2.15)$$

Assume that $\mathcal{R}_0 < 1$ and $\lambda \in \mathbb{C}$ with $\mathfrak{Re}\lambda \ge 0$, we have $1 = |G(\lambda)| \le G(\mathfrak{Re}\lambda) \le G(0) = \mathcal{R}_0 < 1$ which is a contradiction. Hence, the equation $G(\lambda) = 1$ does not have a root with a nonnegative real part when $\mathcal{R}_0 < 1$. Thus, we conclude that the disease-free equilibrium point $E_0$ is locally asymptotically stable whenever $\mathcal{R}_0 < 1$. Now, assume that $\mathcal{R}_0 > 1$ and $\lambda \in \mathbb{R}_+$ then $G(0) = \mathcal{R}_0 > 1$ and $\lim_{\lambda \to +\infty} G(\lambda) = 0$. Since $G(\cdot)$ is a decreasing function, there exists a unique $\lambda_0 > 0$ such that $G(\lambda_0) = 1$. Therefore $E_0$ is unstable whenever $\mathcal{R}_0 > 1$. This completes the proof. $\square$

Biologically speaking, this result means that the disease can be eradicated (when $\mathcal{R}_0 < 1$) when initial sizes of each subpopulation in the system (2.1) are within the basin of attraction of the stable point $E_0$. We now examine the existence of the endemic equilibrium points. For this, let $E^*(\tau) = (S^*, V^*, i^*(\tau), 0)$ be any arbitrary equilibrium point of system (2.1). To find conditions for the existence of equilibrium points for which disease is endemic in the population we set the derivatives concerning time in (2.2) equal to zero (i.e. by solving the abstract algebraic equation $\mathcal{A}E^* + H(E^*) = 0$), that is

$$\left\{ \begin{array}{l} \Lambda - \mathcal{T}(\beta i^*)S^* + \mathcal{T}(\alpha i^*) - (\xi + \mu)S^* = 0, \\ \xi S^* - m\mathcal{T}(\beta i^*)V^* - \mu V^* = 0, \\ \partial_{\tau}i^*(\tau) = -[\mu + \delta(\tau) + \alpha(\tau)]i^*(\tau), \\ i^*(0) = [S^* + mV^*]\mathcal{T}(\beta i^*). \end{array} \right. \qquad (2.16)$$
Then, the equilibrium point $E^*$ with positive components is such that

$$S^* = \frac{\Lambda + \mathcal{T}(\alpha i^*)}{\mathcal{T}(\beta i^*) + \mu + \xi}, \quad V^* = \frac{\xi S^*}{m \mathcal{T}(\beta i^*) + \mu}, \quad i^*(\tau) = i^*(0)\pi(\tau),$$

and $i^*(0)$ is the positive real solution of the quadratic equation

$$a_2(i^*(0))^2 + a_1i^*(0) + a_0 = 0, \qquad (2.17)$$

where

$$\begin{aligned} a_2 &= m[\mathcal{T}(\beta\pi)]^2[1 - \mathcal{T}(\alpha\pi)], \quad a_0 = \mu(\mu + \xi)(1 - \mathcal{R}_0), \\ a_1 &= m\mathcal{T}(\beta\pi)[\mu - \Lambda\mathcal{T}(\beta\pi)] + (\mu + m\xi)\mathcal{T}(\beta\pi)[1 - \mathcal{T}(\alpha\pi)]. \end{aligned}$$

We see that the coefficients $a_2$ and $a_0$ are positives (respectively negatives) if and only if $\mathcal{T}(\alpha\pi) < 1$ and $\mathcal{R}_0 < 1$ (respectively $\mathcal{T}(\alpha\pi) > 1$ and $\mathcal{R}_0 > 1$) respectively. According to Descartes' rule of signs, the existence of the possible positive real roots of equation (2.17) is summarized in Table 1. It follows that under the condition $\mathcal{R}_0 < 1$, it is possible to have one or two endemic equilibria. In this case, it is possible to have a backward bifurcation phenomenon in the system (2.1) i.e. the locally stable disease-free equilibrium $E_0$ coexists with a locally stable endemic equilibrium when $\mathcal{R}_0 < 1$. To check this, the discriminant of quadratic equation (2.17) given by $a_1^2 - 4a_2a_0$ is set to zero and solved for the critical value of $\mathcal{R}_0$ denoted $\mathcal{R}_0^\#$ that is

$$\mathcal{R}_0^\# = 1 - \frac{a_1^2}{4a_2\mu(\mu + \xi)}. \qquad (2.18)$$

Thus, the backward bifurcation phenomenon would occur for the values of threshold $\mathcal{R}_0$ satisfying the relation $\mathcal{R}_0^\# < \mathcal{R}_0 < 1$. It would be very interesting to analyze the system's bifurcations to see if it is possible to have a bistability phenomenon, and if so, to determine its cause.

<table><caption>Table 1: Number of possible positive real roots of the equation (2.17) according to the sign of the coefficients $a_k$, $k = 0, 1, 2$.</caption><thead><tr><th>a<sub>2</sub></th><th>a<sub>1</sub></th><th>a<sub>0</sub></th><th>R<sub>0</sub></th><th>Number of positive solutions of (2.17)</th></tr></thead><tbody><tr><td>-</td><td>-</td><td>-</td><td>R<sub>0</sub> > 1</td><td>0</td></tr><tr><td>+</td><td>-</td><td>-</td><td>R<sub>0</sub> > 1</td><td>1</td></tr><tr><td>+</td><td>+</td><td>-</td><td>R<sub>0</sub> > 1</td><td>1</td></tr><tr><td>-</td><td>+</td><td>-</td><td>R<sub>0</sub> > 1</td><td>0 or 2</td></tr><tr><td>-</td><td>-</td><td>+</td><td>R<sub>0</sub> < 1</td><td>1</td></tr><tr><td>-</td><td>+</td><td>+</td><td>R<sub>0</sub> < 1</td><td>1</td></tr><tr><td>+</td><td>+</td><td>+</td><td>R<sub>0</sub> < 1</td><td>0</td></tr><tr><td>+</td><td>-</td><td>+</td><td>R<sub>0</sub> < 1</td><td>0 or 2</td></tr></tbody></table>

## 2.4 Bifurcation analysis

Herein, we study the existence of bifurcations of system (2.1) around the disease-free equilibrium $E_0$ when $\mathcal{R}_0 = 1$. To this end, we use a recent result introduced by Martcheva and Inaba [30] to detect the presence of bifurcations in partial differential equations and give the necessary and sufficient
conditions for them to occur. Let us set $\beta(\tau) = \bar{\beta}\beta_0(\tau)$ where the function $\beta_0(\tau)$ is normalized as $[S^0 + mV^0]\mathcal{T}(\beta_0\pi) = 1$, which suggests that $\mathcal{R}_0 = 1$ is equivalent to $\bar{\beta} = 1$. In what follows, we consider $\bar{\beta}$ as the bifurcation parameter. Let us introduce the important threshold

$$ \Theta := - \left[ \left( 1 + m \frac{\xi}{\mu} \right) \frac{S^0}{\xi + \mu} + \frac{m^2}{\mu} V^0 \right] \mathcal{T}(\beta_0 \pi) + \left( 1 + m \frac{\xi}{\mu} \right) \frac{\mathcal{T}(\alpha \pi)}{\xi + \mu}. \quad (2.19) $$

Then, we have the following main result:

**Theorem 2.3.** Let Assumption 2.1 be satisfied. If $\Theta > 0$, the system (2.1) undergoes a backward bifurcation at $(E_0, \bar{\beta} = 1)$. Otherwise, the system (2.1) exhibits a forward bifurcation at $(E_0, \bar{\beta} = 1)$ when $\Theta < 0$.

*Proof.* It is based on the Lyapunov-Schmidt method. For this purpose, let us set $z = (z_1, z_2, z_3, 0) \in X_0$ and introduce the nonlinear map $F(z, \bar{\beta}) := \mathcal{A}z + H(z)$ where the linear operator $\mathcal{A}$ and map $H(z)$ are defined in (2.3) and (2.4) respectively. The linearized operator $\mathcal{B} := D_z F(E_0, \bar{\beta})$ acting on $X$ is given by

$$ \mathcal{B}z = \begin{pmatrix} -\bar{\beta}S^0\mathcal{T}(\beta_0z_3) + \mathcal{T}(\alpha z_3) - (\xi + \mu)z_1 \\ \xi z_1 - \bar{\beta}mV^0\mathcal{T}(\beta_0z_3) - \mu z_2 \\ -z_3' - p(\tau)z_3 \\ -z_3(0) + \bar{\beta}[S^0 + mV^0]\mathcal{T}(\beta_0z_3) \end{pmatrix}, \quad (2.20) $$

where its domain $D(\mathcal{B}) = D(\mathcal{A})$. Let us define by $\lambda \in \mathbb{C}$ the eigenvalue of the operator $\mathcal{B}$. Solving the third equation of $\mathcal{B}z = \lambda z$ for $z \in X_0$, we obtain $z_3(\tau) = z_3(0)\pi(\tau)e^{-\lambda\tau}$. Inserting the expression of $z_3(\tau)$ into the four equation in $\mathcal{B}z = \lambda z$ and canceling $z_3(0)$, we get the characteristic equation

$$ \bar{\beta}[S^0 + mV^0] \int_0^\infty \beta_0(\tau)\pi(\tau)e^{-\lambda\tau}d\tau = 1. \quad (2.21) $$

Note that the equation (2.21) has $\lambda = 0$ as a simple isolated eigenvalue when $\bar{\beta} = 1$, and each of the remaining solutions has a negative real part thanks to the proof of Theorem 2.2. To find the eigenvector $\hat{v}_0$ associated with eigenvalue zero, we need to solve the equation $\mathcal{B}z = 0$. We can choose $z_3(\tau) = \pi(\tau)$ and the boundary condition of the four equations of $\mathcal{B}z = 0$ is trivially satisfied at $\bar{\beta} = 1$. From the first equation of $\mathcal{B}z = 0$, we get

$$ z_1^0 := z_1 = \frac{-\bar{\beta}S^0\mathcal{T}(\beta_0\pi) + \mathcal{T}(\alpha\pi)}{\xi + \mu}, \quad (2.22) $$

and the second equation gives

$$ z_2^0 := z_2 = -\frac{\bar{\beta}m}{\mu}V^0\mathcal{T}(\beta_0\pi) + \frac{\xi}{\mu}z_1^0. \quad (2.23) $$

Therefore, we set the vector $\hat{v}_0 = (z_1^0, z_2^0, \pi(\tau), 1)$. Now, we seek the adjoint operator $\mathcal{B}^*$. To do so, let us consider the vector $\Psi = (\Psi_1, \Psi_2, \Psi_3, \Psi_4) \in X^* := \mathbb{R}^2 \times L^1(0, \infty) \times \mathbb{R}$. Then, we have the relation

$$ \langle \Psi, \mathcal{B}z \rangle = [-\bar{\beta}S^0\mathcal{T}(\beta_0z_3) + \mathcal{T}(\alpha z_3) - (\xi + \mu)z_1] \Psi_1 + [\xi z_1 - \bar{\beta}mV^0\mathcal{T}(\beta_0z_3) - \mu z_2] \Psi_2 \\ + \int_0^\infty [-z_3' - p(\tau)z_3] \Psi_3(\tau)d\tau + [-z_3(0) + \bar{\beta}[S^0 + mV^0]\mathcal{T}(\beta_0z_3)] \Psi_4. $$

Integrating by part, assuming $\Psi_3(\infty) = 0$, we obtain

$$ \int_0^\infty [-z_3' - p(\tau)z_3] \Psi_3(\tau)d\tau = z_3(0)\Psi_3(0) + \int_0^\infty [\Psi_3' - p(\tau)\Psi_3] z_3(\tau)d\tau. $$
Hence, we have

$$ \begin{aligned} \langle \Psi, \mathcal{B}z \rangle &= [\xi \Psi_2 - (\xi + \mu)\Psi_1] z_1 - \mu \Psi_2 z_2 + (\Psi_3(0) - \Psi_4)z_3(0) \\ & \quad + \int_0^\infty [\Psi'_3 - p\Psi_3 + (-\bar{\beta}S^0\beta_0 + \alpha)\Psi_1 - \bar{\beta}mV^0\beta_0\Psi_2 + \bar{\beta}(S^0 + mV^0)\beta_0\Psi_4]z_3(\tau)d\tau \\ &= \langle \mathcal{B}^*\Psi, z \rangle, \end{aligned} $$

which should hold for $\Psi \in D(\mathcal{B}^*)$ and $z \in D(\mathcal{B}) \subset X_0$. We take the domain of the operator $\mathcal{B}^*$ as follows

$$ D(\mathcal{B}^*) = \{ \Psi \in \mathbb{R}^2 \times W^{1,\infty}(0, \infty) \times \mathbb{R} : \Psi_3(\infty) = 0, \Psi_4 = \Psi_3(0) \}. $$

Therefore, the adjoint operator $\mathcal{B}^*$ is defined for $\Psi \in D(\mathcal{B}^*) \subset X^* = \mathbb{R}^2 \times L^\infty(0, \infty) \times \mathbb{R}$ by

$$ \mathcal{B}^* \Psi = \begin{pmatrix} \xi \Psi_2 - (\xi + \mu) \Psi_1 \\ -\mu \Psi_2 \\ \Psi'_3 - p \Psi_3 + (-\bar{\beta} S^0 \beta_0 + \alpha) \Psi_1 - \bar{\beta} m V^0 \beta_0 \Psi_2 + \bar{\beta} [S^0 + m V^0] \beta_0 \Psi_4 \\ \Psi_3(0) \end{pmatrix}. \quad (2.24) $$

To find the vector $\hat{v}_0^*$, the unique positive vector satisfying the relation $\langle \mathcal{B}z, \hat{v}_0^* \rangle = 0$ for all $z \in X$, it suffices to solve the equation $\mathcal{B}^*\Psi = 0$ with $\bar{\beta} = 1$. It is easy to see that $\Psi_1 = \Psi_2 = 0$. Solving the differential equation in $\mathcal{B}^*\Psi = 0$, we get $\Psi_3(\tau) = [S^0 + mV^0]\Psi_3(0) \int_\tau^\infty \beta_0(s) \frac{\pi(s)}{\pi(\tau)}ds$. Hence, the eigenfunction $\hat{v}_0^* = (0, 0, \Psi_3(\tau), \Psi_3(0))$ where $\Psi_3(0) > 0$. Computing the second derivative for $F$ with respect to $z$ and $\bar{\beta}$, we obtain

$$ D_{z\bar{\beta}}^2 F(E_0, 1) \hat{v}_0 = \begin{pmatrix} -S^0 \mathcal{T}(\beta_0 \pi) z_1^0 \\ -m V^0 \mathcal{T}(\beta_0 \pi) z_2^0 \\ 0 \\ [S^0 + m V^0] \mathcal{T}(\beta_0 \pi) \end{pmatrix}. \quad (2.25) $$

Therefore, one obtains

$$ b = \langle D_{z\bar{\beta}}^2 F(E_0, 1) \hat{v}_0, \hat{v}_0^* \rangle = [S^0 + mV^0] \Psi_3(0) \mathcal{T}(\beta_0 \pi) > 0. \quad (2.26) $$

On the order hand, the second derivative $D_{zz}^2 F(E_0, 1)[\hat{v}_0, \hat{v}_0]$ can be computed with the following formula:

$$ D_{zz}^2 F(E_0, 1)[\hat{v}_0, \hat{v}_0] = \partial_{hk}^2 F(E_0 + (h+k)\hat{v}_0, 1) = \begin{pmatrix} -2\mathcal{T}(\beta_0\pi)z_1^0 \\ -2m\mathcal{T}(\beta_0\pi)z_2^0 \\ 0 \\ 2[z_1^0 + m z_2^0]\mathcal{T}(\beta_0\pi) \end{pmatrix}. \quad (2.27) $$

Thus, we have

$$ a = \langle D_{zz}^2 F(E_0, 1)[\hat{v}_0, \hat{v}_0], \hat{v}_0^* \rangle = 2[z_1^0 + m z_2^0] \mathcal{T}(\beta\pi) \Psi_3(0) = 2\Theta \mathcal{T}(\beta\pi) \Psi_3(0). \quad (2.28) $$

It follows that $a < 0$ (respectively $a > 0$) if and only if $\Theta < 0$ (respectively $\Theta > 0$). According to [30, Theorem 2.1] (see Appendix A), the bifurcation is backward if and only if $\Theta > 0$ and forward if and only if $\Theta < 0$. This completes the proof. $\square$

The appearance of a backward bifurcation in the epidemiological model (2.1) can lead to the persistence of disease in the population even if the basic reproduction number is less than unity. In such a scenario, the condition $\mathcal{R}_0 < 1$ becomes only a necessary but not sufficient condition for disease eradication. This important property of the model is caused by the treatment rate of infectious individuals since $\Theta < 0$ when $\alpha(\tau) = 0$ for all $\tau \in \mathbb{R}_+$. In this case, the coefficients of the quadratic equation (2.17) $a_2 > 0$ and $a_0 < 0$ when $\mathcal{R}_0 > 1$. Consequently, the system (2.1) admits a unique endemic equilibrium. Moreover, the global asymptotic stability of the disease-free equilibrium of the system (2.1) is given below for the special case $\alpha(\tau) = 0$.
**Theorem 2.4.** Assume that the therapeutic treatment is not effective (i.e. $\alpha(\tau) = 0$ for all $\tau \in \mathbb{R}_+$). Then, the system (2.1) always admits a disease-free equilibrium point $E_0$ which is globally asymptotically stable whenever $\mathcal{R}_0 \le 1$.

*Proof.* To establish the global stability of the disease-free equilibrium $E_0$, we use the Lyapunov function approach. For this, we consider the function $q(\tau) = \int_{\tau}^{\infty} \beta(s) \frac{\pi(s)}{\pi(\tau)} ds$. Note that $q \in L_{+}^{\infty}(0, \infty)$, $q(0) = \mathcal{T}(\beta\pi)$ and its satisfies the following differential equation

$$q'(\tau) - (\mu + \delta(\tau))q(\tau) + \beta(\tau) = 0. \qquad (2.29)$$

Let us introduce the Volterra-type function, which takes the form $h(z) = z - \ln z - 1$. Obviously, $h'(z) = 1 - (1/z)$ and $h(z) \ge 0$ for all $z > 0$. Thus $h(.)$ is decreasing function on $(0, 1)$ and increasing on $(1, +\infty)$. In addition to the fact that $h''(z) > 0$ for all $z > 0$, then $h(.)$ has a unique global minimum at $z_0 = 1$ satisfying $h(z_0) = 0$. Let us consider the following Lyapunov function:

$$W[t] = S_0 h \left( \frac{S(t)}{S^0} \right) + V^0 h \left( \frac{V(t)}{V_0} \right) + [S^0 + mV^0] \mathcal{T}(qi(t, .)). \qquad (2.30)$$

It follows that the function $W[.]$ is nonnegative and well defined concerning the disease-free equilibrium $E_0$, which is a global minimum. Differentiating the function $W[.]$ along the solution of the system (2.1), we have

$$\begin{aligned} \frac{dW[t]}{dt} &= \left(1 - \frac{S^0}{S}\right) \frac{dS}{dt} + \left(1 - \frac{V^0}{V}\right) \frac{dV}{dt} + [S^0 + mV^0] \int_0^\infty q(\tau) \partial_t i(t, \tau) d\tau \\ &= \left(1 - \frac{S^0}{S}\right) \left(\Lambda - \mathcal{T}(\beta i)S + \mathcal{T}(\alpha i) - (\xi + \mu)S\right) + \left(1 - \frac{V^0}{V}\right) \left(\xi S - m \mathcal{T}(\beta i)V - \mu V\right) \\ &\quad + [S^0 + mV^0] \int_0^\infty q(\tau) \left(-\partial_\tau i(t, \tau) - (\mu + \delta(\tau))\right) d\tau. \end{aligned}$$

By using the fact that $\Lambda = (\xi + \mu)S^0$ and $\xi S^0 = \mu V^0$, we get after a few algebraic calculations

$$\begin{aligned} \frac{dW[t]}{dt} ={}& \mu S^0 \left(2 - \frac{S^0}{S} - \frac{S}{S^0}\right) + \xi S^0 \left(3 - \frac{V}{V^0} - \frac{S^0}{S} - \frac{S}{S_0} \frac{V^0}{V}\right) + [S^0 + mV^0] \mathcal{T}(\beta i) \\ & - [S + mV] \mathcal{T}(\beta i) + [S^0 + mV^0] \int_0^\infty q(\tau) \left(-\partial_\tau i(t, \tau) - (\mu + \delta(\tau))\right) d\tau. \end{aligned} \qquad (2.31)$$

Using the integration by parts formula, we have

$$\int_0^\infty q(\tau) \partial_\tau i(t, \tau) d\tau = [q(\tau)i(t, \tau)]_0^\infty - \int_0^\infty q'(\tau) i(t, \tau) d\tau = -q(0)i(t, 0) - \int_0^\infty q'(\tau) i(t, \tau) d\tau. \qquad (2.32)$$

Combining the relations (2.31) and (2.32), we have after calculation

$$\begin{aligned} \frac{dW[t]}{dt} ={}& \mu S^0 \left(2 - \frac{S^0}{S} - \frac{S}{S^0}\right) + \xi S^0 \left(3 - \frac{V}{V^0} - \frac{S^0}{S} - \frac{S}{S_0} \frac{V^0}{V}\right) + (1 - \mathcal{R}_0)i(t, 0) \\ & + [S^0 + mV^0] \int_0^\infty [q'(\tau) - (\mu + \delta(\tau))q(\tau) + \beta(\tau)]i(t, \tau)d\tau. \end{aligned} \qquad (2.33)$$

Finally, we get,

$$\frac{dW[t]}{dt} = \mu S^0 \left(2 - \frac{S^0}{S} - \frac{S}{S^0}\right) + \xi S^0 \left(3 - \frac{V}{V^0} - \frac{S^0}{S} - \frac{S}{S_0} \frac{V^0}{V}\right) + (1 - \mathcal{R}_0)i(t, 0). \qquad (2.34)$$
According to the arithmetic–geometric mean inequality, both terms $2 - \frac{S^0}{S} - \frac{S^0}{S^0}$ and $3 - \frac{V}{V^0} - \frac{S^0}{S} - \frac{S^0}{S^0} - \frac{V^0}{V}$ are nonnegatives and the equality holds if and only if $S = S^0$ and $V = V^0$. Therefore, the function $W[.]$ is nonincreasing since $\frac{dW[t]}{dt} \le 0$ when $\mathcal{R}_0 \le 1$. By the boundedness of function $W[.]$, the alpha limit set of $(S(t), V(t), i(t, .))$ must be contained in the maximal compact invariant subset defined by $\Gamma = \{(S, V, i) : \frac{dW[t]}{dt} = 0\}$. However, the equality $\frac{dW[t]}{dt} = 0$ holds if and only if $S = S^0$, $V = V^0$ and $i(t, 0) = 0$. Integrating the $i(t, \tau)$–equation of system (2.1) using the characteristics method, it follows that $i(t, \tau) = 0$ for all $t > x$. Hence, we get $i(t, \tau) \to 0$ as $t \to +\infty$. Therefore, the set $\Gamma$ is reduced to a singleton $\{E_0\}$ and using the Lasalle’s invariant principle [23], we can conclude that the disease–free equilibrium $E_0$ is globally asymptotically stable when $\mathcal{R}_0 \le 1$. $\square$

Epidemiologically speaking, Theorem 2.4 means that in the absence of treatment, a flow of infectious individuals for all ages of infection will not generate an outbreak of the disease unless $\mathcal{R}_0 > 1$ where the disease will persist. Therefore, reducing the basic reproduction number below unity becomes a necessary and sufficient condition to eradicate the disease in the population. Treatment of infected individuals can lead to non–non–eradication of the disease when the basic reproduction rate is less than one. It is therefore important to develop an optimal combined control strategy aimed at reducing the spread of the disease through appropriate vaccination and treatment protocols at the lowest possible cost.

# 3 Optimal control Problem

## 3.1 Control problem statement

Given the seriousness of the damage caused by infectious diseases in many countries, it is important to know how much and when vaccination and treatment should be applied, in order to control the dynamics of disease transmission effectively and at a minimum cost. For this purpose, we consider that the vaccination and treatment are represented by Lebesgue measurable functions on finite time horizon $[0, T]$, denoted $\xi(t)$ and $\alpha(t, \tau)$ respectively and where $T > 0$ is the final time of intervention strategies. In addition, we replace the upper limit of the integral concerning the infection age by a finite number $\tau_{max} > 0$ for practical consideration, i.e. the integral function $\mathcal{T}(f) = \int_0^{\tau_{max}} f(\tau)d\tau$. Let us introduce the set $Q = [0, T] \times [0, \tau_{max}]$. Given the limited resources and time available to implement strategies to control this infectious disease, policies must be limited to a predefined objective. For this reason, we consider by $\mathcal{D} := \mathcal{L}^\infty(0, T; [0, \xi_+]) \times \mathcal{L}^\infty(Q; [0, \alpha_+])$ the set of admissible control pair which is the space of measurable and bounded functions pair $(\xi(\cdot), \alpha(\cdot, \cdot))$ defined by $\xi(\cdot) : [0, T] \to [0, \xi_+]$ and $\alpha(\cdot, \cdot) : Q \to [0, \alpha_+]$. The positive constants $\xi_+$ and $\alpha_+$ represent the maximum rates at which individuals may be vaccinated and treated, respectively. Hence, implementing both time-dependent controls in the system (2.1), we obtain the following controlled system:

$$ 
\left\{ 
\begin{array}{l}
\displaystyle \frac{dS}{dt} = \Lambda - \mathcal{T}(\beta i(t, \cdot))S + \mathcal{T}(\alpha(t, \cdot)i(t, \cdot)) - (\xi(t) + \mu)S, \\
\displaystyle \frac{dV}{dt} = \xi(t)S - m\mathcal{T}(\beta i(t, \cdot))V - \mu V, \\
\displaystyle \partial_t i(t, \tau) + \partial_\tau i(t, \tau) = -[\mu + \delta(\tau) + \alpha(t, \tau)]i(t, \tau), \\
\displaystyle i(t, 0) = [S + mV]\mathcal{T}(\beta i(t, \cdot)).
\end{array}
\right.
\qquad (3.1) $$

Our main objective is to minimize the total number of infected individuals and the necessary cost of vaccination and treatment. To achieve this goal, we work together with the system (3.1), the following objective functional:

$$ \mathcal{J}(\xi, \alpha) = \int_Q \chi(\tau)i(t, \tau)dtd\tau + \frac{C_1}{2}\|\xi\|_{L^2(0,T)}^2 + \frac{C_2}{2}\|\alpha\|_{L^2(Q)}^2, \qquad (3.2) $$

where $\chi(\cdot) \in L_+^\infty(0, \tau_{max})$, $C_1 > 0$ and $C_2 > 0$ are balancing coefficients transforming the integral into a cost spent over the interval $[0, T]$. Quadratic expression of the control in (3.2) is included to
indicate the nonlinearity of the implementation cost, as it is more costly to increase the control efficiency when it is already high. As mentioned before, our optimal control problem reads as follows: find an admissible control pair $(\xi^*(.), \alpha^*(., .)) \in \mathcal{D}$ steering the optimal trajectory $(S^*(.), V^*(.), i^*(., .))$ satisfying the optimization problem

$$ \mathcal{J}(\xi^*, \alpha^*) = \min_{(\xi, \alpha) \in \mathcal{D}} \mathcal{J}(\xi, \alpha). \qquad (\text{OCP}) $$

## 3.2 First-order necessary optimality conditions

Here, we start by showing the existence of an optimal solution to the optimization problem (OCP) subject to the age-structured system (3.1).

**Theorem 3.1.** Under Assumption 2.1, there exists at least one optimal control pair $(\xi^*(.), \alpha^*(.)) \in \mathcal{D}$ at which corresponds the state variable $(S^*(.), V^*(.), i^*(., a))$, solution of the optimization problem (OCP).

*Proof.* The objective functional $\mathcal{J}(\xi, \alpha) \ge 0$ since the state variables and controls are all nonnegative. Then, it follows that $d = \inf\{\mathcal{J}(\xi, \alpha) : (\xi, \alpha) \in \mathcal{D}\}$ is finite and nonnegative. So, there is a minimizing sequence $(\xi_n, \alpha_n) \in \mathcal{D}$ such that for $n \ge 1$ we have

$$ d \le \mathcal{J}(\xi_n, \alpha_n) \le d + \frac{1}{n}. \qquad (3.3) $$

As the sequence $(\xi_n, \alpha_n)$ is bounded, there exists a subsequence still denoted $(\xi_n, \alpha_n)$ that converges to some $(\xi^*, \alpha^*)$ for the weak-\star topology of $L^\infty(0, T) \times L^\infty(Q)$. The set $\mathcal{D}$ is a closed convex subset of $L^\infty(0, T) \times L^\infty(Q)$, so it is weakly closed. Therefore, $(\xi^*, \alpha^*) \in \mathcal{D}$. Denote by $(S_n, V_n, i_n)$ the solution of system (3.1) corresponding to control pair $(\xi_n, \alpha_n)$. The sequence $(S_n, V_n)$ is uniformly bounded and equicontinuous on $[0, T]$. Using Arzela-Ascoli's Theorem, we can extract a subsequence still denoted $(S_n, V_n)$ which converges uniformly to the limit $(S^*, V^*)$ in $\mathcal{C}(0, T)$. Let us denote by

$$ \pi_{\alpha_n}(t, \tau) = e^{-\int_0^\tau [\mu + \delta(s) + \alpha_n(t - \tau + s, s)] ds}. \qquad (3.4) $$

It is easy to see that the function $\pi_{\alpha_n}$ is Lipschitz in the following sense, for $(\alpha^1, \alpha^2)$, there exists a constant $k \ge 0$ such that

$$ |\pi_{\alpha^1}(t, \tau) - \pi_{\alpha^2}(t, \tau)| \le k \mathcal{T}(|\alpha^1(t, \cdot) - \alpha^2(t, \cdot)|). \qquad (3.5) $$

As consequence, we have the convergence $\pi_{\alpha_n}(t, \tau) \to \pi_{\alpha^*}(t, \tau)$ as $n \to \infty$ almost everywhere in $Q$. By using the method of characteristics, we get the explicit expression of $i_n(t, \tau)$-equation

$$ i_n(t, \tau) = i_0(\tau - t) \frac{\pi_{\alpha_n}(t, \tau)}{\pi_{\alpha_n}(t, \tau - t)} \mathbb{1}_{\{\tau>t\}} + i_n(t - \tau, 0) \pi_{\alpha_n}(t, \tau) \mathbb{1}_{\{\tau<t\}}. \qquad (3.6) $$

This sequence is bounded since the sequences $(S_n)$, $(V_n)$, $(\mathcal{T}(i_n))$ and $(\alpha_n)$ are all bounded. Then, we can extract a subsequence still denoted $(i_n(t, \tau))$ that converges weakly to $i^*(t, \tau)$ in $L^2(Q)$ defined as follows:

$$ i^*(t, \tau) = i_0(\tau - t) \frac{\pi_{\alpha^*}(t, \tau)}{\pi_{\alpha^*}(t, \tau - t)} \mathbb{1}_{\{\tau>t\}} + i^*(t - \tau, 0) \pi_{\alpha^*}(t, \tau) \mathbb{1}_{\{\tau<t\}}. \qquad (3.7) $$

Since the sequence $(\mathcal{T}(i_n))$ is bounded, it converges to $\mathcal{T}(i^*)$ by the uniqueness of the limit. Moreover, passing to the limit in the differential equations satisfied by the subsequences $(S_n)$ and $(V_n)$ in the controlled system (3.1), we obtain:

$$ \left\{ \begin{array}{l} \frac{dS}{dt} = \Lambda - \mathcal{T}(\beta i^*) S^* + \mathcal{T}(\alpha(t, \cdot) i^*) - (\xi^*(t) + \mu)S^*, \\ \frac{dV}{dt} = \xi^*(t)S^* - m\mathcal{T}(\beta i^*)V^* - \mu V^*. \end{array} \right. \qquad (3.8) $$

Passing to the limit as $n \to \infty$ in (3.3) and lower semicontinuity of objective functional $\mathcal{J}(., .)$, we obtain $\lim_{n \to \infty} \mathcal{J}(\xi_n, \alpha_n) = \mathcal{J}(\xi^*, \alpha^*) = d$. Hence, $((S^*, V^*, i^*), (\xi^*, \alpha^*))$ is an optimal solution of the optimization problem (OCP). This achieves the proof. $\square$
Now, we derive the first-order necessary conditions for optimality of a control pair must satisfy to be optimal. Using the framework of Barbu [6], we derive the optimal control as a combination of the state and adjoint variables. The adjoint variables are determined by first introducing the sensitivity functions. For this purpose, let us denote by $\Psi := (S, V, i)$ and define the solution map $(\xi, \alpha) \mapsto \Psi(\xi, \alpha)$ of the system (3.1). The sensitivity functions $\varphi_S(t), \varphi_V(t)$ and $\varphi_i(t, \tau)$ associated to state variables $S(t), V(t)$ and $i(t, \tau)$, are obtained by the Gâteaux derivatives $\varepsilon^{-1}[\Psi((\xi, \alpha) + \varepsilon(q, p)) - \Psi(\xi, \alpha)] \to (\varphi_S, \varphi_V, \varphi_i)$ as $\varepsilon \to 0^+$ in $(L^\infty(0, T))^2 \times L^\infty(Q)$, where $(\xi, \alpha) + \varepsilon(q, p) \in \mathcal{D}$ and $(q, p) \in L^\infty(0, T) \times L^\infty(Q)$. By using the explicit expression of state variables of the system (3.1), it is easy to establish that the map $(\xi, \alpha) \mapsto \Psi(\xi, \alpha)$ is Lipschitz in $L^\infty$. Then, according to the result in Barbu [6, page 17], the existence of sensitivity functions is guaranteed. In addition, the sensitivity functions satisfy the following system of differential equations:

$$ 
\left\{
\begin{aligned}
\frac{d\varphi_S}{dt} &= -\mathcal{T}(\beta\varphi_i) + \mathcal{T}(\alpha\varphi_i) + \mathcal{T}(pi) - qS - [\mathcal{T}(\beta i) + \xi + \mu]\varphi_S, \\
\frac{d\varphi_V}{dt} &= \xi\varphi_S + qS - m\mathcal{T}(\beta i)\varphi_V - m\mathcal{T}(\beta\varphi_i)V - \mu\varphi_V, \\
\partial_t\varphi_i + \partial_\tau\varphi_i &= -[\mu + \delta + \alpha]\varphi_i - pi, \\
\varphi_i(t, 0) &= [S + mV]\mathcal{T}(\beta\varphi_i) + [\varphi_S + m\varphi_V]\mathcal{T}(\beta i).
\end{aligned}
\right.
\quad (3.9) $$

Next, we introduce the adjoint variables $\mathcal{Z}_S(t), \mathcal{Z}_V(t)$ and $\mathcal{Z}_i(t, \tau)$ corresponding to the state variables $S(t), V(t)$ and $i(t, \tau)$ of system (3.1), respectively. The adjoint system satisfied by the adjoint variables is then derived by using the adjoint operator associated with the sensitivity system (3.9) together with appropriate transversality and boundary conditions (see, for instance, [24] for more details). Specifically, the adjoint system is given by:

$$ 
\left\{
\begin{aligned}
\frac{d\mathcal{Z}_S}{dt} &= [\mathcal{T}(\beta i) + \mu]\mathcal{Z}_S + \xi(t)[\mathcal{Z}_S - \mathcal{Z}_V] - \mathcal{Z}_i(t, 0)\mathcal{T}(\beta i), \\
\frac{d\mathcal{Z}_V}{dt} &= [m\mathcal{T}(\beta i) + \mu]\mathcal{Z}_V - m\mathcal{T}(\beta i)\mathcal{Z}_i(t, 0), \\
\partial_t\mathcal{Z}_i + \partial_\tau\mathcal{Z}_i &= -\chi(\tau) + [\beta(\tau)S - \alpha(t, \tau)]\mathcal{Z}_S + m\beta(\tau)\mathcal{Z}_V V \\
&\quad -[S + mV]\beta(\tau)\mathcal{Z}_i(t, 0) + [\mu + \delta(\tau) + \alpha(t, \tau)]\mathcal{Z}_i,
\end{aligned}
\right.
\quad (3.10) $$

endowed with the following transversality and boundary conditions for $(t, \tau) \in Q$:

$$ \mathcal{Z}_S(T) = 0, \quad \mathcal{Z}_V(T) = 0, \quad \mathcal{Z}_i(T, \tau) = 0, \quad \mathcal{Z}_i(t, \tau_{\max}) = 0. \quad (3.11) $$

The existence of the adjoint solutions can be proved thanks to a fixed point argument mapping principle, see for instance [6]. In what follows, we employ tangent normal cone techniques in nonlinear functional analysis (see [6, 3, 4]) to deduce the first-order necessary conditions for optimality of optimal control. To do so, let us denote by $\mathcal{T}_{\mathcal{D}}(\xi, \alpha)$ and $\mathcal{N}_{\mathcal{D}}(\xi, \alpha)$, the tangent and normal cone of space $\mathcal{D}$ at control pair $(\xi, \alpha)$ respectively. We have the following result:

**Theorem 3.2.** Let $(\xi, \alpha) \in \mathcal{D}$ be an optimal admissible solution of the optimization problem (OCP) to which the trajectory $(S(t), V(t), i(t, \tau))$ is associated, and let $(\mathcal{Z}_S(t), \mathcal{Z}_V(t), \mathcal{Z}_i(t, \tau))$ be an adjoint vector satisfying the system (3.10)-(3.11). Then, we have the following optimal control structure:

$$ \xi(t) = \mathcal{P}_1 \left( \frac{[\mathcal{Z}_S(t) - \mathcal{Z}_V(t)]S(t)}{C_1} \right) \quad \text{and} \quad \alpha(t, \tau) = \mathcal{P}_2 \left( \frac{[-\mathcal{Z}_S(t) + \mathcal{Z}_i(t, \tau)]i(t, \tau)}{C_2} \right), \quad (3.12) $$

where the projection maps $\mathcal{P}_j(x) = \min\{\max\{0, x\}, x_{j \max}\}$ for $j = 1, 2$ with $x_{1 \max} = \xi_+$ and $x_{2 \max} = \alpha_+$.

*Proof.* For any element of the tangent cone $(q, p) \in \mathcal{T}_{\mathcal{D}}(\xi, \alpha)$ and let us denote by $(\xi^\varepsilon, \alpha^\varepsilon) := (\xi, \alpha) + \varepsilon(q, p)$, then we have $(\xi^\varepsilon, \alpha^\varepsilon) \in \mathcal{D}$ for $\varepsilon > 0$ small enough. Let $(S^\varepsilon, V^\varepsilon, i^\varepsilon)$ be the solution of the optimization problem (OCP) corresponding to the control pair $(\xi^\varepsilon, \alpha^\varepsilon)$. Since $\mathcal{J}(.,.)$ is minimal on $(\xi, \alpha)$, it follows that $\mathcal{J}(\xi^\varepsilon, \alpha^\varepsilon) \ge \mathcal{J}(\xi, \alpha)$ and then

$$ \int_Q \chi(\tau) (i^\varepsilon(t, \tau) - i(t, x)) dtd\tau + \frac{C_1}{2} \int_0^T (\xi^\varepsilon(t) - \xi^2(t)) dt + \frac{C_2}{2} \int_Q (\alpha^\varepsilon(t, \tau) - \alpha^2(t, \tau)) dtd\tau \ge 0. \quad (3.13) $$
Passing to the limit as $\varepsilon \to 0^+$ in the inequality above, we obtain

$$ \int_Q \chi(\tau)\varphi_i(t, \tau)dtd\tau + C_1 \int_0^T q(t)\xi(t)dt + C_2 \int_Q p(t, \tau)\alpha(t, \tau)dtd\tau \geq 0. \quad (3.14) $$

Let us multiply the first three equations in (3.9) by the variables $Z_S(t)$, $Z_V(t)$, $Z_i(t, \tau)$ and the equations in (3.10) by the variables $\varphi_S(t)$, $\varphi_V(t)$ and $\varphi_i(t, \tau)$, respectively. Then, we obtain, thanks to the relationship between the sensitive and adjoint operators (see, for instance, [24]), the relation

$$ \varphi_S \frac{dZ_S}{dt} + \varphi_V \frac{dZ_V}{dt} + \int_0^{\tau_{\max}} \varphi_i (\partial_t Z_i + \partial_\tau Z_i) d\tau + Z_S \frac{d\varphi_S}{dt} + Z_V \frac{d\varphi_V}{dt} + \int_0^{\tau_{\max}} Z_i (\partial_t \varphi_i + \partial_\tau \varphi_i) d\tau = 0. $$

By expanding the above equation and using integration by parts, the initial and boundary conditions, we obtain the following relationship

$$ \int_0^{\tau_{\max}} \chi(\tau)\varphi_i(t, \tau)dtd\tau - \int_0^{\tau_{\max}} [Z_S(t) - Z_i(t, \tau)] p(t, \tau)i(t, \tau)d\tau + [Z_S(t) - Z_V(t)] q(t)S(t) = 0. \quad (3.15) $$

Integrating (3.15) over $[0, T]$, we deduce that

$$ \int_Q \chi(\tau)\varphi_i(t, \tau)dtd\tau = \int_Q [Z_S(t) - Z_i(t, \tau)] p(t, \tau)i(t, \tau)dtd\tau + \int_0^T [-Z_S(t) + Z_V(t)] q(t)S(t)dt. \quad (3.16) $$

Consequently, inequality (3.2) becomes

$$ \int_Q [(-Z_S + Z_i)i - C_2\alpha]p(t, \tau)dtd\tau + \int_0^T [(Z_S - Z_V)S - C_1\xi]q(t)dt \leq 0 \quad (3.17) $$

for all $(q, p) \in \mathcal{T}_D(\xi, \alpha)$. Hence, according to the structure of normal cone (see, for instance, [6, 3] for more details) it follows that the expression $([Z_S - Z_V]S - C_1\xi, [-Z_S + Z_i]i - C_2\alpha) \in \mathcal{N}_D(\xi, \alpha)$. By the standard optimality arguments<fnref n="1" /> and taking into account the boundaries of each control strategies, we obtain the representation given in (3.12) which gives the desired result. This achieves the proof. □

**Remark 3.1.** At the final time $T > 0$, the optimal controls defined in (3.12) vanish, that are $\xi(T) = 0$ and $\alpha(T, \tau) = 0$ for all $\tau \in [0, \tau_{\max}]$ since the adjoint state variables are all equal to zero at the final time $T$ by (3.11).

Note that Theorem 3.1 does not guarantee the uniqueness of the optimal control pair. However, using the optimal control structure (3.12), this uniqueness can be obtained using the standard procedure based on Ekeland's variational principle [6, 12]. This principle is used in order to generate a sequence of controls and their corresponding states that converge to the optimal control and its corresponding state via the use of the convergence of a minimizing sequence of approximate objective functionals. More precisely, we have the following result

**Theorem 3.3.** There is a positive real constant $K_T$ that depends on the final time $T$ such that for $K_T < 1$, the optimization problem (OCP) has a unique optimal control pair $(\xi, \alpha) \in \mathcal{D}$ characterized by (3.12).

<footnote n="1">¹Refer to Barbu [6], the normal cone $\mathcal{N}_D(\xi, \alpha)$ is characterized as follows: for $(\omega_1, \omega_2) \in \mathcal{N}_D(\xi, \alpha)$, we have $\omega_1 > 0$ (resp. $\omega_2 > 0$) when $\xi = \xi_+$ (resp. $\alpha = \alpha_+$); $\omega_1 < 0$ (resp. $\omega_2 < 0$) when $\xi = 0$ (resp. $\alpha = 0$) and finally $\omega_1 = 0$ (resp. $\omega_2 = 0$) when $\xi \in [0, \xi_+]$ (resp. $\alpha \in [0, \alpha_+]$).</footnote>
Before starting the proof of the Theorem above, we embed the objective functional $\mathcal{J}(\xi, \alpha)$ into $L^1(0, T) \times L^1(Q)$ by defining the following functional:

$$ \hbar(\xi, \alpha) = \begin{cases} \mathcal{J}(\xi, \alpha), & \text{if } (\xi, \alpha) \in \mathcal{D} \\ +\infty, & \text{otherwise.} \end{cases} \qquad (3.18) $$

Let us introduce the following technical lemma that we shall use to establish the uniqueness result of the optimal control pair of (OCP), whose proof is described in Appendix B.

**Lemme 3.1.** We have the following properties:

1. For $T > 0$ sufficiently small, there are positive constants $\mathcal{C}_{1T}$ and $\mathcal{C}_{2T}$ such that:

(i) The map $(\xi, \alpha) \in \mathcal{D} \to (S, V, i)$ is Lipschitz in the following ways:

$$ \|(S^1, V^1, i^1) - (S^2, V^2, i^2)\|_{L^\infty} \le \mathcal{C}_{1T} \|(\xi^1, \alpha^1) - (\xi^2, \alpha^2)\|_{L^\infty}, $$

where $(S^k, V^k, i^k)$ is a solution of (3.1) associated to control pair $(\xi^k, \alpha^k)$ for $k \in \{1, 2\}$.

(ii) For $(\xi^k, \alpha^k) \in \mathcal{D}$, the adjoint system (3.10) admits a weak solution $(\mathcal{Z}_S^k, \mathcal{Z}_V^k, \mathcal{Z}_i^k) \in L^\infty(0, T) \times L^\infty(0, T) \times L^\infty(Q)$ such that for $k \in \{1, 2\}$ we have

$$ \|(\mathcal{Z}_S^1, \mathcal{Z}_V^1, \mathcal{Z}_i^1) - (\mathcal{Z}_S^2, \mathcal{Z}_V^2, \mathcal{Z}_i^2)\|_{L^\infty} \le \mathcal{C}_{2T} \|(\xi^1, \alpha^1) - (\xi^2, \alpha^2)\|_{L^\infty}. $$

2. The functional $\hbar(.,.)$ is lower semi-continuous with respect to $(\xi, \alpha)$ in $L^1(0, T) \times L^1(Q)$.

*Proof.* (of Theorem 3.3) Since the functional $\hbar(.,.)$ is lower semi-continuous in $L^1$, then according to Ekeland's variational principle [12, 6], it follows that for any given $\varepsilon > 0$, there exists a control pair $(\xi_\varepsilon, \alpha_\varepsilon) \in L^1(0, T) \times L^1(Q)$ such that:

$$ (i) \quad \hbar(\xi_\varepsilon, \alpha_\varepsilon) \le \inf_{(\xi, \alpha) \in \mathcal{D}} \hbar(\xi, \alpha) + \varepsilon, $$

$$ (ii) \quad \hbar(\xi_\varepsilon, \alpha_\varepsilon) = \inf_{(\xi, \alpha) \in \mathcal{D}} \{\hbar(\xi, \alpha) + \sqrt{\varepsilon} \|(\xi, \alpha) - (\xi_\varepsilon, \alpha_\varepsilon)\|_{L^1}\}. $$

Note that, the perturbed functional $\hbar_\varepsilon(\xi, \alpha) = \hbar(\xi, \alpha) + \sqrt{\varepsilon} \|(\xi, \alpha) - (\xi_\varepsilon, \alpha_\varepsilon)\|_{L^1}$ attains its infimum at the control pair $(\xi_\varepsilon, \alpha_\varepsilon)$. From item (ii) and a similar argument as that in Subsection 3.2 give the characterization of control pair $(\xi_\varepsilon, \alpha_\varepsilon)$ by

$$ \xi_\varepsilon(t) = \mathcal{P}_1 \left( \frac{[\mathcal{Z}_S^\varepsilon - \mathcal{Z}_V^\varepsilon]S^\varepsilon - \sqrt{\varepsilon}\vartheta_1^\varepsilon}{C_1} \right) \quad \text{and} \quad \alpha_\varepsilon(t, \tau) = \mathcal{P}_2 \left( \frac{[-\mathcal{Z}_S^\varepsilon + \mathcal{Z}_i^\varepsilon]i^\varepsilon - \sqrt{\varepsilon}\vartheta_2^\varepsilon}{C_2} \right), \qquad (3.19) $$

where $(S^\varepsilon, V^\varepsilon, i^\varepsilon)$ and $(\mathcal{Z}_S^\varepsilon, \mathcal{Z}_V^\varepsilon, \mathcal{Z}_i^\varepsilon)$ are the solutions of controlled and adjoint system respectively, corresponding to the control pair $(\xi_\varepsilon, \alpha_\varepsilon)$. The functions $\vartheta_1^\varepsilon \in L^\infty(0, T)$, $\vartheta_2^\varepsilon \in L^\infty(Q)$ and $|\vartheta_j^\varepsilon| \le 1$ for $j \in \{1, 2\}$ and $(t, \tau) \in Q$. We are now ready to prove our main result of this subsection concerning the uniqueness of an optimal control pair solution of the optimization problem (OCP). To do so, let us consider the following map $\mathcal{M}: \mathcal{D} \to \mathcal{D}$ defined by

$$ \mathcal{M}(\xi, \alpha) = \left( \mathcal{P}_1 \left( \frac{[\mathcal{Z}_S - \mathcal{Z}_V]S}{C_1} \right), \mathcal{P}_2 \left( \frac{[-\mathcal{Z}_S + \mathcal{Z}_i]i}{C_2} \right) \right). $$

Let $(S^k, V^k, i^k)$ and $(\mathcal{Z}_S^k, \mathcal{Z}_V^k, \mathcal{Z}_i^k)$ are the state and adjoint variables corresponding to the control pair $(\xi^k, \alpha^k)$ with $k \in \{1, 2\}$. Using the Lipschitz properties of the state and adjoint variable established in
Lemma 3.1, we have

$$
\begin{aligned}
\|\mathcal{M}(\xi^1, \alpha^1) - \mathcal{M}(\xi^2, \alpha^2)\|_{L^\infty} &= \left\| \mathcal{P}_1 \left( \frac{[\mathcal{Z}_S^1 - \mathcal{Z}_V^1]S^1}{C_1} \right) - \mathcal{P}_1 \left( \frac{[\mathcal{Z}_S^2 - \mathcal{Z}_V^2]S^2}{C_1} \right) \right\|_{L^\infty(0, T)} \\
&\quad + \left\| \mathcal{P}_2 \left( \frac{[-\mathcal{Z}_S^1 + \mathcal{Z}_i^1]i^1}{C_2} \right) - \mathcal{P}_2 \left( \frac{[-\mathcal{Z}_S^2 + \mathcal{Z}_i^2]i^2}{C_2} \right) \right\|_{L^\infty(Q)} \\
&\le C_1^{-1} \|[\mathcal{Z}_S^1 - \mathcal{Z}_V^1]S^1 - [\mathcal{Z}_S^2 - \mathcal{Z}_V^2]S^2\|_{L^\infty(0,T)} \\
&\quad + C_2^{-1} \|[[-\mathcal{Z}_S^1 + \mathcal{Z}_i^1]i^1 - [-\mathcal{Z}_S^2 + \mathcal{Z}_i^2]i^2\|_{L^\infty(Q)} \\
&\le K_1 \|(S^1, i^1) - (S^2, i^2)\|_{L^\infty} + K_2 \|(\mathcal{Z}_S^1, \mathcal{Z}_V^1, \mathcal{Z}_i^1) - (\mathcal{Z}_S^2, \mathcal{Z}_V^2, \mathcal{Z}_i^2)\|_{L^\infty} \\
&\le (K_1 \mathcal{C}_{1T} + K_2 \mathcal{C}_{2T}) (\|\xi^1 - \xi^2\|_{L^\infty(0,T)} + \|\alpha^1 - \alpha^2\|_{L^\infty(Q)}) \\
&\le \mathcal{C}_{3T} (\|(\xi^1, \alpha^1) - (\xi^2, \alpha^2)\|_{L^\infty}),
\end{aligned}
$$

where the constants $K_1$ and $K_2$ depend on the $L^\infty$ bounds on the state and adjoint state variables and

$$ \mathcal{C}_{3T} = K_1 \mathcal{C}_{1T} + K_2 \mathcal{C}_{2T}, \qquad (3.20) $$

where $\mathcal{C}_{1T}$ and $\mathcal{C}_{2T}$ are the Lipschitz constants obtained in Lemma 3.1. Clearly $\mathcal{M}$ is a contraction function if $\mathcal{C}_{3T} < 1$. Hence, $\mathcal{M}$ has a unique fixed point $(\xi^*, \alpha^*) \in \mathcal{D}$ by the Banach contraction theorem when $\mathcal{C}_{3T} < 1$. We show that this fixed point is an optimal control pair by using the approximating minimizers sequence $(\xi_\epsilon, \alpha_\epsilon)$ from Ekeland's principle and corresponding state variables $S^\epsilon, V^\epsilon$ and $i^\epsilon$, and adjoint variables $\mathcal{Z}_S^\epsilon, \mathcal{Z}_V^\epsilon$ and $\mathcal{Z}_i^\epsilon$. From Lemma 3.1 and the contraction property of the application $\mathcal{M}$, we have

$$
\begin{aligned}
\|(\xi^*, \alpha^*) - (\xi_\epsilon, \alpha_\epsilon)\|_{L^\infty} &= \left\| \mathcal{M}(\xi^*, \alpha^*) - \left( \mathcal{P}_1 \left( \frac{[\mathcal{Z}_S^\epsilon - \mathcal{Z}_V^\epsilon]S^\epsilon - \sqrt{\epsilon}\vartheta_1^\epsilon}{C_1} \right), \mathcal{P}_2 \left( \frac{[-\mathcal{Z}_S^\epsilon + \mathcal{Z}_i^\epsilon]i^\epsilon - \sqrt{\epsilon}\vartheta_2^\epsilon}{C_2} \right) \right) \right\|_{L^\infty} \\
&\le \|\mathcal{M}(\xi^*, \alpha^*) - \mathcal{M}(\xi_\epsilon, \alpha_\epsilon)\|_{L^\infty} \\
&\quad + \left\| \mathcal{M}(\xi_\epsilon, \alpha_\epsilon) - \left( \mathcal{P}_1 \left( \frac{[\mathcal{Z}_S^\epsilon - \mathcal{Z}_V^\epsilon]S^\epsilon - \sqrt{\epsilon}\vartheta_1^\epsilon}{C_1} \right), \mathcal{P}_2 \left( \frac{[-\mathcal{Z}_S^\epsilon + \mathcal{Z}_i^\epsilon]i^\epsilon - \sqrt{\epsilon}\vartheta_2^\epsilon}{C_2} \right) \right) \right\|_{L^\infty} \\
&\le \mathcal{C}_{3T} \|(\xi^*, \alpha^*) - (\xi_\epsilon, \alpha_\epsilon)\|_{L^\infty} + \sqrt{\epsilon}(C_1^{-1} + C_2^{-1}).
\end{aligned}
$$

Then for $\mathcal{C}_{3T} < 1$, we obtain

$$ \|(\xi^*, \alpha^*) - (\xi_\epsilon, \alpha_\epsilon)\|_{L^\infty} \le \frac{\sqrt{\epsilon}(C_1^{-1} + C_2^{-1})}{1 - \mathcal{C}_{3T}}, $$

which gives passing to the limit as $\epsilon \to 0$ the convergence $(\xi_\epsilon, \alpha_\epsilon) \to (\xi^*, \alpha^*)$. Since $\hbar$ is lower semi-continuous and using property (i) of Ekeland's principle, the inequality $\hbar(\xi_\epsilon, \alpha_\epsilon) \le \inf_{(\xi, \alpha) \in \mathcal{D}} \hbar(\xi, \alpha) + \epsilon$ implies (as $\epsilon \to 0$) that $\hbar(\xi^*, \alpha^*) \le \inf_{(\xi, \alpha) \in \mathcal{D}} \hbar(\xi, \alpha)$. Therefore, we have $\hbar(\xi^*, \alpha^*) = \inf_{(\xi, \alpha) \in \mathcal{D}} \hbar(\xi, \alpha)$. Hence, it suffices to take $K_T = \mathcal{C}_{3T}$ and we obtain the indicated result. $\square$

# 4 Numerical simulations

In this section, we run numerical simulations to illustrate the effect of optimal control strategies on disease transmission dynamics. Generally, it is not possible to solve optimal control problems analytically. Therefore, we use a numerical method to approximate the optimal solutions and display the results. Among the practical approaches of optimization algorithms, we use the Forward-Backward sweep method [24, 4, 15, 33] to approximate numerically the solution of problem (OCP). The model (3.1) and the adjoint
<figure><img src="image_3.png" /><figcaption>Figure 1: Simulation of system (3.1) with the optimal control variables for the same control costs $C_1 = C_2 = 50$.</figcaption></figure>

system (3.10) are approximated by a semi-implicit finite-difference numerical scheme [4]. However, these models can be approximated by much more sophisticated numerical schemes, such as the finite element method [29, 38]. Our numerical simulations are performed for the control time $T = 50$ and with the maximum age of infection $\tau_{max} = 30$. Baseline parameters $\Lambda = 700$, $\mu = 0.04$, $m = 0.25$ and $\delta(\tau) = 0.05$. Moreover, we use the age-dependent infection rate $\beta(\tau) = 0.00005(1 + \tau/(1 + 5\tau))$. The initial and boundary conditions are $S(0) = 10^4$, $V(0) = 0$ and $i_0(\tau) = 60(\tau + 20)e^{-0.4(\tau+50)}$. Let $I(t) = \int_0^{\tau_{max}} i(t, \tau)d\tau$ be the total number of infected individuals, and we further set $\xi_+ = 0.75$ and $\alpha_+ = 0.65$.

The numerical simulations of the optimization problem (OCP) are illustrated in Figure 1–2 for various weight factors $C_1$ and $C_2$. We observe that the application of optimal controls keeps susceptible individuals at a relatively low level and increases the population of vaccinated individuals for almost the entire intervention period. In addition, it is optimal to start vaccination at the maximum rate. However, when the cost of treatment is higher than that of vaccination (see figure 2), more effort should be devoted to vaccination than when control costs are identical (see figure 1). The application of the controls decreases with time and cancels out at the final instant in both figures, which is consistent with the theoretical result of the remark 3.1. As a result, the total population of infected people, which had remained relatively low, increased at the end of the intervention campaign, but not enough. This increase in the number of infected people is because controls are decreasing. Optimal treatment differs according to the age of infection. It is higher in people who have been infected for a short time, and decreases over time. This means that to effectively fight this disease, we need to treat people from the beginning of their infection.

## 5 Conclusion

In this study, we have formulated and analyzed an SVI age-structured model describing the dynamics of an infectious disease in a population that takes into account imperfect vaccine and treatment as control measures. We assumed that the vaccination of susceptible individuals does not confer permanent immunity to recipients, and the recovery of infected individuals confers temporary immunity on the
<figure><img src="image_4.png" /><figcaption>Figure 2: Simulation of system (3.1) with the optimal control variables for the control costs $C_1 = 50$ and $C_2 = 200$.</figcaption></figure>

recovered individuals. We provided the theoretical results of the model. We studied the well-posedness of solutions, which includes the nonnegativity, boundedness, and asymptotic smoothness of solutions. We highlighted the expression for the basic reproduction number $R_0$ and carried out the stability of equilibria of the model. We showed that when $R_0 < 1$, the model exhibits the phenomenon of backward bifurcation, where the stable equilibrium point coexists with a stable endemic equilibrium point. As a consequence, the condition $R_0 < 1$ is only a necessary but not sufficient condition for the disease to disappear. We shown that this bifurcation is caused by the rate of treatment of infected people. Assuming that vaccination and treatment rates are time-dependent functions, we have formulated and solved an optimal control problem aiming to minimize the number of infected individuals and the costs associated with vaccination and treatment. This problem is solved numerically using the Backward–Forward sweep method. Numerical results illustrate the fact that optimal vaccination and treatment can significantly reduce disease prevalence and slow the spread of infection in the population, while keeping the number of infected individuals relatively low. In addition, optimal control strategies remain dependent on control costs and the age of infection.

## Acknowledgments

The authors would like to thank the anonymous reviewers for their constructive comments.

## A Appendix A

We present in this appendix the recent result established by Martcheva and Inaba [30] to detect the presence of backward and forward bifurcations and the driving necessary and sufficient conditions for their occurrence in the infinite dynamical system.

**Theorem A.1** (Martcheva and Inaba [30]). Let $Y$ and $Z$ be Banach spaces, $x \in Y$ and $q \in \mathbb{R}$ is a parameter. We consider the following abstract differential equation

$$ \frac{dx}{dt} = F(x, q), \quad F: Y \times \mathbb{R} \to Z. \qquad (A.1) $$
Without loss of generality, we assume that $0$ is an equilibrium point of the system (that is $F(0, q) = 0$ for all $q \in \mathbb{R}$) and assume

1. $A := D_x F(0, q_0)$ is the linearization around the equilibrium 0 evaluated at a critical value of parameter $q_0$, such that $A$ is a closed operator with a simple isolated eigenvalue zero and remaining eigenvalues having negative real part. Let $\hat{v}_0$ be the unique (up to a constant) positive solution of $Av = 0$.

2. $F(\tau, q) \in C^2(U_0 \times I_0, Z)$ for some neighbourhood $U_0$ of 0 and interval $I_0$ containing $q_0$.

3. Assume $Z^*$ is the dual of $Z$ and $\langle \cdot, \cdot \rangle$ is the pairing between $Z$ and $Z^*$. Assume $\hat{v}_0^* \in Z^*$ is the unique (up to a constant) positive vector satisfying $\langle Ax, \hat{v}_0^* \rangle = 0$ for all $x \in Y$, that is $\dim(\ker A^*) = 1$ where $A^*$ is the adjoint of $A$ and $\ker A^* = \operatorname{span}\{\hat{v}_0^*\}$.

4. Assume $\langle D_{xq}^2 F(0, q_0)\hat{v}_0, \hat{v}_0^* \rangle \neq 0$ where $D_{xq}^2 F(0, q_0)$ is the second derivative of $F$ with respect to $x$ and $q$.

Then, the direction of the bifurcation is determined by the numbers $a = \langle D_{xx}^2 F(0, q_0)[\hat{v}_0, \hat{v}_0], \hat{v}_0^* \rangle$ and $b = \langle D_{xq}^2 F(0, q_0)\hat{v}_0, \hat{v}_0^* \rangle$, where $D_{xx}^2 F(0, q_0)[h_1, h_2]$ is the second derivative of $F$ concerning $x$ applied to the functions $h_1$ and $h_2$. If $b > 0$, then the bifurcation is backward if and only if $a > 0$ and forward if and only if $a < 0$.

**Remark A.1.** The components of the eigenvector $\hat{v}_0$, which corresponds to positive entries in the disease-free steady state, may be negative. This is in general the case with the partial differential equations as well [30].

B Proof of Lemma 3.1

1. (i) Let us consider the functions defined by

$$h_r^{(1)}(t, i, \xi) = e^{-\int_r^t [\mathcal{T}(\beta i(\sigma, .)) + \xi(\sigma) + \mu] d\sigma} \quad \text{and} \quad h_r^{(2)}(t, i) = e^{-\int_r^t [\mathcal{T}(\beta i(\sigma, .)) + \mu] d\sigma} \quad (B.1)$$

for $r > 0$ and $t \in [0, T]$. Then by direct computation based on the following arguments: for all $x, y \in \mathbb{R}_+$, we have $|e^{-x} - e^{-y}| \le |x - y|$, there exists the positive constants $K_k$ with $k \in \{1, 2, 3\}$ such that the functions $h_r^{(1)}$ and $h_r^{(2)}$ satisfy for $(t, i^k, \xi^k)$ with $k \in \{1, 2\}$ and $t \in [0, T]$, the following inequalities,

$$\begin{aligned} |h_r^{(1)}(t, i^1, \xi^1) - h_r^{(1)}(t, i^2, \xi^2)| &\le K_1 \int_0^t |\mathcal{T}(i^1(\sigma, .)) - \mathcal{T}(i^2(\sigma, .))| d\sigma + K_2 \int_0^t |\xi^1(\sigma) - \xi^2(\sigma)| d\sigma, \\ |h_r^{(1)}(t, i^1) - h_r^{(1)}(t, i^2)| &\le K_3 \int_0^t |\mathcal{T}(i^1(\sigma, .)) - \mathcal{T}(i^2(\sigma, .))| d\sigma. \end{aligned}$$

Using the method of integrating factors on the ordinary differential equations and the method of characteristics on the first, second, and third equations of the system (3.1) respectively, we obtain the following expression of state variables:

$$\begin{aligned} S(t) &= S_0 h_0^{(1)}(t, i, \xi) + \int_0^t [\Lambda + \mathcal{T}(\alpha(r, .)i(r, .))] h_r^{(1)}(t, i, \xi) dr; \\ V(t) &= V_0 h_0^{(2)}(t, i) + \int_0^t \xi(r) S(r) h_r^{(2)}(t, i) dr; \end{aligned} \quad (B.2)$$

$$i(t, \tau) = i_0(\tau - t) \frac{\pi_{\alpha}(t, \tau)}{\pi_{\alpha}(t, \tau - t)} \mathbb{1}_{\{\tau > t\}} + i(t - \tau, 0) \pi_{\alpha}(t, \tau) \mathbb{1}_{\{\tau \le t\}}.$$
From the above representation of solutions of system (3.1) and using the Lipschitz properties of the functions $h_r^{(1)}$ and $h_r^{(2)}$ for all $r \ge 0$ for the variables, the proof of these estimates follow similarly as in corresponding result found in [34, 33].

(ii) Follows like in Item (i).

2. Let $(\xi_n, \alpha_n) \to (\xi, \alpha)$ as $n \to \infty$, $(S_n, V_n, i_n)$ and $(S, V, i)$ be the state of system (3.1) corresponding to $(\xi_n, \alpha_n)$ and $(\xi, \alpha)$, respectively. By Riesz theorem, there is a subsequence still denoted $(\xi_n, \alpha_n)$, such that $(\xi_n^2, \alpha_n^2) \to (\xi^2, \alpha^2)$ almost everywhere in $[0, T] \times Q$ as $n \to \infty$. Thus, Lebesgue's dominated convergence theorem yields that $\|\xi_n\|_{L^1(0,T)}^2 \to \|\xi\|_{L^1(0,T)}^2$ and $\|\alpha_n\|_{L^1(Q)}^2 \to \|\alpha\|_{L^1(Q)}^2$ as $n \to \infty$. On the other hand, it follows that $\int_Q \chi(\tau)i_n(t, \tau)dtd\tau \to \int_Q \chi(\tau)i(t, \tau)dtd\tau$ as $n \to \infty$. By Fatou's Lemma, it follows that $\underline{\underline{\hbar}}(\xi, \alpha) \le \liminf_{n \to \infty} \underline{\underline{\hbar}}(\xi_n, \alpha_n)$. Hence, the functional $\underline{\underline{\hbar}}(\cdot, \cdot)$ is lower semi-continuous. This achieves the proof $\square$

## References

[1] B. Ainseba, Z. Feng, M. Iannelli, and F. Milner. Control strategies for tb epidemics. *SIAM. J. Appl. Math*, 77(1):82–107, 2017.

[2] B. Ainseba and M. Iannelli. Optimal screening in structured sir epidemics. *Math. Model. Nat. Phenom*, 7(3):12–27, 2012.

[3] S. Anita. *Analysis and control of Age-Dependent population dynamics*. Kluwer, Boston, MA, 2000.

[4] S.ANIȚA, V. Arnăutu, and V. Capasso. *An Introduction to Optimal Control Problems in Life Sciences and Economics*. Springer Science, New York, 2011.

[5] V. Barbu. *Optimal Control of Variational Inequalities*. Research Notes in Mathematics 100. Pitman Publishing, 1984.

[6] V. Barbu. *Mathematical methods in optimization of differential systems*, volume 310. Springer science & Business Media, 1994.

[7] F. Brauer. Backward bifurcations in simple vaccination models. *J. Math. Anal. Appl*, 298(2):418–431, 2004.

[8] G. Chowell and F. Brauer. *The Basic Reproduction Number of Infectious Diseases: Computation and Estimation Using Compartmental Epidemic Models*, pages 1–30. Springer Netherlands, Dordrecht, 2009.

[9] R. D. Demasse, J.-J. Tewa, S. Bowong, and Y. Emvudu. Optimal control for an age-structured model for the transmission of hepatitis B. *J. Math. Biol.*, (73):305–333, 2016.

[10] O. Diekmann, J. Heesterbeek, and J. Metz. On the definition and the computation of the basic reproduction ratio $R_0$ in models for infectious-diseases in heterogeneous populations. *J. Math. Biol.*, 28:365–82, 02 1990.

[11] A. Ducrot, Z. Liu, and P. Magal. Essential growth rate for bounded linear perturbation of non-densely defined cauchy problems. *J Math Anal Appl*, 341(1):501–518, 2008.

[12] I. Ekeland. On the variational principle. *J. Math. Anal. Appl.*, 47(2):324–353, 1974.

[13] K. Fister and S. Lenhart. Optimal control of a competitive system with age-structure. *J. Math. Anal. Appl.*, 291(2):526–537, 2004.
[14] K. R. Fister and S. Lenhart. Optimal harvesting in an age-structured predator–prey model. *Appl. Math. Optim.*, 54(1):1–15, 2006.

[15] Y. Fotso Fotso, S. Touzeau, F. Grognard, B. Tsanou, and S. Bowong. Optimal control of coffee berry borers: Synergy between bio-insecticide and traps. *J. Optim. Theory. Appl.*, 196:882–899, 2023.

[16] Y. Fotso Fotso, S. Touzeau, B. Tsanou, F. Grognard, and S. Bowong. Mathematical modelling of a pest in an age-structured crop model: The coffee berry borer case. *Appl. Math. Model.*, 110:193–206, 2022.

[17] Y. Fotso Fotso, S. Touzeau, B. Tsanou, F. Grognard, and S. Bowong. Optimal impulsive control of coffee berry borers in a berry age-structured epidemiological model. *ARIMA*, 38, 2024.

[18] A. Gumel. Causes of backward bifurcations in some epidemiological models. *J. Math. Anal. Appl*, 395(1):355–365, 2012.

[19] Z.-K. Guo, H. Xiang, and H.-F. Huo. Analysis of an age-structured tuberculosis model with treatment and relapse. *J. Math. Biol*, 82:1–37, 2021.

[20] H. W. Hethcote. Optimal ages of vaccination for measles. *Math. Biosci.*, 89(1):29 – 52, 1988.

[21] H. Inaba. *Age-structured population dynamics in demography and epidemiology*. Springer, Singapore, 2017.

[22] H.-D. Kwon, J. Lee, and S.-D. Yang. Optimal control of an age-structured model of HIV infection. *Appl. Math. Comput*, 219(5):2766–2779, 2012.

[23] J. P. La Salle. *The Stability of Dynamical Systems*. Society for Industrial and Applied Mathematics, 1976.

[24] S. Lenhart and J. T. Workman. *Optimal control applied to biological models*. Taylor & Francis Group, London SW15 2NU, UK, 2007.

[25] C.-L. Li, C.-H. Li, and C.-Y. Cheng. Analysis of an epidemiological model with age of infection, vaccination, quarantine and asymptomatic transmission. *J. Frankl. Inst.*, 2022.

[26] H. Liu, J. Yu, and G. Zhu. Global behaviour of an age-infection-structured HIV model with impulsive drug-treatment strategy. *J. Theor. Biol*, 253(4):749–754, 2008.

[27] Z. Luo, Z.-R. He, and W.-T. Li. Optimal birth control for predator–prey system of three species with age-structure. *Appl. Math. Comput.*, 155(3):665 – 685, 2004.

[28] P. Magal and S. Ruan. *Theory and Applications of Abstract Semilinear Cauchy Problems*, volume 201. Applied Mathematical Sciences Springer, Cham, 2018.

[29] M. Marin, A. Hobiny, and I. Abbas. The effects of fractional time derivatives in porothermoelastic materials using finite element method. *Mathematics*, 9(14):1606, 2021.

[30] M. Martcheva and H. Inaba. A Lyapunov–Schmidt method for detecting backward bifurcation in age-structured population models. *J. Biol. Dyn.*, 14(1):543–565, 2020.

[31] M. Martcheva and X.-Z. Li. Competitive exclusion in an infection-age structured model with environmental transmission. *J. Math. Anal. Appl.*, 408(1):225 – 246, 2013.
[32] M. Martcheva and H. Thieme. Progression age enhanced backward bifurcation in an epidemic model with super-infection. *J. Math. Biol.*, 46:385-424, 06 2003.

[33] E. Numfor. Optimal treatment in a multi-strain within-host model of HIV with age structure. *J. Math. Anal. Appl.*, (480), 2019.

[34] E. Numfor, S. Bhattacharya, S. Lenhart, and M. Martcheva. Optimal control in coupled within-host and between-host models. *Math. Model. Nat. Phenom.*, 9(4):171-203, 2014.

[35] E. J. Park, M. Iannelli, M. Y. Kim, and S. Anita. Optimal harvesting for periodic age-dependent population dynamics. *SIAM J. Appl. Math*, (58):1648-1666, 1999.

[36] A. Perasso. An introduction to the basic reproduction number in mathematical epidemiology. *ESAIM: ProcS*, 62:123-138, 2018.

[37] Z. Qiu and Z. Feng. Transmission dynamics of an influenza model with age of infection and antiviral treatment. *J. Dyn. Differ. Equ.*, 22(4):823-851, 2010.

[38] T. Saeed and I. Abbas. Finite element analyses of nonlinear dpl bioheat model in spherical tissues using experimental data. *Mech. Based Des. Struct. Mach.*, 50(4):1287-1297, 2022.

[39] H. Thieme. Semiflows generated by lipschitz perturbations of non-densely defined operators. *Differ. Integral Equ.*, 3(6):1035-1066, 1990.

[40] J. Wang, R. Zhang, and T. Kuniya. The dynamics of an SVIR epidemiological model with infection age. *IMA J. Appl. Math.*, 81(2):321-343, 12 2015.

[41] D. Yan and H. Cao. The global dynamics for an age-structured tuberculosis transmission model with the exponential progression rate. *Appl. Math. Model*, 75:769-786, 2019.

[42] J. Yang, M. Zhou, and X. Li. Backward bifurcation of an age-structured epidemic model with partial immunity: The lyapunov-schmidt approach. *Appl. Math. Lett.*, 133:108292, 2022.