Global stability of a SEIR discrete delay
differential-difference system with protection phase

Mostafa Adimy, Abdennasser Chekroun, Toshikazu Kuniya, Hanene Meghelli

► To cite this version:

Mostafa Adimy, Abdennasser Chekroun, Toshikazu Kuniya, Hanene Meghelli. Global stability of a SEIR discrete delay differential-difference system with protection phase. Mathematical Methods in the Applied Sciences, 2023, 2023, pp.21. 10.1002/mma.9533 . hal-04273580

HAL Id: hal-04273580

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

Submitted on 7 Nov 2023

**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 CC BY 4.0 - Attribution - International License
Global stability of a SEIR discrete delay
differential-difference system with protection
phase

Mostafa Adimy$^a$, Abdennasser Chekroun$^{b*}$, Toshikazu Kuniya$^c$
and Hanene Meghelli$^b$

June 23, 2023

Abstract

We consider an epidemiological model with the four classical compart-
ments of susceptible, exposed, infected and recovered population. We add
a new compartment that is supposed to describe, for a limited time, indi-
viduals that are protected from the epidemic through vaccination or med-
ication, for instance. We model the protection phase by an age-structured
partial differential equation. The age is the time since an individual en-
tered the protection phase. The model is then reduced by integration on
the characteristics to a differential-difference system with delay. The dis-
crete delay represents the limited duration of the protection phase. After
establishing the basic properties of the model, we show that the disease-
free equilibrium (DFE) is globally asymptotically stable when the basic
reproduction number is less than one and is unstable when this number is
greater than one. Furthermore, we show that even if there is no mortality
during the protection phase and the basic reproduction number is greater
than one, the endemic equilibrium is globally asymptotically stable. The
proofs of the global asymptotic stability of both equilibria are based on
carefully constructed Lyapunov functions. To complete this study on the
global dynamics, we discuss some results on weak and strong uniform per-
sistence of the disease. Finally, numerical simulations are performed to
illustrate and complete our main results.

Keywords: SIR epidemic model, delay differential-difference system, global asymptotic stability, Lyapunov functional.

MSC: 34D23, 37N25, 92D30, 34K06

<footnote>a Univ Lyon, Inria, Université Claude Bernard Lyon 1, CNRS UMR5208,
Institut Camille Jordan, F-69603 Villeurbanne, France.</footnote>

<footnote>b Laboratoire d'Analyse Nonlinéaire et Mathématiques Appliquées, University of Tlemcen, Tlemcen 13000, Algeria.</footnote>

<footnote>c Graduate School of System Informatics, Kobe University, 1-1 Rokkodai-cho, Nada-ku, Kobe 657-8501, Japan.</footnote>

<footnote>*corresponding author : chekroun@math.univ-lyon1.fr</footnote>
# 1 Introduction

Ordinary (ODE) and partial differential equations (PDE) have long played an important role in understanding, modeling and forecasting the spread of epidemics, [17, 26, 31, 37]. Many studies have been conducted on the dynamics of the SIR (S: susceptible, I: infected and R: recovered) and SEIR (E: exposed) epidemic models (see, for example, [8, 13, 12, 29, 32] and the references cited therein).

Sometimes simple ODE and PDE models cannot capture the rich variety of dynamics observed during an epidemic. In many cases, differential equations or partial differential equations with or without time delays are used, [25, 28]. The interpretation of time delays depends on the problem. It may model the incubation period or exposure to the disease, the duration of infection, the duration of immunity after infection or vaccination, etc (see, for instance, [11, 12, 15, 19, 20, 23, 30, 31, 33, 37, 39]). In most of these works, necessary and sometimes sufficient conditions for the global stability of disease-free and endemic steady-states have been addressed and in some cases the existence of possible periodic solutions has been obtained.

Recently in [1, 2, 3, 4, 5, 6, 7, 8], we have proposed and studied a new class of models based on a combination of differential or partial differential and difference equations with discrete and distributed delay. In [1, 3, 7, 8], we studied systems with and without spatial structure applied to blood cell regulation and production in the bone marrow (hematopoiesis). We considered systems of differential-difference equations with delay where the delay represents the duration of the cell cycle. We studied the existence and stability of steady-states, the existence of a Hopf bifurcation for the positive equilibrium, the existence and non-existence of monotone traveling wave front solutions in the case of an unbounded space, etc. In [2, 5], we applied these techniques to a generalization of the Kermack-McKendrick SIR model where we also combined a differential system (or PDE) with a difference equation with time delay. We focused on epidemic models of susceptible, infected and recovered (SIR) individuals considering a time-limited phase of protection, e.g. by vaccination or drugs. In these models, the delay represents the duration of the protection phase. We studied the global asymptotic stability of the two steady-states: disease-free and endemic. We designed quadratic and logarithmic Lyapunov functions to establish this global stability. We proved that the global stability is completely determined by the basic reproduction number. In [6], we studied the existence and non-existence of non-trivial traveling wave solutions for a general class of diffusive Kermack-McKendrick SIR epidemic models. We obtained almost complete information on the threshold and minimum wave speed that describe the transition between the existence and non-existence of such non-trivial traveling waves that indicate whether the epidemic can propagate or not. In most of the situations we encountered, we made a rather extensive study of the models. However, a complete theory of coupled differential-difference systems remains to be done. In this sense, an important effort has already been brought in the Lyapunov-Krasovskii type functional approach for this type of systems (see [10, 24, 34, 35] and references therein).

For many infections, there is an important latency period during which individuals have been infected but are not yet infectious themselves. During this period, individuals are in a compartment called E (for exposed). In this pa-
per, we consider a more general class of Kermack-McKendrick models than the one studied in [2], by adding the E compartment. We then obtain a model of susceptible, exposed, infected and recovered (SEIR) individuals with a time-limited protection phase described by a deference equation. The main result of this work is the global asymptotic stability of the two steady-states: disease-free and endemic. The proof of the global stability is based on carefully constructed Lyapunov functions.

The organization of the paper is as follows. We present the mathematical model in Section 2. This model is based on a combination of differential and difference equations with delay. In this same section, we discuss results concerning the existence, uniqueness, positivity and uniform boundedness of the solutions. Section 3 is devoted to the basic reproduction number and the expression of the endemic steady-state as a function of this number. In Section 4, we prove the local and global asymptotic stability of the disease-free steady-state. For the global stability, we use a Lyapunov functional. In Section 5, we deal with the weak and strong uniform persistence of the disease. We also construct a Lyapunov functional to show the global attractiveness of the endemic steady-state. In Section 6, we illustrate and complete our results with a discussion and numerical simulations.

## 2 Presentation of the model

In this work, we study a SEIR model with a protection phase. The model is represented by a differential-difference system with delay. More precisely, the population is divided into five compartments of individuals: the susceptible $S(t)$, the infected $I(t)$, the exposed $E(t)$, the recovered $R(t)$ and the protected individuals $P(t)$ (see Figure 1). We structure the protected individuals according to the time $a$ since they entered this compartment, and we denote their density by $p(t, a)$, with $a \in (0, \tau)$ (see [2, 5] for more details). Then, we have, for $t > 0$,

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

The model is given, for $t > 0$, by

$$ 
\left\{
\begin{aligned}
S'(t) &= \Lambda - (h + a_S)S(t) - \beta S(t)I(t) + (1 - \alpha)e^{-a_p \tau}u(t - \tau), \\
E'(t) &= \beta S(t)I(t) - (\mu + a_E)E(t), \\
I'(t) &= \mu E(t) - (r + a_I)I(t), \\
R'(t) &= rI(t) - a_R R(t), \\
u(t) &= hS(t) + \alpha e^{-a_p \tau}u(t - \tau),
\end{aligned}
\right.
\qquad (1) $$

where $u(t) = p(t, 0)$, $t > 0$, represents the new protected individuals. We note

$$\phi(t) := e^{-a_p t} p(0, -t), \quad -\tau \le t \le 0.$$

Hence,

$$p(t, \tau) = e^{-a_p \tau} \begin{cases} u(t - \tau), & t > \tau, \\ \phi(t - \tau), & 0 \le t \le \tau. \end{cases}$$
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic illustration of the model.</figcaption></figure>

We obtain

$$P(t) = \int_{0}^{\tau} e^{-a_p a} p(t-a, 0) da = \int_{0}^{\tau} e^{-a_p a} u(t-a) da, \quad t > 0.$$

By derivating $P$, we obtain, for $t > 0$,

$$P'(t) = -a_p P(t) + u(t) - e^{-a_p \tau} u(t - \tau).$$

This system is completed by the following initiates conditions,

$$S(0) = S_0, E(0) = E_0, I(0) = I_0, R(0) = R_0 \text{ and } u(\theta) = \phi(\theta), \text{ for } -\tau \le \theta \le 0, \tag{2}$$

with $\phi$ a continuous function on the interval $[-\tau, 0]$. All parameters are positive and their interpretation is as follows. $\Lambda$ represents births and immigration, $h$ is the protection rate, $a_S$, $a_p$, $a_E$, $a_I$, and $a_R$ are the mortality rates, $\beta$ is the contact rate, $\mu$ is the infection rate, $r$ is the recovery rate, $\tau$ is the temporary protection period, and $\alpha \in (0, 1)$ is the fraction of individuals who, at the end of the previous period of protection, enter a new phase of protection, for example, through a vaccine booster or a renewal of treatment (see, [2, 5, 6]).

**Remark 1.** In the case $\alpha = 0$, which means that no individual re-enters the protection phase (no individual updates its vaccine). Then, the system becomes a classical delayed differential equation. Indeed, we obtain $u(t) = hS(t)$. Then, $S$ satisfies

$$S'(t) = \Lambda - (h + a_S)S(t) - \beta S(t)I(t) + h e^{-a_p \tau} S(t - \tau).$$
The total population $N = S + E + I + R + P$, satisfies, for $t > 0$,

$N'(t) \leq \Lambda - \rho N(t)$, where $\rho := \min\{a_S, a_E, a_I, a_R, a_p\} > 0$.

We obtain

$$ \limsup_{t \to +\infty} N(t) \leq \frac{\Lambda}{\rho}. $$

This means that each component is uniformly bounded.

We can show the existence, uniqueness and positivity of the solutions of System (1)-(2) by the method of steps. More precisely, we have the following theorem which is easy to prove and its proof is omitted.

**Theorem 2.1.** Assume that the initial condition (2) satisfies the compatibility condition

$$ \phi(0) = hS_0 + \alpha e^{-a_p\tau}\phi(-\tau), $$

with $\phi$ continuous on the interval $[-\tau, 0]$. Then, the system (1)-(2) admits a unique solution $(S, E, I, R, u)$, where $(S, E, I, R)$ has a continuous first derivative for all $t > 0$ and the function $u$ is continuous for all $t \geq -\tau$. Moreover, all solutions with nonnegative initial conditions are nonnegative.

*Proof.* The existence, uniqueness and regularity of the solution can be proved by steps in the intervals $[(k-1)\tau, k\tau]$, $k = 1, 2, \dots$. Now assume that the initial condition (2) is nonnegative and prove that the solution of the system (1) is nonnegative. It is clear that for all $t \in [0, \tau]$, $u(t - \tau) = \phi(t - \tau) \geq 0$. If $S(0) = 0$, then $S'(0) = \Lambda + (1 - \alpha)e^{-a_p\tau}\phi(-\tau) > 0$. In any case, there exists $t_0 \in (0, \tau)$ such that $S(t) > 0$ for $t \in (0, t_0)$. Suppose by contradiction that there exists $t_1 \in [t_0, \tau]$ such that $S(t_1) = 0$. We choose $t_1$ as the first point such that $S(t_1) = 0$. Then, we obtain $S'(t_1) = \Lambda + (1 - \alpha)e^{-a_p\tau}\phi(t_1 - \tau) > 0$ and this is clearly absurd. We proved that $S(t)$ is nonnegative on the interval $[0, \tau]$. By steps, we can extend this reasoning to the whole interval $[0, +\infty)$. The functions $S(t)$ and $u(t)$ being positive, we can consider the system in $(E, I, R)$ as an ordinary linear differential equation with time varying coefficients. Therefore, if the initial condition is nonnegative, the solution is also nonnegative. $\square$

# 3 Steady-states and basic reproduction number

In this section, we establish the existence of steady-states and derive the basic reproduction number associated with the system (1). An equilibrium $(\bar{S}, \bar{E}, \bar{I}, \bar{R}, \bar{u})$ of the system (1) satisfies

$$ \left\{ \begin{array}{l} \Lambda - (h + a_S)\bar{S} - \beta\bar{S}\bar{I} + (1 - \alpha)e^{-a_p\tau}\bar{u} = 0, \\ \beta\bar{S}\bar{I} - (\mu + a_E)\bar{E} = 0, \\ \mu\bar{E} - r\bar{I} - a_I\bar{I} = 0, \\ r\bar{I} - a_R\bar{R} = 0, \\ \bar{u} = h\bar{S} + \alpha e^{-a_p\tau}\bar{u}. \end{array} \right. \qquad (3) $$

From the third and fourth equations of (3), we have

$$ \bar{E} = \frac{r + a_I}{\mu} \bar{I} \quad \text{and} \quad \bar{R} = \frac{r}{a_R} \bar{I}. \qquad (4) $$
Moreover, the last equation of (3) gives

$$ \bar{u} = \frac{h}{1 - \alpha e^{-a_p \tau}} \bar{S}. \qquad (5) $$

By replacing (4) in the second equation of (3), we obtain

$$ \beta \bar{S} \bar{I} - (\mu + a_E)(r + a_I) \frac{\bar{I}}{\mu} = 0. $$

This implies that

$$ \bar{I} = 0 \quad \text{or} \quad \bar{S} = \frac{(\mu + a_E)(r + a_I)}{\beta \mu}. $$

If $\bar{I} = 0$, then $\bar{E} = \bar{R} = 0$ and the first equation of (3) gives

$$ \Lambda - (h + a_S)\bar{S} + (1 - \alpha)e^{-a_p \tau} \frac{h}{1 - \alpha e^{-a_p \tau}} \bar{S} = 0. $$

Then, we obtain

$$ \bar{S} = S^0 := \frac{\Lambda(1 - \alpha e^{-a_p \tau})}{(h + a_S)(1 - \alpha e^{-a_p \tau}) - h(1 - \alpha)e^{-a_p \tau}} = \frac{\Lambda}{a_S + h \frac{1 - e^{-a_p \tau}}{1 - \alpha e^{-a_p \tau}}}, \quad (6) $$

and

$$ u^0 = \frac{h}{1 - \alpha e^{-a_p \tau}} S^0. $$

This means that $(S^0, 0, 0, 0, u^0)$ is always an equilibrium of the system (1). This is the disease-free equilibrium of the system (1).

Now, we assume that $\bar{I} > 0$. Then,

$$ \bar{S} = \frac{(\mu + a_E)(r + a_I)}{\beta \mu} \quad \text{and} \quad \bar{I} = \frac{\Lambda - (h + a_S)\bar{S} + (1 - \alpha)e^{-a_p \tau} \frac{h}{1 - \alpha e^{-a_p \tau}} \bar{S}}{\beta \bar{S}}. \qquad (7) $$

$\bar{I}$ exists if and only if

$$ \Lambda - \left[ \frac{(h + a_S)(1 - \alpha e^{-a_p \tau}) - (1 - \alpha)e^{-a_p \tau}h}{1 - \alpha e^{-a_p \tau}} \right] \bar{S} > 0. $$

The above inequality is equivalent to

$$ \frac{\Lambda(1 - \alpha e^{-a_p \tau})\beta\mu}{[a_S(1 - \alpha e^{-a_p \tau}) + h(1 - e^{-a_p \tau})](\mu + a_E)(r + a_I)} > 1. $$

Therefore, we obtain a threshold $\mathcal{R}_0 > 1$ for the existence of a unique endemic equilibrium given by $(\bar{S}, \bar{E}, \bar{I}, \bar{R}, \bar{u}) = (S^*, E^*, I^*, R^*, u^*)$, satisfying (4), (5), and (7). This threshold is given by

$$ \begin{aligned} \mathcal{R}_0 &:= \frac{\Lambda(1 - \alpha e^{-a_p \tau})\beta\mu}{[a_S(1 - \alpha e^{-a_p \tau}) + h(1 - e^{-a_p \tau})](\mu + a_E)(r + a_I)}, \\ &= \frac{\Lambda\beta\mu}{[a_S + h \frac{1 - e^{-a_p \tau}}{1 - \alpha e^{-a_p \tau}}](\mu + a_E)(r + a_I)}, \qquad (8) \\ &= \frac{\beta\mu}{(\mu + a_E)(r + a_I)} S^0, \end{aligned} $$
where $S^0$ is given by (6). We defined the basic reproduction number $\mathcal{R}_0$ associate to the differential-difference system (1) as the threshold $\mathcal{R}_0 > 1$ that allows the existence of the endemic equilibrium. In fact, we have simply proved the following theorem.

**Theorem 3.1.** 1. The disease-free equilibrium $(S^0, 0, 0, 0, u^0)$ of System (1) always exists and it is defined by

$$S^0 = \frac{\Lambda}{a_S + h \frac{1-e^{-a_p \tau}}{1-\alpha e^{-a_p \tau}}} \quad \text{and} \quad u^0 = \frac{h}{1-\alpha e^{-a_p \tau}} S^0.$$

2. The endemic equilibrium of System (1) exists if and only if $\mathcal{R}_0 > 1$, where $\mathcal{R}_0$ is given by (8), and this endemic equilibrium is defined by

$$S^* = \frac{(\mu + a_E)(r + a_I)}{\beta \mu},$$

$$E^* = \frac{r + a_I}{\mu} I^*,$$

$$I^* = \frac{\Lambda - (h + a_S)S^* + h \frac{(1-\alpha)e^{-a_p \tau}}{1-\alpha e^{-a_p \tau}} S^*}{\beta S^*} = \frac{1}{\beta} \left[ \frac{\Lambda}{S^*} - \left( a_S + h \frac{1 - e^{-a_p \tau}}{1 - \alpha e^{-a_p \tau}} \right) \right],$$

$$R^* = \frac{r}{a_R} I^*,$$

$$u^* = \frac{h}{1 - \alpha e^{-a_p \tau}} S^*.$$

In the rest of this paper, we will prove the well-known property associated with $\mathcal{R}_0$ as a threshold, which says that epidemic can invade if $\mathcal{R}_0 > 1$, whereas it cannot if $\mathcal{R}_0 < 1$.

# 4 Local and global asymptotic stability of the disease-free equilibrium

## 4.1 Local asymptotic stability of the disease-free equilibrium

We deal with the local asymptotic stability of $(S^0, 0, 0, 0, u^0)$. In the case $\mathcal{R}_0 < 1$, $(S^0, 0, 0, 0, u^0)$ is the only equilibrium. We notice that in the system (1), $S, E, I$ and $u$ are independent of $R$. So, in the following we can neglect the component $R$, and consider the system

$$\left\{
\begin{aligned}
S'(t) &= \Lambda - (h + a_S)S(t) - \beta S(t)I(t) + (1 - \alpha)e^{-a_p \tau}u(t - \tau), \\
E'(t) &= \beta S(t)I(t) - (\mu + a_E)E(t), \\
I'(t) &= \mu E(t) - (r + a_I)I(t), \\
u(t) &= hS(t) + \alpha e^{-a_p \tau}u(t - \tau).
\end{aligned}
\right.
\qquad (9)$$

As in [7], we linearize the differential-difference system (9) and derive the characteristic equation. The linearized system of (9) about any equilibrium $(\bar{S}, \bar{E}, \bar{I}, \bar{u})$
is given by

$$
\begin{equation}
\begin{cases}
S'(t) &= -(h + a_S)S(t) - \beta \bar{S} I(t) - \beta \bar{I} S(t) + (1 - \alpha)e^{-a_p \tau} u(t - \tau), \\
E'(t) &= \beta \bar{S} I(t) + \beta \bar{I} S(t) - (\mu + a_E)E(t), \\
I'(t) &= \mu E(t) - (r + a_I)I(t), \\
u(t) &= hS(t) + \alpha e^{-a_p \tau} u(t - \tau).
\end{cases}
\end{equation}
$$

We consider the functions S, E, I and u, given by

$$
S(t) = c_1 e^{\lambda t}, \quad E(t) = c_2 e^{\lambda t}, \quad I(t) = c_3 e^{\lambda t} \text{ and } u(t) = c_4 e^{\lambda t}, \quad \lambda \in \mathbb{C}, \quad c_i \in \mathbb{R}_+^*.
$$

Then, we obtain the following system

$$
\begin{pmatrix}
\lambda + h + a_S + \beta \bar{I} & 0 & \beta \bar{S} & -(1-\alpha)e^{-a_p \tau} e^{-\lambda \tau} \\
-\beta \bar{I} & \lambda + \mu + a_E & -\beta \bar{S} & 0 \\
0 & -\mu & \lambda + r + a_I & 0 \\
-h & 0 & 0 & 1 - \alpha e^{-a_p \tau} e^{-\lambda \tau}
\end{pmatrix}
\begin{pmatrix}
c_1 \\
c_2 \\
c_3 \\
c_4
\end{pmatrix}
=
\begin{pmatrix}
0 \\
0 \\
0 \\
0
\end{pmatrix}.
$$

The characteristic equation is given by

$$
\det A_{\lambda}^{\tau} = 0,
$$

with

$$
A_{\lambda}^{\tau} = \begin{pmatrix} \lambda + h + a_S + \beta \bar{I} & 0 & \beta \bar{S} & -(1-\alpha)e^{-a_p \tau} e^{-\lambda \tau} \\ -\beta \bar{I} & \lambda + \mu + a_E & -\beta \bar{S} & 0 \\ 0 & -\mu & \lambda + r + a_I & 0 \\ -h & 0 & 0 & 1 - \alpha e^{-a_p \tau} e^{-\lambda \tau} \end{pmatrix}.
$$

For the disease-free equilibrium ($S^0$, $0, 0, u^0$), the characteristic equation becomes

$$
\Delta_0(\lambda) = \left[-h(1-\alpha)e^{-a_p \tau - \lambda \tau} + (\lambda + h + a_S)(1-\alpha e^{-a_p \tau - \lambda \tau})\right] \\
\times \left[(\lambda + \mu + a_E)(\lambda + r + a_I) - \beta \mu S^0\right] = 0. \tag{10}
$$

We obtain the following theorem.

**Proposition 1.** (i) Assume that $\mathcal{R}_0 > 1$. Then, there exists a positive real root of (10) and the disease-free steady-state $(S^0, 0, 0, u^0)$ is unstable.

(ii) Assume that $\mathcal{R}_0 < 1$. Then, all the roots of the characteristic equation (10) have negative real parts, and the disease-free steady-state $(S^0, 0, 0, u^0)$ is locally asymptotically stable.

*Proof.* From the characteristic equation (10), we have either

$$
-h(1-\alpha)e^{-a_p \tau - \lambda \tau} + (\lambda + h + a_S)(1-\alpha e^{-a_p \tau - \lambda \tau}) = 0, \quad (11)
$$

or

$$
(\lambda + \mu + a_E)(\lambda + r + a_I) - \beta\mu S^0 = 0. \quad (12)
$$

We suppose that $\lambda = \xi + i\omega$, with $\xi \ge 0$, is a root of Equation (11). Then, we have

$$
e^{-(\lambda+a_p)\tau} = \frac{\lambda+h+a_S}{\alpha\lambda+h+\alpha a_S}.
$$
This means that

$$ \left| \frac{\lambda + h + a_S}{\alpha\lambda + h + \alpha a_S} \right|^2 < 1. $$

By developing the modulus, we obtain

$$ (1 - \alpha^2)\omega^2 + (1 - \alpha)(\xi + a_S)[(1 + \alpha)(\xi + a_S) + 2h] < 0. $$

Since $\alpha \in (0, 1)$, we obtain a contradiction. Then, no root with positive real part satisfies Equation (11). Therefore, the local stability of $(S^0, 0, 0, u^0)$ is determined by the sign of the real part of $\lambda \in \mathbb{C}$ satisfying (12). This equation is equivalent to

$$ \lambda^2 + (\mu + a_E + r + a_I)\lambda + (\mu + a_E)(r + a_I) - \beta\mu S^0 = 0. $$

We have $\mu + a_E + r + a_I > 0$ and

$$ (\mu + a_E)(r + a_I) - \beta\mu S^0 = (\mu + a_E)(r + a_I)(1 - \mathcal{R}_0). $$

Then, according to the criterion of Routh-Hurwitz, the disease-free equilibrium $(S^0, 0, 0, u^0)$ is locally asymptotically stable if $\mathcal{R}_0 < 1$ and unstable if $\mathcal{R}_0 > 1$. $\square$

## 4.2 Global asymptotic stability of the disease-free equilibrium

To prove the global asymptotic stability of the disease-free equilibrium $(S^0, 0, 0, 0, u^0)$, we construct a Lyapunov functional. This is stated in the following theorem.

**Theorem 4.1.** Assume that $\mathcal{R}_0 < 1$. Then, the disease-free equilibrium $(S^0, 0, 0, 0, u^0)$ of the system (1) is globally asymptotically stable.

*Proof.* Let $(S^0, 0, 0, 0, u^0)$ be the disease-free equilibrium. By putting, for $t \ge 0$,

$$ \tilde{S}(t) = S(t) - S^0 \quad \text{and} \quad \tilde{u}_t = u_t - u^0, $$

where the notation $u_t, t \ge 0$, means the function $u_t : \theta \in [-\tau, 0] \mapsto u_t(\theta) := u(t + \theta)$, the system (1) becomes

$$ \left\{ \begin{aligned} \tilde{S}'(t) &= -(h + a_S)\tilde{S}(t) - \beta\tilde{S}(t)I(t) - \beta S^0 I(t) + (1 - \alpha)e^{-a_p\tau}\tilde{u}(t - \tau), \\ E'(t) &= \beta\tilde{S}(t)I(t) + \beta S^0 I(t) - (\mu + a_E)E(t), \\ I'(t) &= \mu E(t) - (r + a_I)I(t), \\ R'(t) &= rI(t) - a_R R(t), \\ \tilde{u}(t) &= h\tilde{S}(t) + \alpha e^{-a_p\tau}\tilde{u}(t - \tau). \end{aligned} \right. \quad (13) $$

We construct the following Lyapunov functional, $t \ge 0$,

$$ \begin{gathered} V : \mathbb{R} \times \mathbb{R}_+^3 \times C([-\tau, 0], \mathbb{R}) \to \mathbb{R}^+ \\ (\tilde{S}(t), E(t), I(t), R(t), \tilde{u}_t) \mapsto \tilde{V}(t) := V(\tilde{S}(t), E(t), I(t), R(t), \tilde{u}_t), \end{gathered} $$
defined by

$$
\begin{align*}
\tilde{V}(t) = {}& 2 \frac{(h+a_S) - h\alpha e^{-2a_p\tau} - a_S\alpha^2 e^{-2a_p\tau}}{(e^{-a_p\tau}(1-\alpha))^2} \tilde{S}^2(t) + \int_{-\tau}^{0} \tilde{u}_t^2(s) ds \\
& + 2 \frac{(h+a_S) - h\alpha e^{-2a_p\tau} - a_S\alpha^2 e^{-2a_p\tau}}{(e^{-a_p\tau}(1-\alpha))^2} S^0 \left[ E(t) + \frac{\mu+a_E}{\mu} (I + \eta R(t)) \right],
\end{align*}
$$

where η is a constant satisfying

$$
0 < \eta < \frac{\mu\beta(1 - \mathcal{R}_0)}{r(\mu + a_E)}.
$$

We set

$$
\epsilon_1 := 2 \frac{(h + a_S) - h\alpha e^{-2a_p\tau} - a_S\alpha^2 e^{-2a_p\tau}}{(e^{-a_p\tau}(1-\alpha))^2} > 0. \quad (14)
$$

We compute the derivative of V along the solution trajectory. We get, for $t \ge 0$,

$$
\begin{align*}
\tilde{V}'(t) &= \epsilon_1 \tilde{S}(t) \tilde{S}'(t) + \tilde{u}^2(t) - \tilde{u}^2(t - \tau) \\
&\quad + \epsilon_1 S^0 \left( E'(t) + \frac{\mu + a_E}{\mu} I'(t) + \frac{\mu + a_E}{\mu} \eta R'(t) \right), \\
&= \epsilon_1 \tilde{S}(t) \left[ -(h + a_S) \tilde{S}(t) + (1-\alpha)e^{-a_p\tau} \tilde{u}(t-\tau) \right] + \tilde{u}^2(t) - \tilde{u}^2(t - \tau) \\
&\quad - \epsilon_1 \beta \tilde{S}^2(t) I(t) - \epsilon_1 \beta \tilde{S}(t) I(t) S^0 + \epsilon_1 \beta S^0 \tilde{S}(t) I(t) + \epsilon_1 S^0 \beta S^0 I(t) \\
&\quad - \epsilon_1 S^0 (\mu + a_E) E(t) + \epsilon_1 S^0 (\mu + a_E) E(t) - \epsilon_1 S^0 \frac{\mu + a_E}{\mu} (r + a_I) I(t) \\
&\quad + \epsilon_1 S^0 \frac{\mu + a_E}{\mu} r \eta I(t) - \epsilon_1 S^0 \frac{\mu + a_E}{\mu} \eta a_R R(t).
\end{align*}
$$

Then, we have

$$
\begin{align*}
\tilde{V}'(t) ={}& -\epsilon_1(h+a_S)\tilde{S}^2(t) + \epsilon_1(1-\alpha)e^{-a_p\tau}\tilde{u}(t-\tau)\tilde{S}(t) + h^2\tilde{S}^2(t) \\
& + (\alpha e^{-a_p\tau})^2\tilde{u}^2(t-\tau) + 2h\alpha e^{-a_p\tau}\tilde{S}(t)\tilde{u}(t-\tau) - \tilde{u}^2(t-\tau) \\
& -\epsilon_1\beta\tilde{S}^2(t)I(t) + \epsilon_1 S^0\beta S^0 I(t) - \epsilon_1 S^0 \frac{\mu+a_E}{\mu}(r+a_I)I(t) \\
& + \epsilon_1 S^0 \frac{\mu+a_E}{\mu}r\eta I(t) - \epsilon_1 S^0 \frac{(\mu+a_E)}{\mu}\eta a_R R(t).
\end{align*}
$$

Therefore, we get

$$
\begin{align*}
\tilde{V}'(t) ={}& [-\epsilon_1(h+a_S) + h^2] \tilde{S}^2(t) + [\epsilon_1 e^{-a_p\tau}(1-\alpha) + 2h\alpha e^{-a_p\tau}] \tilde{u}(t-\tau) \tilde{S}(t) \\
& + \tilde{u}^2(t-\tau) [ (\alpha e^{-a_p\tau})^2 - 1 ] \\
& - \epsilon_1 \beta \tilde{S}^2(t) I(t) + \epsilon_1 \beta S^0 I(t) \frac{(\mu+a_E)(r+a_I)}{\mu\beta} \left( \mathcal{R}_0 - 1 + \frac{\mu+a_E}{\mu\beta} r\eta \right) \\
& - \epsilon_1 S^0 \frac{(\mu+a_E)}{\mu} \eta a_R R(t).
\end{align*}
$$

We denote

$$
\begin{align*}
A &:= \epsilon_1 (h + a_S) - h^2, \\
B &:= \epsilon_1 e^{-a_p \tau} (1 - \alpha) + 2h \alpha e^{-a_p \tau}, \\
C &:= -(\alpha e^{-a_p \tau})^2 + 1 > 0.
\end{align*}
$$
We claim that $A > 0$. In fact,

$$A = 2 \frac{(h + a_S)^2 - \alpha e^{-2a_p \tau} (h + \alpha a_S)(h + a_S) - h^2 (e^{-a_p \tau} (1 - \alpha))^2}{(e^{-a_p \tau} (1 - \alpha))^2}.$$

Then, we claim that $(h+a_S)^2 - \alpha e^{-2a_p \tau} (h+\alpha a_S)(h+a_S) - h^2 (e^{-a_p \tau} (1-\alpha))^2 > 0$.
Indeed,

$$ \begin{aligned} & (h + a_S)^2 - \alpha e^{-2a_p \tau} (h + \alpha a_S)(h + a_S) - h^2 (e^{-a_p \tau} (1 - \alpha))^2 \\ &= h^2 + a_S^2 + 2ha_S - \alpha e^{-2a_p \tau} (h^2 + ha_S + \alpha ha_S + \alpha a_S^2) - h^2 (e^{-a_p \tau} (1 - \alpha))^2, \\ &= a_S^2 (1 - \alpha^2 e^{-2a_p \tau}) + 2ha_S - ha_S \alpha e^{-2a_p \tau} - \alpha^2 ha_S e^{-2a_p \tau} \\ &\quad + h^2 (1 - \alpha e^{-2a_p \tau} - \alpha^2 e^{-2a_p \tau} + 2\alpha e^{-2a_p \tau}), \\ &= a_S^2 (1 - \alpha^2 e^{-2a_p \tau}) + ha_S (2 - \alpha e^{-2a_p \tau} - \alpha^2 e^{-2a_p \tau}) \\ &\quad + h^2 (1 - e^{-2a_p \tau} - \alpha^2 e^{-2a_p \tau} + \alpha e^{-2a_p \tau}). \end{aligned} $$

Since $0 < \alpha < 1$, then

$$a_S^2(1 - \alpha^2 e^{-2a_p \tau}) + ha_S(2 - \alpha e^{-2a_p \tau} - \alpha^2 e^{-2a_p \tau}) + h^2(1 - e^{-2a_p \tau} - \alpha^2 e^{-2a_p \tau} + \alpha e^{-2a_p \tau}) > 0.$$

Using the notations above, we obtain

$$ \begin{aligned} \tilde{V}'(t) &= -A\tilde{S}^2(t) + B\tilde{S}(t)\tilde{u}(t-\tau) - C\tilde{u}^2(t-\tau) - \epsilon_1\beta\tilde{S}^2(t)I(t), \\ &\quad +\epsilon_1\beta S^0 I(t) \frac{(\mu + a_E)(r + a_I)}{\mu\beta} \left(\mathcal{R}_0 - 1 + \frac{\mu + a_E}{\mu\beta}r\eta\right) \\ &\quad -\epsilon_1 S^0 \frac{\mu + a_E}{\mu}\eta a_R R(t), \\ &\leq -A\tilde{S}^2(t) + B\tilde{S}(t)\tilde{u}(t-\tau) - C\tilde{u}^2(t-\tau), \\ &\quad +\epsilon_1\beta S^0 I(t) \frac{(\mu + a_E)(r + a_I)}{\mu\beta} \left(\mathcal{R}_0 - 1 + \frac{\mu + a_E}{\mu\beta}r\eta\right) \\ &\quad -\epsilon_1 S^0 \frac{\mu + a_E}{\mu}\eta a_R R(t). \end{aligned} $$

This leads to the following inequality

$$ \begin{aligned} \tilde{V}'(t) \leq & -C \left[ \left( \tilde{u}(t-\tau) - \frac{B}{2C} \tilde{S}(t) \right)^2 + \frac{4AC - B^2}{4C^2} \tilde{S}^2(t) \right] \\ & + \epsilon_1 \beta S^0 I(t) \frac{(\mu + a_E)(r + a_I)}{\mu \beta} \left( \mathcal{R}_0 - 1 + \frac{\mu + a_E}{\mu \beta} r \eta \right) \\ & - \epsilon_1 S^0 \frac{\mu + a_E}{\mu} \eta a_R R(t). \end{aligned} \quad (15) $$

The next step is to analyze the sign of $-A\tilde{S}^2 + B\tilde{S}\tilde{u}(t-\tau) - C\tilde{u}^2(t-\tau)$. We compute the discriminant

$$ \begin{aligned} \Delta(\epsilon_1) &= B^2 - 4AC, \\ &= (\epsilon_1 e^{-a_p \tau} (1 - \alpha) + 2h\alpha e^{-a_p \tau})^2 - 4(-\epsilon_1(h + a_S) + h^2)(\alpha^2 e^{-2a_p \tau} - 1). \end{aligned} $$
Next, we compute $\Delta(\epsilon_1)$ and obtain

$$
\begin{align*}
\Delta(\epsilon_1) &= \epsilon_1^2 e^{-2a_p \tau} (1-\alpha)^2 + 4h^2 \alpha^2 e^{-2a_p \tau} - (-4\epsilon_1(h+a_S) + 4h^2)(\alpha^2 e^{-2a_p \tau} - 1) \\
&\quad + 4\epsilon_1 e^{-2a_p \tau} (1-\alpha)h\alpha, \\
&= \epsilon_1^2 e^{-2a_p \tau} (1-\alpha)^2 + 4 \left[ (1-\alpha)h\alpha e^{-2a_p \tau} + (\alpha^2 e^{-2a_p \tau} - 1)(h+a_S) \right] \epsilon_1 \\
&\quad - 4h^2(\alpha^2 e^{-2a_p \tau} - 1) + 4h^2 \alpha^2 e^{-2a_p \tau}, \\
&= \epsilon_1^2 e^{-2a_p \tau} (1-\alpha)^2 + 4 \left[ h\alpha e^{-2a_p \tau} + a_S \alpha^2 e^{-2a_p \tau} - (h+a_S) \right] \epsilon_1 + 4h^2.
\end{align*}
$$

The quantity $\epsilon_1$ defined by (14) is the value that minimize the function $\Delta$. In another side, the discriminant of $\Delta$ given above by a second degree polynomial function is calculated as follow

$$
\begin{align*}
\Delta_{\Theta} &= 16 \left[ (h(\alpha e^{-2a_p \tau} - 1) + a_S((\alpha e^{-a_p \tau})^2 - 1))^2 - h^2(e^{-a_p \tau}(1-\alpha))^2 \right], \\
&= 16 \left[ h(\alpha e^{-2a_p \tau} - 1) + a_S((\alpha e^{-a_p \tau})^2 - 1) + h(e^{-a_p \tau}(1-\alpha)) \right] \\
&\quad \times \left[ h(\alpha e^{-2a_p \tau} - 1) + a_S((\alpha e^{-a_p \tau})^2 - 1) - h(e^{-a_p \tau}(1-\alpha)) \right].
\end{align*}
$$

Obviously, we have

$$
h(\alpha e^{-2a_p \tau} - 1) + a_S((\alpha e^{-a_p \tau})^2 - 1) - h(e^{-a_p \tau}(1-\alpha)) < 0.
$$

Moreover, we get

$$
\begin{align*}
& h(\alpha e^{-2a_p \tau} - 1) + a_S((\alpha e^{-a_p \tau})^2 - 1) + h(e^{-a_p \tau}(1-\alpha)) \\
&= h\alpha e^{-2a_p \tau} - h + a_S((\alpha e^{-a_p \tau})^2 - 1) + he^{-a_p \tau} - \alpha h e^{-a_p \tau}, \\
&= h(e^{-a_p \tau} - 1) + \alpha h e^{-a_p \tau} (e^{-a_p \tau} - 1) + a_S((\alpha e^{-a_p \tau})^2 - 1) < 0.
\end{align*}
$$

Hence, the positivity of $\Delta_{\Theta}$ implies directly that, for $\epsilon_1$ given by (14), we have
$\Delta(\epsilon_1) < 0$. We know that $C > 0$. So, (15) is equivalent to

$$
\begin{align}
\tilde{V}'(t) \le {}& \frac{B^2 - 4AC}{4C} \tilde{S}^2(t) + \epsilon_1 \beta S^0 I(t) \frac{(\mu + a_E)(r + a_I)}{\mu \beta} \left( \mathcal{R}_0 - 1 + \frac{\mu + a_E}{\mu \beta} r \eta \right) \nonumber \\
& - \epsilon_1 S^0 \frac{(\mu + a_E)}{\mu} \eta a_R R(t) =: -\gamma_1 \tilde{S}^2(t) - \gamma_2 I(t) - \gamma_3 R(t), \tag{16}
\end{align}
$$

with

$$
\gamma_1 := \frac{-B^2 + 4AC}{4C} > 0,
$$

$$
\gamma_2 := \epsilon_1 \beta S^0 \frac{(\mu + a_E)(r + a_I)}{\mu \beta} \left( -\mathcal{R}_0 + 1 - \frac{\mu + a_E}{\mu \beta} r \eta \right) > 0,
$$

$$
\gamma_3 := \epsilon_1 S^0 \frac{(\mu + a_E)}{\mu} \eta a_R > 0.
$$

As a first consequence, we can say that the function $\tilde{V}$ is non-increasing and

$$
\tilde{V}(t) := V(\tilde{S}(t), E(t), I(t), R(t), \tilde{u}_t) \xrightarrow{t \to +\infty} \inf_{s \ge 0} [V(\tilde{S}(s), E(s), I(s), R(s), \tilde{u}_s)] := V^* \ge 0.
$$

In addition, by integrating (16), we obtain

$$
\begin{align*}
& \gamma_1 \int_0^t \tilde{S}^2(s) ds + \gamma_2 \int_0^t I(s) ds + \gamma_3 \int_0^t R(s) ds \\
& \le V(\tilde{S}(0), E(0), I(0), R(0), \tilde{u}_0) - V(\tilde{S}(t), E(t), I(t), R(t), \tilde{u}_t). \tag{17}
\end{align*}
$$
By monotonicity, the limits of the expressions in the inequality (17) exist when $t$ tends to $+\infty$. We then obtain

$$ \lim_{t \to +\infty} \int_0^t \tilde{S}^2(s) ds + \lim_{t \to +\infty} \int_0^t I(s) ds + \lim_{t \to +\infty} \int_0^t R(s) ds \\ \le \frac{1}{\min\{\gamma_1, \gamma_2, \gamma_3\}} [V(\tilde{S}(0), E(0), I(0), R(0), \tilde{u}_0) - V^*]. $$

We can see directly from the system (13), that the functions $\tilde{S}'(t)$, $I'(t)$ and $R'(t)$ are uniformly bounded. Then, $\tilde{S}(t)$, $I(t)$ and $R(t)$ are uniformly continuous. We conclude by applying the lemma of Barbalat [26], to the function $t \mapsto \int_0^t \tilde{S}^2(s) ds + \int_0^t I(s) ds + \int_0^t R(s) ds$, that

$$ \lim_{t \to +\infty} \tilde{S}^2(t) + I(t) + R(t) = 0. $$

As a consequence,

$$ \lim_{t \to +\infty} \tilde{S}(t) = 0, \quad \lim_{t \to +\infty} I(t) = 0 \quad \text{and} \quad \lim_{t \to +\infty} R(t) = 0. $$

Since $\lim_{t \to +\infty} \tilde{S}(t) = 0$, then from Lemma 3.5 of [16], we get that

$$ \lim_{t \to +\infty} \tilde{u}(t) = 0. $$

So, the expression of the function $V$ implies that

$$ \lim_{t \to +\infty} E(t) = V^*. \qquad (18) $$

We claim that $V^* = 0$. In fact, the function $I$ is bounded and differentiable. Then, the fluctuations lemma, see Lemma A.14 of [36], ensures the existence of a sequence $t_n \to +\infty$ such that

$$ \lim_{t \to +\infty} I'(t_n) = 0. $$

Then, the third equation of the system (13) implies that

$$ \lim_{t \to +\infty} E(t_n) = 0. $$

Therefore, we get that $V^* = 0$. From (18), we deduce that

$$ \lim_{t \to +\infty} E(t) = 0. $$

The global asymptotic stability of the disease-free equilibrium is ensured. $\square$

# 5 Persistence and global attractivity of the endemic equilibrium

In this section, we want to examine the long-term behavior of the system in the case where $\mathcal{R}_0 > 1$. More precisely, we want to know if the infection persists or not. We start by proving the weak persistence of the disease.
**Lemma 5.1.** Assume that $\mathcal{R}_0 > 1$. Then, there exists a constant $\epsilon > 0$ such that, for any initial condition $(S_0, E_0, I_0, R_0, \phi) \in \mathbb{R}_+ \times (\mathbb{R}_+^2 \setminus \{0, 0\}) \times \mathbb{R}_+ \times C([-\tau, 0], \mathbb{R}_+)$, we have

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

*Proof.* We suppose by contradiction that $\limsup_{t \to +\infty} I(t) \le \epsilon$, for any small $\epsilon > 0$. Then, there exists a sufficiently large $t_\epsilon > 0$ such that $I(t) \le \epsilon$, for all $t \ge t_\epsilon$. Hence, we get for all $t \ge t_\epsilon$,

$$ S'(t) \ge \Lambda - \mu S(t) - \beta \epsilon S(t) + (1 - \alpha)e^{-a_p \tau} u(t - \tau). $$

We put $\liminf_{t \to +\infty} S(t) = S_\infty$ and $\liminf_{t \to +\infty} u(t) = u_\infty$. Then, there exists a sequence $t_m$ tending to infinity as $m$ tends to infinity, such that $S(t_m) \to S_\infty$ and $S'(t_m) \to 0$. This yields to

$$ 0 \ge \Lambda - (h + a_S)S_\infty - \beta S_\infty \epsilon + (1 - \alpha)e^{-a_p \tau} u_\infty. $$

We also obtain

$$ u_\infty \ge \frac{hS_\infty}{1 - \alpha e^{-a_p \tau}}. $$

Then, we have

$$ S_\infty \ge \frac{\Lambda}{a_S + h \frac{1 - e^{-a_p \tau}}{1 - \alpha e^{-a_p \tau}} + \beta \epsilon}. $$

By using the fact that $\mathcal{R}_0 > 1$, we can consider $\epsilon$ and $\lambda > 0$ sufficiently smalls such that

$$ \mathcal{R}_0^{\epsilon, \lambda} := \frac{\beta \mu}{(\lambda + \mu + a_E)(\lambda + r + a_I)} S_\epsilon^0 > 1, $$

where

$$ S_\epsilon^0 := \frac{\Lambda}{a_S + h \frac{1 - e^{-a_p \tau}}{1 - \alpha e^{-a_p \tau}} + \beta \epsilon}. $$

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

$$ -e^{-\lambda t_\epsilon} I(t_\epsilon) > \mu \int_{t_\epsilon}^{+\infty} e^{-\nu t} E(t) dt - (\lambda + r + a_I) \int_{t_\epsilon}^{+\infty} e^{-\nu t} I(t) dt. $$

We also obtain

$$ 0 > -e^{-\lambda t_\epsilon} E(t_\epsilon) > \beta S_\epsilon^0 \int_{t_\epsilon}^{+\infty} e^{-\nu t} I(t) dt - (\lambda + \mu + a_E) \int_{t_\epsilon}^{+\infty} e^{-\nu t} E(t) dt. $$

Therefore, using the fact that $\mathcal{R}_0^{\epsilon, \lambda} > 1$, we arrive to

$$ \begin{aligned} 0 &> -e^{-\lambda t_\epsilon} I(t_\epsilon), \\ &> \frac{1}{(\lambda + \mu + a_E)} [\beta \mu S_\epsilon^0 - (\lambda + \mu + a_E)(\lambda + r + a_I)] \int_{t_\epsilon}^{+\infty} e^{-\nu t} I(t) dt > 0. \end{aligned} $$

This leads to a contradiction. Then, we get the weak uniform persistence. $\square$

Using an idea similar to the one in the article [22], we can prove the strong persistence of the disease. This is done in the following theorem.
**Theorem 5.2.** Assume that $\mathcal{R}_0 > 1$. Then, there exists a constant $\epsilon_1 > 0$ such that

$$ \liminf_{t \to +\infty} I(t) > \epsilon_1 $$

for any initial condition $(S_0, E_0, I_0, R_0, \phi) \in \mathbb{R}_+ \times (\mathbb{R}_+^2 \setminus \{0, 0\}) \times \mathbb{R}_+ \times C([-\tau, 0], \mathbb{R}_+)$.

*Proof.* The proof follows the same ideas as those corresponding to the similar result in [2] and is omitted. $\square$

Next, we study the global stability of the endemic equilibrium $(S^*, E^*, I^*, R^*, u^*)$ of System (1). The strategy to prove the global asymptotic stability of the endemic equilibrium $(S^*, E^*, I^*, R^*, u^*)$ of the system (1), consists in constructing a well adapted Lyapunov functional.

**Theorem 5.3.** Assume that $\mathcal{R}_0 > 1$ and $a_p = 0$. Then, the endemic equilibrium $(S^*, E^*, I^*, R^*, u^*)$ of the system (1) is globally attractive.

*Proof.* Since $S, E, I$ and $u$ are independent on $R$, we can omit the equation of $R$. Assume that $\mathcal{R}_0 > 1$. Then, the endemic equilibrium $(S^*, E^*, I^*, u^*)$ exists. We construct the following Lyapunov functional

$$ W : \mathbb{R}_+^3 \times C([-\tau, 0], \mathbb{R}_+) \to \mathbb{R}_+ \\ (S(t), E(t), I(t), u_t) \to \tilde{W}(t) := W(S(t), E(t), I(t), u_t), $$

defined, for $t \ge 0$, by

$$ \tilde{W}(t) := W_1(t) + W_2(t), \\ W_1(t) := S^*g\left(\frac{S(t)}{S^*}\right) + E^*g\left(\frac{E(t)}{E^*}\right) + \frac{\beta S^*(I^*)^2}{\mu E^*}g\left(\frac{I(t)}{I^*}\right), \\ W_2(t) := \kappa \int_{-\tau}^{0} u^*g\left(\frac{u_t(s)}{u^*}\right) ds, $$

where $g(x) := x - 1 - \ln x$, $x > 0$ and $\kappa$ is a positive constant to be defined below. The differentiation of $W_1$ along the solution trajectory is

$$ W_1'(t) = \left(1 - \frac{S^*}{S}\right) \left[ \Lambda - a_S S - \beta S I - u(t) + e^{-a_p \tau} u(t - \tau) \right] \\ + \left(1 - \frac{E^*}{E}\right) \left[ \beta S I - (\mu + a_E)E \right] + \frac{\beta S^* I^*}{\mu E^*} \left(1 - \frac{I^*}{I}\right) [\mu E - (r + a_I)I]. $$

That is,

$$ W_1'(t) = \left(1 - \frac{S^*}{S}\right) \left\{ a_S(S^* - S) + \beta S^* I^* - \beta S I + u^* - u(t) + e^{-a_p \tau} [u(t - \tau) - u^*] \right\} \\ + \left(1 - \frac{E^*}{E}\right) \beta S I + \left(1 - \frac{E}{E^*}\right) (\mu + a_E)E^* \\ + \frac{\beta S^* I^*}{\mu E^*} \left[ \left(1 - \frac{I^*}{I}\right) \mu E + \left(1 - \frac{I}{I^*}\right) (r + a_I)I^* \right]. $$
Then, we obtain

$$
\begin{aligned}
W_1'(t) = & - \frac{a_S(S - S^*)^2}{S} + \left(1 - \frac{S^*}{S} - \frac{SI}{S^*I^*} + \frac{I}{I^*}\right) \beta S^* I^* \\
& + \left(1 - \frac{S^*}{S}\right) \{u^* - u(t) + e^{-a_p \tau} [u(t - \tau) - u^*]\} \\
& + \left(1 - \frac{E^*}{E}\right) \beta SI + \left(1 - \frac{E}{E^*}\right) \beta S^* I^* \\
& + \frac{\beta S^* I^*}{\mu E^*} \left[\left(1 - \frac{I^*}{I}\right) \mu E + \left(1 - \frac{I}{I^*}\right) \mu E^*\right].
\end{aligned}
$$

This means that

$$
\begin{aligned}
W_1'(t) &= - \frac{a_S(S - S^*)^2}{S} + \left(3 - \frac{S^*}{S} - \frac{SIE^*}{S^*I^*E} - \frac{EI^*}{E^*I}\right) \beta S^* I^* \\
&\quad + \left(1 - \frac{S^*}{S}\right) \{u^* - u(t) + e^{-a_p \tau} [u(t - \tau) - u^*]\}, \\
&= - \frac{a_S(S - S^*)^2}{S} - \left[g\left(\frac{S^*}{S}\right) + g\left(\frac{SIE^*}{S^*I^*E}\right) + g\left(\frac{EI^*}{E^*I}\right)\right] \beta S^* I^* \\
&\quad + \left(1 - \frac{S^*}{S}\right) \{u^* - u(t) + e^{-a_p \tau} [u(t - \tau) - u^*]\}.
\end{aligned}
$$

Let $\hat{u}_t := u_t - u^*$, $t \ge 0$. We then have

$$
\begin{aligned}
& \left(1 - \frac{S^*}{S}\right) \{u^* - u(t) + e^{-a_p \tau} [u(t - \tau) - u^*]\} \\
&= \left(1 - \frac{u^* - \alpha e^{-a_p \tau} u^*}{u(t) - \alpha e^{-a_p \tau} u(t - \tau)}\right) [-\hat{u}(t) + e^{-a_p \tau} \hat{u}(t - \tau)], \\
&= \frac{\hat{u}(t) - \alpha e^{-a_p \tau} \hat{u}(t - \tau)}{u(t) - \alpha e^{-a_p \tau} u(t - \tau)} [-\hat{u}(t) + e^{-a_p \tau} \hat{u}(t - \tau)], \\
&= \frac{\hat{u}(t) - \alpha \delta \hat{u}(t - \tau)}{u(t) - \alpha \delta u(t - \tau)} [-\hat{u}(t) + \delta \hat{u}(t - \tau)],
\end{aligned}
$$

where $\delta := e^{-a_p \tau} \le 1$. On the other hand, we have

$$
\begin{aligned}
W_2'(t) &= \kappa \left[ u^* g \left( \frac{u(t)}{u^*} \right) - u^* g \left( \frac{u(t-\tau)}{u^*} \right) \right], \\
&= \kappa \left[ u(t) - u^* - u^* \ln \frac{u(t)}{u^*} - u(t-\tau) + u^* + u^* \ln \frac{u(t-\tau)}{u^*} \right], \\
&= \kappa \left[ \hat{u}(t) - \hat{u}(t-\tau) + u^* \ln \frac{u(t-\tau)}{u(t)} \right], \\
&= \kappa \left[ \hat{u}(t) - \hat{u}(t-\tau) + u^* \frac{u(t-\tau)}{u(t)} - u^* - u^* g \left( \frac{u(t-\tau)}{u(t)} \right) \right], \\
&= \kappa \left( 1 - \frac{u^*}{u(t)} \right) \left[ \hat{u}(t) - \hat{u}(t-\tau) \right] - \kappa u^* g \left( \frac{u(t-\tau)}{u(t)} \right), \\
&= \kappa \frac{\hat{u}(t)}{u(t)} \left[ \hat{u}(t) - \hat{u}(t-\tau) \right] - \kappa u^* g \left( \frac{u(t-\tau)}{u(t)} \right).
\end{aligned}
$$
Hence, the derivative of $\tilde{W}$ is given, for $t > 0$, by

$$ 
\begin{aligned}
\tilde{W}'(t) &= W_1'(t) + W_2'(t), \\
&= -\frac{a_S(S-S^*)^2}{S} - \left[ g\left(\frac{S^*}{S}\right) + g\left(\frac{SIE^*}{S^*I^*E}\right) + g\left(\frac{EI^*}{E^*I}\right) \right] \beta S^* I^* \\
&\quad - \kappa u^* g\left(\frac{u(t-\tau)}{u(t)}\right) - \frac{\hat{u}(t) - \alpha\delta\hat{u}(t-\tau)}{u(t) - \alpha\delta u(t-\tau)} [\hat{u}(t) - \delta\hat{u}(t-\tau)] \\
&\quad + \kappa \frac{\hat{u}(t)}{u(t)} [\hat{u}(t) - \hat{u}(t-\tau)].
\end{aligned}
$$

Recall that $a_p = 0$. Then, we have $\delta = 1$. Taking $\kappa = 1$, we obtain

$$ 
\begin{aligned}
\tilde{W}'(t) &= -\frac{a_S(S-S^*)^2}{S} - \left[ g\left(\frac{S^*}{S}\right) + g\left(\frac{SIE^*}{S^*I^*E}\right) + g\left(\frac{EI^*}{E^*I}\right) \right] \beta S^* I^* \\
&\quad - u^* g\left(\frac{u(t-\tau)}{u(t)}\right) + \left[ -\frac{\hat{u}(t) - \alpha\hat{u}(t-\tau)}{u(t) - \alpha u(t-\tau)} + \frac{\hat{u}(t)}{u(t)} \right] [\hat{u}(t) - \hat{u}(t-\tau)], \\
&= -\frac{a_S(S-S^*)^2}{S} - \left[ g\left(\frac{S^*}{S}\right) + g\left(\frac{SIE^*}{S^*I^*E}\right) + g\left(\frac{EI^*}{E^*I}\right) \right] \beta S^* I^* \\
&\quad - u^* g\left(\frac{u(t-\tau)}{u(t)}\right) + \alpha \frac{u(t)\hat{u}(t-\tau) - u(t-\tau)\hat{u}(t)}{[u(t) - \alpha u(t-\tau)] u(t)} [\hat{u}(t) - \hat{u}(t-\tau)], \\
&= -\frac{a_S(S-S^*)^2}{S} - \left[ g\left(\frac{S^*}{S}\right) + g\left(\frac{SIE^*}{S^*I^*E}\right) + g\left(\frac{EI^*}{E^*I}\right) \right] \beta S^* I^* \\
&\quad - u^* g\left(\frac{u(t-\tau)}{u(t)}\right) - \frac{\alpha u^*}{hS(t)u(t)} [u(t) - u(t-\tau)]^2, \\
&\le 0.
\end{aligned}
$$

Hence, $W$ is a Lyapunov function. The last part of the proof of the theorem 4.1 can be adapted to show the convergence of any solution to the endemic equilibrium $(S^*, E^*, I^*, u^*)$. $\square$

**Proposition 2.** Assume that $\mathcal{R}_0 > 1$ and $\alpha = 0$. Then, the endemic equilibrium $(S^*, E^*, I^*, R^*, u^*)$ of the system (1) is globally attractive.

*Proof.* We consider the same functional $W$ used in the proof of the last theorem, with $\alpha = 0$ and $\kappa = e^{-a_p \tau}$,

$$ \tilde{W}(t) := S^* g\left(\frac{S(t)}{S^*}\right) + E^* g\left(\frac{E(t)}{E^*}\right) + \frac{\beta S^*(I^*)^2}{\mu E^*} g\left(\frac{I(t)}{I^*}\right) + \kappa \int_{-\tau}^{0} u^* g\left(\frac{u_t(s)}{u^*}\right) ds. $$

We obtain

$$ 
\begin{aligned}
\tilde{W}'(t) &= -\frac{a_S(S-S^*)^2}{S} - \left[ g\left(\frac{S^*}{S}\right) + g\left(\frac{SIE^*}{S^*I^*E}\right) + g\left(\frac{EI^*}{E^*I}\right) \right] \beta S^* I^* \\
&\quad - e^{-a_p \tau} u^* g\left(\frac{u(t-\tau)}{u(t)}\right) - \frac{(1-e^{-a_p \tau})}{u(t)} [\hat{u}(t)]^2, \\
&\le 0.
\end{aligned}
$$

Hence, $W$ is a Lyapunov function.

<figure><img src="image_3.png" /></figure>
# 6 Discussion and numerical simulations

In this paper, we present and study a SEIR epidemiological model with a class of individuals temporarily protected, through vaccination or medication for instance. In the literature, there are many epidemiological models with a compartment of individuals temporarily protected from infection (see for instance [11, 15, 19, 20, 23, 31]). In most cases, these systems are modeled directly by delay differential equations, or they are obtained from age-structured PDEs and reduced by integration using the method of characteristics to delay differential systems, where the delay generally corresponds to the duration of the protection phase. A restrictive assumption is often made to reduce the model to a delay differential system: the entire protected population becomes susceptible at the end of the protection phase. However, for many infectious diseases, such as hepatitis A and B, human papillomavirus, influenza, COVID 19, etc., people who renew their vaccine remain protected, and only those who don't renew their vaccine become susceptible. In this case, the model can no longer be reduced to a classical delay differential system. The resulting model as we have shown, is a coupled system between differential equations and a difference equation. With the exception of our works [2, 5], this approach has never been used in epidemiological models. A novelty that our model can bring compared to what exists in the literature is that we can reduce the basic reproduction number $R_0$ by increasing the new parameter $\alpha$, and thus reduce the epidemic (see Figures 2 and 3). This observation is interesting from an epidemiological point of view, as it is generally easier to increase the rate of renewal of protection for individuals already protected (because this population is already under medical care) than to increase the rate of protection $h$ or the duration $\tau$.

In addition, we consider here a more general class of Kermack-McKendrick models than the one considered in [2], adding the compartment of exposed individuals $E$. Indeed, for many infections, such as influenza [21] and tuberculosis [17], there is a more or less long latency period during which individuals have been infected but are not yet infectious themselves. During the period $1/\mu$, infected individuals are in the compartment $E$ but are not yet infectious (see also [32] for other cases of infectious diseases with a latency period).

The model we consider in this work, is a system composed of differential equations and a difference equation with delay. In addition to the disease-free equilibrium, it admits an endemic steady-state. We focus on the global asymptotic stability of the two steady-states. The construction of Lyapunov functions allows us to prove the global asymptotic stability of both equilibria. More precisely, by using the basic reproduction number $R_0$ as a threshold, we obtain a necessary and sufficient condition for the global asymptotic stability of the disease-free equilibrium. Furthermore, we show that the endemic equilibrium is globally attractive if $R_0 > 1$ and there is no mortality during the protection phase. To complete the study of the global dynamics of the system, we show results on the uniform weak and strong persistence of the disease.

We present some numerical simulations to illustrate the results obtained in the previous sections. For the simulations, we set the parameters as follows

$$ \Lambda = 2, \mu = 0.15, r = 0.15, a_S = a_E = a_I = 0.1, \beta = 0.1 \text{ and } h = 0.5. \quad (19) $$

We choose to modify the parameters $\alpha$ and $T := a_p\tau$. For some figures, we also change the values of $\beta$ and $h$. Figure 2 shows the variation of $R_0$ as functions
of the parameters $\alpha$ and $T := a_p\tau$. We see that $\mathcal{R}_0$ is a decreasing function with respect to $\alpha$ and $T$. Moreover, by increasing the value of the parameter $h$, the region $\mathcal{R}_0(\alpha, T) < 1$ increases. From an epidemiological point of view, this means that the epidemic could be stopped, or at least slowed down, by increasing one of the three parameters $h$, $\tau$ or $\alpha$. Combining the three parameters offers more possibilities for controlling the epidemic, as in reality some parameters, such as $\tau$ or sometimes even $h$, can no longer be increased. Figure 3 shows the variation of the equilibrium of infected individuals $I^*$ as a function of the parameters $\alpha$ and $T := a_p\tau$. We can see that it is a decreasing function with respect to each of these two parameters. This shows that even when $\mathcal{R}_0$ remains above 1, the parameters $\tau$ and $\alpha$ reduce the epidemic by decreasing the number of infected individuals $I^*$ at equilibrium. Figure 4 shows the global convergence of the solution to the disease-free equilibrium $(S^0, 0, 0, u^0)$ and Figure 5 illustrates the global stability of the endemic steady-state $(S^*, E^*, I^*, u^*)$.

<figure><img src="image_4.png" /><figcaption>Figure 2: Representation of the curve $\mathcal{R}_0(\alpha, T) = 1$, with $T := a_p\tau$, by varying the parameter $h$.</figcaption></figure>

Although we proved the global stability of the endemic equilibrium when $\mathcal{R}_0 > 1$ only for the particular cases $\alpha = 0$ or $a_p = 0$. These are technical conditions for obtaining a Lyapunov function. It is possible that the endemic equilibrium remains globally stable even in the cases $\alpha > 0$ and $a_p > 0$ (this is confirmed by numerical simulations and by the persistence property proved in the general case $\mathcal{R}_0 > 1$). The assumption $a_p = 0$ means that there is no mortality during the protection phase, and the assumption $\alpha = 0$ means that at the end of the protection phase, all individuals become susceptible.

## Conflict of interest

This work does not have any conflicts of interest.
<figure><img src="image_5.png" /><figcaption>Figure 3: Representation of the equilibrium of infected individuals *I*\* as a function of the parameters α and T := *a*<sub>p</sub>τ.</figcaption></figure>

**Acknowledgement(s)**

A. Chekroun thanks grant, PRFU: C00L03UN29012022002, from DGRSDT of Algeria. T. Kuniya is supported by JSPS Grant-in-Aid for Early-Career Scientists [grant number 19K14594] and the Japan Agency for Medical Research and Development [grant number JP20fk0108535].

**References**

[1] M. Adimy, A. Chekroun, A. El Abdllaoui, and A. Marzorati. Discrete maturity and delay differential-difference model of hematopoietic cell dynamics with applications to acute myelogenous leukemia. *Journal of Biological Systems*, (2022) 1-31.

[2] M. Adimy, A. Chekroun, and C. P. Ferreira. Global dynamics of a differential-difference system: a case of Kermack-McKendrick SIR model with age-structured protection phase. *Mathematical Biosciences and Engineering, 17* (2020), 1329-1354.

[3] M. Adimy, A. Chekroun, and B. Kazmierczak. Traveling waves in a coupled reaction-diffusion and difference model of hematopoiesis. *Journal of Differential Equations, 262* (2017), 4085-4128.
<figure><img src="image_6.png" /><figcaption>Figure 4: The disease-free equilibrium is globally asymptotically stable for the parameter values given by (19), with in addition $\beta = 0.05$, $\alpha = 0.2$, $a_p = 0.1$ and $\tau = 4$. In this case, $R_0 = 0.4 < 1$. The initial condition is set to $S_0 = 0.75$, $I_0 = 0.15$, $E_0 = 0.1$ and $u(\theta) = 0.2$, for $\theta \in [-\tau, 0]$.</figcaption></figure>

[4] M. Adimy, A. Chekroun, and T. Kuniya. Coupled reaction-diffusion and difference system of cell-cycle dynamics for hematopoiesis process with Dirichlet boundary conditions. *Journal of Mathematical Analysis and Applications*, **479** (2019), 1030-1068.

[5] M. Adimy, A. Chekroun, and T. Kuniya. Global asymptotic stability for a distributed delay differential-difference system of a Kermack-McKendrick SIR model. *Applicable Analysis*, (2022), 1-13.

[6] M. Adimy, A. Chekroun, and T. Kuniya. Traveling waves of a differential-difference diffusive Kermack-McKendrick epidemic model with age-structured protection phase. *Journal of Mathematical Analysis and Applications*, **505** (2022), 125464.

[7] M. Adimy, A. Chekroun, and T.-M. Touaoula. Age-structured and delay differential-difference model of hematopoietic stem cell dynamics. *Discrete and Continuous Dynamical Systems - Series B*, **20** (2015), 2765-2791.

[8] M. Adimy, A. Chekroun, and T.-M. Touaoula. A delay differential-difference system of hematopoietic stem cell dynamics. *Comptes Rendus Mathematique*, **353** (2015), 303-307.

[9] L. J. Allen and P. V. den Driessche. The basic reproduction number in some discrete-time epidemic models. *Journal of Difference Equations and Applications*, **14** (2008), 1127-1147.
<figure><img src="image_7.png" /><figcaption>Figure 5: The global stability of the endemic equilibrium for the parameter values given by (19), with in addition $\beta = 0.2$, $\alpha = 0.2$, $a_p = 0.1$ and $\tau = 4$. In this case, $R_0 = 1.6 > 1$. The initial condition is set to $S_0 = 0.75$, $I_0 = 0.15$, $E_0 = 0.1$ and $u(\theta) = 0.2$ for $\theta \in [-\tau, 0]$.</figcaption></figure>

[10] B. Balachandran, T. Kalmar-Nagy, and D. E. Gilsinn. *Delay Differential Equations: Recent Advances and New Directions*. Springer-Verlag New York Inc., New York, 2009.

[11] E. Beretta, T. Hara, W. Ma, and Y. Takeuchi. Global asymptotic stability of an SIR epidemic model with distributed time delay. *Nonlinear Analysis, 47* (2001), 4107-4115.

[12] E. Beretta and Y. Takeuchi. Global stability of an SIR epidemic model with time delays. *Journal of Mathematical Biology, 33* (1995), 250-260.

[13] F. Brauer and C. Castillo-Chavez. *Mathematical models in population biology and epidemiology*. Texts in Applied Mathematics, Springer, New York, 2nd edition, 2012.

[14] V. Capasso and G. Serio. A generalization of the Kermack-McKendrick deterministic epidemic model. *Mathematical Biosciences, 42* (1978), 43-61.

[15] K. L. Cooke and P. van den Driessche. Analysis of an SEIRS epidemic model with two delays. *Journal of Mathematical Biology, 35* (1996), 240-260.

[16] M. A. Cruz and J. K. Hale. Stability of functional differential equations of neutral type. *Journal of Differential Equations, 7* (1970), 334-355.
[17] K. Das, B.S.N. Murthy, S.A. Samad, and M.H.A. Biswas. Mathematical transmission analysis of SEIR tuberculosis disease model. *Sensors International*, **7** (2021), 100120.

[18] O. Diekmann, J. A. P. Heesterbeek, and J. A. J. Metz. On the definition and the computation of the basic reproduction ratio $R_0$ in models for infectious diseases in heterogeneous populations. *Journal of Mathematical Biology*, **28** (1990), 365-382.

[19] Y. Enatsu, Y. Nakata, and Y. Muroya. Global stability of SIR epidemic models with a wide class of nonlinear incidence rates and distributed delays. *Discrete & Continuous Dynamical Systems - B*, **15** (2011), 61-74.

[20] Y. Enatsu, Y. Nakata, and Y. Muroya. Lyapunov functional techniques for the global stability analysis of a delayed SIRS epidemic model. *Nonlinear Analysis: Real World Applications*, **13** (2012), 2120-2133.

[21] F. Etbaigha, A.R. Willms, and Z. Poljak. An SEIR model of influenza A virus infection and reinfection within a farrow-to-finish swine farm. *PLoS One*, **13**(9), (2018), e0202493.

[22] H. I. Freedman and P. Moson. Persistence definitions and their connections. *Proceedings of the American Mathematical Society*, **109** (1990), 1025-1033.

[23] S. Gao, L. Chen, and Z. Teng. Pulse vaccination of an SEIR epidemic model with time delay. *Nonlinear Analysis: Real World Applications*, **9** (2008), 599-607.

[24] K. Gu and Y. Liu. Lyapunov-Krasovskii functional for uniform stability of coupled differential-functional equations. *Automatica*, **45** (2009), 798-804.

[25] J. K. Hale and S. M. Verduyn Lunel. *Introduction to Functional Differential Equations*. Springer, 1993.

[26] M. Hou, G. Duan, and M. Guo. New versions of Barbalat's lemma with applications. *Journal of Control Theory and Applications*, **8** (2010), 545-547.

[27] W. O. Kermack and A. G. McKendrick. A contribution to the mathematical theory of epidemics. *Proc. R. Soc. Lond. A*, **115** (1927), 700-721.

[28] Y. Kuang. *Delay Differential Equations: With Applications in Population Dynamics*. Academic Press, 1993.

[29] M. Y. Li. *An Introduction to Mathematical Modeling of Infectious Diseases*. Mathematics of Planet Earth, Mathematics of Planet Earth Collection. Springer International Publishing, 2018.

[30] M. Y. Li, Z. Shuai, and C. Wang. Global stability of multi-group epidemic models with distributed delays. *Journal of Mathematical Analysis and Applications*, **361** (2010), 38-47.

[31] Q. Liu, D. Jiang, T. Hayat, and B. Ahmad. Analysis of a delayed vaccinated SIR epidemic model with temporary immunity and lévy jumps. *Nonlinear Analysis: Hybrid Systems*, **27** (2018), 29-43.
[32] M. Martcheva. *An Introduction to Mathematical Epidemiology*. Springer Publishing Company, Incorporated, 1st edition, 2015.

[33] C. C. McCluskey. Global stability of an SIR epidemic model with delay and general nonlinear incidence. *Mathematical Biosciences and Engineering, 7* (2010), 837-850.

[34] P. Pepe. On the asymptotic stability of coupled delay differential and continuous time difference equations. *Automatica, 47* (2005), 107-112.

[35] P. Pepe, Z.-P. Jiang, and E. Fridman. A new Lyapunov-Krasovskii methodology for coupled delay differential and difference equations. *International journal of control, 81* (2008), 107-115.

[36] H. L. Smith and H. R. Thieme. *Dynamical Systems and Population Persistence*. Americal Mathematical Society, 2011.

[37] S. Tipsri and W. Chinviriyasit. The effect of time delay on the dynamics of an SEIR model with nonlinear incidence. *Chaos, Solitons & Fractals, 75* (2015), 153-172.

[38] P. Van den Driessche and J. Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. *Mathematical Biosciences, 180* (2002), 29-48.

[39] R. Xu, Z. Wang, and F. Zhang. Global stability and hopf bifurcations of an SEIR epidemiological model with logistic growth and time delay. *Applied Mathematics and Computation, 269* (2015), 332-342.