# Global Asymptotic Stability of a Hybrid Differential-Difference System Describing SIR and SIS Epidemic Models with a Protection Phase and a Nonlinear Force of Infection

Mostafa Adimy, Abdennasser Chekroun, Charlotte Dugourd-Camus, Hanene Meghelli

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

► To cite this version:

Mostafa Adimy, Abdennasser Chekroun, Charlotte Dugourd-Camus, Hanene Meghelli. Global Asymptotic Stability of a Hybrid Differential-Difference System Describing SIR and SIS Epidemic Models with a Protection Phase and a Nonlinear Force of Infection. Qualitative Theory of Dynamical Systems, 2023, 23 (34), pp.32. 10.1007/s12346-023-00891-z. hal-04271594

HAL Id: hal-04271594

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Global asymptotic stability of a hybrid
differential-difference system describing SIR and
SIS epidemic models with a protection phase and
a nonlinear force of infection

Mostafa Adimy<fnref n="1" />*†, Abdennasser Chekroun<fnref n="2" />†,
Charlotte Dugourd-Camus<fnref n="1" />†, Hanene Meghelli<fnref n="2" />†

Abstract

We study the local and global asymptotic stability of the two steady-states,
disease-free and endemic, of hybrid differential-difference SIR and SIS epidemic
models with a nonlinear force of infection and a temporary phase of protection
against the disease, e.g. by vaccination or medication. The initial model is an
age-structured system that is reduced using the method of characteristic lines
to a hybrid system, coupled between differential equations and a time continu-
ous difference equation. We first prove that the solutions of the original system
can be obtained from the reduced one. We then focus on the reduced system
to obtain new results on the asymptotic stability of the two steady-states. We
determine the local asymptotic stability of the two steady-states by studying the
associated characteristic equation. We then discuss their global asymptotic sta-
bility in various situations (SIR, SIS, mass action, nonlinear force of infection),
by constructing appropriate Lyapunov functions.

<footnote n="1">¹*Inria, Univ Lyon, Université Lyon 1, CNRS UMR 5208, Institut Camille
Jordan, 43 Bd. du 11 novembre 1918, Villeurbanne, F-69200, France.</footnote>

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

<footnote>*Corresponding author(s). E-mail(s): mostafa.adimy@inria.fr;
Contributing authors: chekroun@math.univ-lyon1.fr;
charlotte.dugourd@inria.fr ; mhanene97@gmail.com;
†These authors contributed equally to this work.</footnote>
**Keywords:** Local and global asymptotic stability, Lyapunov function, age-structured model, SIR and SIS epidemic models, nonlinear force of infection, hybrid differential-difference system.

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

# 1 Introduction

Among the pioneering mathematical works describing infectious diseases are those by D. Bernoulli (who studied the case of smallpox [13]), W.O. Kermack, A.G. McKendrick (at the origin of the compartmental epidemic models [31]) and R. Ross (who modeled the transmission of malaria [42]). Since then, mathematical modeling has become essential for analyzing the spread of infectious diseases. It has continued to produce sophisticated new tools and systems to answer questions linked to understanding the dynamics, control or eradication of epidemic diseases. The following are some of the most well-known models [9, 10, 16, 20, 22, 23, 37, 39, 40, 41]. In [9], the author revisited the classic model of Kermack-McKendrick [31], and focused on the heterogeneity of transmission induced by a variety of factors. Different types of models (discrete and continuous population systems, models with delays, structured population models, etc.) are detailed in the books [14, 21, 39, 41]. In [16], a generalization of the Kermack-McKendrick model [31], using a nonlinear force of infection that can account for saturation effects, was performed. In [23], the authors established the global asymptotic stability of SIRS epidemic models with nonlinear forces of infection and distributed time delays. The paper [38] studied the Lyapunov function and global asymptotic stability for an age-structured epidemic model using semigroup theory. Age-structured epidemic models were also studied in [17, 18]. The authors provided a mathematical analysis including global asymptotic stability using Lyapunov functions. In most of the above references, the models are constituted of compartments: susceptible, exposed, infectious, removed, vaccinated or protected individuals, etc. Acronyms are often used to describe these models, such as SIR, SIS, SEIR, etc.

A qualitative study of the dynamics of these epidemiological models was performed in detail (see, for instance, [8, 14, 16, 29, 33, 39] and the references cited therein). The next-generation matrix approach is one of the most popular techniques, [8, 21], used in compartmental epidemiological models to calculate the famous basic reproduction number $R_0$. Delayed differential systems have also been largely used in this domain (see, for instance, [10, 12, 19, 23, 24, 27, 34, 36, 40, 44, 46]). Cooke and van den Driessche propose in their paper [19] a SEIRS model with two different time delays, namely a latency period and a period of temporary immunity. In the paper [27], a SEIR model with delays takes into account an impulse vaccination strategy. The authors of this article first analyze a deterministic SIR delayed epidemic model with temporary immunity through vaccination, then consider a delayed model with Lévy jumps. Other models with delays and a nonlinear force of infection have been investigated by McCluskey, [40]. The authors of [44] also considered a nonlinear force
of infection, showing that delay can destabilize the system and provide periodic solu-
tions. The papers [10, 12, 23, 24, 34, 46], investigated the global asymptotic stability
of epidemic models with time delays.

There are also many models that investigated the effects of the protection period,
for example through vaccination or treatment, see [7, 15, 26, 37, 45]. Some infectious
disease models structure the vaccination period according to age (age here being the
time at which the vaccine is administered), with a time-limited duration of immunity,
[22, 30, 35]. Age-structured epidemic models are mainly based on age-structured popu-
lation models of the Lotka-McKendrick type, [39] and the references cited therein. The
local and global asymptotic stability of steady-states are a fundamental and difficult
problem for age-structured epidemic models, [22, 35, 38, 39].

Recently, we have analyzed hybrid differential-difference models with discrete and distributed delay [4, 5], focusing on hematopoiesis modeling. These systems are couplings between differential equations and time-continuous delay difference equations. This approach facilitates the study of the asymptotic behavior of the system, in particular they have the possibility to highlight particular behaviors such as the global asymptotic stability, the existence of periodic solutions, see [4]. More recently, we continued analyzing hybrid differential-difference models with discrete and distributed delay, with applications to epidemiology [1, 2]. The main result is the study of the global asymptotic stability of the two steady-states. In our previous articles, [1, 2], we only considered incidence functions in the form of mass action (3).

In this paper, we focus mainly on SIR and SIS epidemic models with a protection phase and a nonlinear force of infection satisfying the conditions (5). For example, vaccinated people may be protected for a limited period, but they can also lose their immunity if they fail to renew their vaccine. Our aim in this work is to focus on the protection phase, which we will describe in more detail than the other compartments. It is also important to point out that, on the one hand, the use of this type of nonlinear hybrid model gives rise to new and interesting theoretical results for the mathematical community and, on the other hand, we can recover the solution of our original model from this hybrid model.

The paper is organized as follows. The model is formulated in Section 2. In Section 3, we discuss the existence of the two steady-states, disease-free and endemic, of System (14) as a function of the basic reproduction number, $\mathcal{R}_0$. In Section 4, the local asymptotic stability of the two steady-states is investigated by developing a characteristic equation. Section 5 is devoted to the global asymptotic stability of the two steady-states by constructing appropriate Lyapunov functions. We treat two cases, an SIS system, $\alpha > 0$, with bilinear incidence (mass action) and an SIR system, $\alpha = 0$, with a nonlinear force of infection $g$ satisfying the conditions (5). The last section concludes with numerical simulations and a general discussion.

## 2 Model formulation

Newly-susceptible individuals include newborns, immigrants, etc., with total rate $b$, or people whose immunity has been lost after being infected, with a rate $\alpha$ (SIS), or vaccinated people who have not updated their vaccine, with a rate $1 - \nu$, after
a supposed fixed period of time $\tau$. Then, the protection renewal rate is $\nu \in (0, 1)$. Infected individuals can develop permanent immunity with a rate $\mu$ and pass into the compartment of recovered individuals (SIR). Each of the four compartments has a mortality rate denoted by $a_S$, $a_p$, $a_I$, $a_R$. We assume that susceptible individuals move to the protected compartment, for example by being vaccinated, at a rate $\eta$.

Let us consider now the compartment of protected individuals. We denote by $a$ the time since which an individual is temporarily protected, where the maximum duration is equal to $\tau > 0$. This means that the protected individuals are distributed by age over the protected compartment, denoted $p(t, a)$, for $a \in [0, \tau]$ (see [1, 2] for more details). Then, the total population of protected individuals is given by

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

Individuals entering the protection phase, $a = 0$, come from the susceptible compartment with a rate $\eta$, or from individuals that have reached the end of their protection phase, $a = \tau$, and are renewing it with a rate $\nu$. Then,

$$p(t, 0) = \eta S(t) + \nu p(t, \tau). \qquad (1)$$

The evolution of individuals in the protected compartment is given by the following age-structured equation

$$\frac{\partial p}{\partial t}(t, a) + \frac{\partial p}{\partial a}(t, a) = -a_p p(t, a), \quad 0 < a < \tau, \qquad (2)$$

with the boundary condition (1). The force of infection is defined by the function $g(I)$, so that the number of infected individuals per unit time is equal to $g(I)S$ (incidence). There are different forms of infection force, depending on the assumptions made about disease transmission. One of these forms is mass action

$$g(I) = \beta I, \quad \text{with } \beta > 0. \qquad (3)$$

The second form is the saturating force of infection, i.e. the force of infection increases with the number of infected individuals, but reaches a threshold when the number of infected individuals becomes very large. A reasonably simple form of this function is

$$g(I) = \frac{\beta I}{1 + aI}, \quad \beta, a > 0. \qquad (4)$$

In our study, we consider the two cases, mass action where the function $g$ satisfies (3) or any function $g$ satisfying

$$g(0) = 0, \quad g'(I) > 0, \quad g''(I) < 0 \quad \text{for all } I \ge 0. \qquad (5)$$
This last is more general than the form (4). We can see that the case of mass action (3) satisfies only part of the conditions (5), since $g'' = 0$. The conditions (5) give rise to the following properties.

**Lemma 1.** Assume that the function $g$ satisfies the conditions (5). Then, we have the following statements:

(i) $g(I) < g'(0)I$, for all $I > 0$;

(ii) $g'(I) \le \frac{g(I)}{I}$, for all $I > 0$;

(iii) For an arbitrary number $I^* > 0$, the following inequality holds:

$$ \left[ 1 - \frac{g(I)}{g(I^*)} \right] \left[ \frac{Ig(I^*)}{I^*g(I)} - 1 \right] < 0, \quad \text{for } I > 0 \text{ and } I \neq I^*. $$

*Proof.* For the proof of (i) and (iii), see Lemma 3.10 of [25]. The statement (ii) is equivalent to $G(I) := Ig'(I) - g(I) \le 0$, for all $I > 0$. It is not difficult to see that the derivative $G'(I) = Ig''(I) \le 0$. We conclude that $G(I) := Ig'(I) - g(I) \le G(0) = 0$. This finish the proof. $\square$

<figure><img src="image_3.png" /><figcaption>**Fig. 1:** Example of a function $g$ satisfying (5) (solid line) and the corresponding functions $I \mapsto g'(I)I$ (dashed line), $I \mapsto g'(0)I$ (dotted line), with $\beta = 0.015$ and $a = 0.6$.</figcaption></figure>

Figure 1 illustrates an example of a function $g$ satisfying statements (i), (ii) and (iii) of Lemma 1.
To complete our model, and in line with what was explained above, we have the flow of new susceptible individuals at each time $t > 0$, given by

$$b + \alpha I(t) + (1 - \nu)p(t, \tau). \qquad (6)$$

The equations for susceptible, infected and recovered individuals are therefore as follows

$$\left\{
\begin{aligned}
S'(t) &= b - (a_S + \eta + g(I(t)))S(t) + \alpha I(t) + (1 - \nu)p(t, \tau), \\
I'(t) &= g(I(t))S(t) - (a_I + \alpha + \mu)I(t), \\
R'(t) &= -a_R R(t) + \mu I(t).
\end{aligned}
\right.
\qquad (7)$$

With the equations (1), (2) and the initial condition

$$S(0) = S_0, \quad p(0, a) = p_0(a), \quad a \in (0, \tau), \quad I(0) = I_0, \quad R(0) = R_0, \qquad (8)$$

the model is now established. In fact, we have combined the two cases, SIS and SIR, into a single system. With the choice of the parameters $\alpha$ and $\mu$ positive or equal to zero, we can easily recover the results for our particular scenarios, SIR and SIS (see Figure 2).

As an example, by choosing $\alpha = 0$ and $\mu > 0$, like the diagram on the left of Figure 2, we could model chickenpox. An SIS model ($\alpha > 0$ and $\mu = 0$, see the diagram on the right of Figure 2) could be used to model influenza. Particular attention will also be paid to the different types of infection force $g$. In the analysis, we do not explicitly treat the case $\mu = 0$ for the SIS model, as this amounts to considering that the mortality rate of infected individuals is $a_I + \mu$.

<figure><img src="image_4.png" /><figcaption>Fig. 2: Schematic representation of the SIR and SIS models with protection phase.</figcaption></figure>

Our first objective is to reduce the model (1), (2) and (7) to a hybrid differential-difference system that is easier to investigate. For this purpose, we use the method of
characteristic lines (See for instance [30] or [39]) to solve Equation (2) and obtain

$$p(t, a) = \begin{cases} e^{-a_p a} p(t - a, 0), & t > a, \\ e^{-a_p t} p_0(a - t), & t < a. \end{cases} \qquad (9)$$

Consider the new state

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

It represents the new individuals in the protected compartment, at $a = 0$. The function $p(t, a)$ is completely determined by the function $u$ and the initial condition $p_0$. In particular, the total population of protected individuals is given, for $t > \tau$, by

$$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.$$

In the model (1), (2) and (7), only the cases $a = 0$ and $a = \tau$ are involved. We define the function

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

Then, we have

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

Furthermore, the boundary condition (1) becomes

$$u(t) = \begin{cases} \eta S(t) + \nu e^{-a_p \tau} u(t - \tau), & t > 0, \\ \phi(t), & -\tau < t < 0. \end{cases}$$

In particular, for $0 < t < \tau$, the total population of protected individuals is given by

$$P(t) = \int_{0}^{t} e^{-a_p a} u(t - a) da + \int_{t}^{\tau} e^{-a_p a} e^{-a_p (t-a)} p_0(a - t) da, \\ = \int_{0}^{t} e^{-a_p a} u(t - a) da + \int_{t}^{\tau} e^{-a_p a} \phi(t - a) da.$$

We conclude that for any $t > 0$, the total population of protected individuals is

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

In summary, we have reduced the model (1), (2) and (7) to the following hybrid system

$$\left\{ \begin{aligned} S'(t) &= b - (a_S + \eta + g(I(t)))S(t) + \alpha I(t) + (1 - \nu)e^{-a_p \tau} u(t - \tau), \\ I'(t) &= g(I(t))S(t) - (a_I + \alpha + \mu)I(t), \\ R'(t) &= -a_R R(t) + \mu I(t), \\ u(t) &= \eta S(t) + \nu e^{-a_p \tau} u(t - \tau), \end{aligned} \right. \qquad (10)$$
with non-negative initial condition

$$S(0) = S_0, \quad I(0) = I_0, \quad R(0) = R_0, \quad u(t) = \phi(t), \quad t \in [-\tau, 0]. \qquad (11)$$

The model (10) is hybrid in the sense that it is composed of a system of differential equations and a time-continuous difference equation.

The assumption that $p_0$ is integrable on $(0, \tau)$ is sufficient for the existence of (weak) solutions. However, for our qualitative study, we consider more regularities on $p_0$. We assume that $p_0$ is continuous on $[0, \tau]$ and satisfies the compatibility condition

$$p_0(0) = \eta S_0 + \nu p_0(\tau).$$

This means that

$$\phi \in C([-\tau, 0]) \quad \text{and} \quad \phi(0) = \eta S_0 + \nu e^{-a_p \tau} \phi(-\tau). \qquad (12)$$

We then get the following result.

**Theorem 2.** Assume that the initial condition (11) is nonnegative and satisfies the compatibility condition (12). Then, the system (10)-(11) admits a unique solution $(S, I, R, u)$, where $(S, I, R)$ has a continuous first derivative for all $t > 0$ and the function $u$ is continuous for all $t \ge -\tau$. Moreover, all solutions are bounded and nonnegative.

*Proof.* The existence, uniqueness, positivity and regularity of the solution can be proved by steps on the intervals $[(k-1)\tau, k\tau]$, $k = 1, 2, \dots$. To prove that the solution is bounded, we first note that

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

Consider the total population $N(t) := S(t) + I(t) + R(t) + P(t)$. Then, we have for $t > 0$,

$$N'(t) = b - a_S S(t) - a_I I(t) - a_R R(t) - a_p P(t) \le b - \xi N(t),$$

with

$$\xi = \min\{a_S, a_I, a_R, a_p\} > 0.$$

So,

$$N(t) \le \max\{N(0), b/\xi\}.$$

Therefore, the total population $N(t)$ is bounded. Thus, each of the sub-populations $S(t)$, $I(t)$, $R(t)$ and $u(t)$ is bounded. □

Using the expression (9), we obtain results similar to those of Theorem 2 for the model (1), (2), (7) and (8). With additional regularity assumptions on the function $p_0$, we can obtain classical solutions. We can also relax the assumptions on the function $p_0$ to obtain weaker solutions.

**Remark 1.** Here are some particular cases covered by the model (10).
• If $\nu = 0$, the hybrid system (10) becomes the following SIR model with delay

$$ egin{cases} S'(t) = b - (a_S + \eta + g(I(t)))S(t) + \alpha I(t) + \eta e^{-a_p \tau} S(t - \tau), \\ I'(t) = g(I(t))S(t) - (a_I + \alpha + \mu)I(t), \\ R'(t) = -a_R R(t) + \mu I(t). \end{cases} $$

• If $\tau = 0$, we obtain the ordinary differential system

$$ egin{cases} S'(t) = b - (a_S + g(I(t)))S(t) + \alpha I(t), \\ I'(t) = g(I(t))S(t) - (a_I + \alpha + \mu)I(t), \\ R'(t) = -a_R R(t) + \mu I(t). \end{cases} $$

Since in System (10) the equations for $S$, $I$ and $u$ do not depend on $R$, we omit in our study the equation for $R$. Then, the model becomes

$$ egin{cases} S'(t) = b - (a_S + \eta + g(I(t)))S(t) + \alpha I(t) + (1 - \nu)e^{-a_p \tau} u(t - \tau), \\ I'(t) = g(I(t))S(t) - (a_I + \alpha + \mu)I(t), \\ u(t) = \eta S(t) + \nu e^{-a_p \tau} u(t - \tau), \end{cases} \quad (14) $$

with the initial condition

$$ S(0) = S_0, \quad I(0) = I_0, \quad u(t) = \phi(t), \quad t \in [-\tau, 0]. \quad (15) $$

In our study, we will often link the results obtained for System (14) to the initial model (1), (2) and (7).

## 3 Steady-states and basic reproduction number

In this section, we discuss the existence of steady-states and establish the basic reproduction number associated with the model (14). We assume that the function $g$ is given by the expression (3) or that it satisfies the conditions (5). A steady-state $\bar{E} = (\bar{S}, \bar{I}, \bar{u})$ is a stationary solution of the system (14). It satisfies the system

$$ egin{cases} b - (a_S + \eta + g(\bar{I}))\bar{S} + \alpha\bar{I} + (1 - \nu)e^{-a_p \tau}\bar{u} = 0, \\ g(\bar{I})\bar{S} - (a_I + \alpha + \mu)\bar{I} = 0, \\ (1 - \nu e^{-a_p \tau})\bar{u} = \eta\bar{S}. \end{cases} \quad (16) $$

The third equation immediately yields

$$ \bar{u} = \frac{\eta}{1 - \nu e^{-a_p \tau}} \bar{S}. $$
Then, the system (16) becomes

$$ egin{cases} b - (a_S + \eta + g(\bar{I}))\bar{S} + \alpha\bar{I} + \frac{\eta(1-\nu)e^{-a_p\tau}}{1-\nu e^{-a_p\tau}}\bar{S} = 0, \\ g(\bar{I})\bar{S} - (a_I + \alpha + \mu)\bar{I} = 0. \end{cases} \quad (17) $$

We can remark that the corresponding steady-state for the system (1)-(2) is

$$ \bar{p}(a) = e^{-a_p a} \bar{u} = e^{-a_p a} \frac{\eta}{1 - \nu e^{-a_p \tau}} \bar{S}, \quad 0 < a < \tau, $$

with the total stationary protected population given by

$$ \bar{P} := \int_{0}^{\tau} \bar{p}(a) da = \frac{\eta(1 - e^{-a_p \tau})}{a_p(1 - \nu e^{-a_p \tau})} \bar{S}. $$

Since the parameter $\nu \in (0,1)$, then all the expressions established above are well defined.

Let $\bar{E}^0 := (\bar{S}^0, 0, \bar{u}^0)$ be the disease-free equilibrium of System (14). Then, thanks to the fact that $g(0) = 0$, we get, from the first equation of System (17), the expression

$$ \bar{S}^0 = \frac{b}{a_S + \eta \frac{1 - e^{-a_p \tau}}{1 - \nu e^{-a_p \tau}}} \quad \text{and} \quad \bar{u}^0 = \frac{\eta}{1 - \nu e^{-a_p \tau}} \bar{S}^0. \quad (18) $$

The disease-free equilibrium $\bar{E}^0$ always exists.

Let us now investigate the existence of an endemic equilibrium $E^* := (S^*, I^*, u^*)$, with $I^* > 0$. The first equation of System (17) gives

$$ \left[ a_S + \eta \frac{1 - e^{-a_p \tau}}{1 - \nu e^{-a_p \tau}} + g(I^*) \right] S^* = b + \alpha I^*. $$

In this last equation, we can identify the expression of $\bar{S}^0$, given by (18). That is

$$ \left[ \frac{b}{\bar{S}^0} + g(I^*) \right] S^* = b + \alpha I^*. $$

As $g(I^*) > 0$, we can deduce from the second equation of (17) that

$$ S^* = \frac{(a_I + \alpha + \mu)I^*}{g(I^*)}. $$

By substituting $S^*$ in the above equation, we obtain

$$ \left[ \frac{b}{\bar{S}^0} + g(I^*) \right] \frac{(a_I + \alpha + \mu)I^*}{g(I^*)} = b + \alpha I^*. $$
We conclude that $I^* > 0$ is a solution of the equation

$$g(I^*) = \frac{b(a_I + \alpha + \mu)I^*}{\bar{S}^0(b - (a_I + \mu)I^*)}, \qquad (19)$$

with $\bar{S}^0$ given by (18). The existence of an endemic steady-state depends only on the existence of $I^* > 0$ solution of Equation (19). It is also interesting to note that the last expression involves the disease-free equilibrium $\bar{S}^0$.

We can see from (19) that in fact $0 < I^* < \frac{b}{a_I + \mu}$. Then, we define the function $f$ by

$$f(I^*) := \frac{b(a_I + \alpha + \mu)I^*}{\bar{S}^0(b - (a_I + \mu)I^*)}, \quad 0 < I^* < \frac{b}{a_I + \mu}.$$

We can notice that $f$ is convex and $g$ is concave, and both are increasing, with $g(0) = f(0) = 0$. Then, a necessary and sufficient condition for having the existence of a unique $I^* > 0$ such that $g(I^*) = f(I^*)$ is $g'(0) > f'(0)$ (see Figure 3 for an illustration). The derivative of the function $f$ is given by

$$f'(I^*) = \frac{b^2(a_I + \alpha + \mu)}{\bar{S}^0(b - (a_I + \mu)I^*)^2}, \quad 0 < I^* < \frac{b}{a_I + \mu}.$$

So,

$$f'(0) = \frac{a_I + \alpha + \mu}{\bar{S}^0}.$$

Therefore, the condition $g'(0) > f'(0)$ is equivalent to

$$\frac{g'(0)\bar{S}^0}{a_I + \alpha + \mu} > 1. \qquad (20)$$

We put

$$\mathcal{R}_0 := \frac{g'(0)\bar{S}^0}{a_I + \alpha + \mu}. \qquad (21)$$

$\mathcal{R}_0$ is the basic reproduction number associated with our model. Observe that in the expression of $\mathcal{R}_0$, the term $\bar{S}^0$ represents the average number of susceptible individuals in a population where there are no infected individuals and $g'(0)/(a_I + \alpha + \mu)$ is the rate of these susceptible individuals that get infected. It corresponds to the average number of newly infected susceptible individuals generated by an infected individual in a population where all individuals are susceptible. According to this threshold, we reformulate the existence theorem of the steady-states.

**Theorem 3.** (i) The hybrid system (14) always admits a disease-free equilibrium $\bar{E}^0 := (\bar{S}^0, 0, \bar{u}^0)$ with

$$\bar{S}^0 = \frac{b}{a_S + \eta \frac{1 - e^{-a_p \tau}}{1 - \nu e^{-a_p \tau}}} \quad \text{and} \quad \bar{u}^0 = \frac{\eta}{1 - \nu e^{-a_p \tau}} \bar{S}^0.$$
<figure><img src="image_5.png" /><figcaption>Fig. 3: Curves of functions f and g for $\mathcal{R}_0 > 1$. The abscissa of the point of intersection of the two curves corresponds to the $I^*$ of the endemic equilibrium.</figcaption></figure>

(ii) A unique endemic equilibrium $E^* = (S^*, I^*, u^*)$ of the hybrid system (14) exists if and only if $\mathcal{R}_0 > 1$, with $\mathcal{R}_0$ given by (21), $I^*$ is the unique positive solution of (19), and

$$S^* = \frac{(a_I + \alpha + \mu)I^*}{g(I^*)}, \quad u^* = \frac{\eta}{1 - \nu e^{-a_p \tau}} S^*.$$

For the case $g(I) = \beta I$, we have explicitly

$$S^* = \frac{a_I + \alpha + \mu}{\beta}, \quad u^* = \frac{\eta}{1 - \nu e^{-a_p \tau}} S^* \quad \text{and} \quad I^* = \frac{b(a_I + \alpha + \mu)}{\beta \bar{S}^0(a_I + \mu)} [\mathcal{R}_0 - 1].$$

As indicated above, the two steady-states of the initial system (1)-(2) can be obtained by using the expression

$$\bar{p}(a) = e^{-a_p a} \bar{u} = e^{-a_p a} \frac{\eta}{1 - \nu e^{-a_p \tau}} \bar{S}, \quad 0 < a < \tau.$$

The total stationary protected population is also given by

$$\bar{P} := \int_{0}^{\tau} \bar{p}(a) da = \frac{\eta(1 - e^{-a_p \tau})}{a_p(1 - \nu e^{-a_p \tau})} \bar{S}.$$

## 4 Local asymptotic stability

In this section, we are interested in the local asymptotic stability of the two steady-states of the hybrid system (14). We assume that the function g is given by the expression (3) or that it satisfies the conditions (5). We are now going to linearize the system (14) around the two steady-states. Since the difference equation (the equation of u) is linear, this can be achieved directly, as with classical systems of delay differential equations (see the book of J.K. Hale and S.M. Verduyn Lunel [28]). In fact, using
the same notations as in [28], the difference equation has the following form

$$D(u_t) = \eta S(t), \qquad (22)$$

where $D: C([-\tau, 0]) \to \mathbb{R}$ is the linear and bounded operator defined by

$$D(\phi) = \phi(0) - \nu e^{-a_p \tau} \phi(-\tau), \quad \phi \in C([-\tau, 0]),$$

and $u_t, t \ge 0$, is the element of $C([-\tau, 0])$ defined for any $\theta \in [-\tau, 0]$ by the translation $u_t(\theta) := u(t + \theta)$. We refer the reader to [28], pages 93 and 274, for more details on the properties of the equation (22). In particular, since $\nu \in (0, 1)$, the operator $D$ is stable (see Definition 3.1 of [28]).

We linearize the system (14) around any steady-state $\bar{E} = (\bar{S}, \bar{I}, \bar{u})$. This gives the following linear system

$$\left\{
\begin{aligned}
s'(t) &= -(a_S + \eta + g(\bar{I}))s(t) - g'(\bar{I})\bar{S}i(t) + \alpha i(t) + (1-\nu)e^{-a_p \tau}v(t-\tau), \\
i'(t) &= g(\bar{I})s(t) + g'(\bar{I})\bar{S}i(t) - (a_I + \alpha + \mu)i(t), \\
v(t) &= \eta s(t) + \nu e^{-a_p \tau}v(t-\tau).
\end{aligned}
\right.
\qquad (23)$$

To establish the characteristic equation, we look for exponential solutions of (23),

$$s(t) = e^{\lambda t} s_0, \quad i(t) = e^{\lambda t} i_0, \quad v(t) = e^{\lambda t} v_0, \quad \lambda \in \mathbb{C}.$$

This is equivalent to finding $\lambda \in \mathbb{C}$ and a non-trivial vector $(s_0, i_0, v_0)$ solution of the system

$$\left\{
\begin{array}{l}
(\lambda + \eta + a_S + g(\bar{I}))s_0 + (g'(\bar{I})\bar{S} - \alpha)i_0 - (1-\nu)e^{-(\lambda+a_p)\tau}v_0 = 0, \\
-\eta s_0 + (1-\nu e^{-(\lambda+a_p)\tau})v_0 = 0, \\
-g(\bar{I})s_0 + (\lambda + a_I + \alpha + \mu - g'(\bar{I})\bar{S})i_0 = 0.
\end{array}
\right.
$$

Then, we obtain the characteristic equation

$$\Delta(\lambda) := \begin{vmatrix} \lambda + A & g'(\bar{I})\bar{S} - \alpha & -(1-\nu)e^{-(\lambda+a_p)\tau} \\ -\eta & 0 & 1-\nu e^{-(\lambda+a_p)\tau} \\ -g(\bar{I}) & \lambda + B & 0 \end{vmatrix} = 0,$$

with

$$A := \eta + a_S + g(\bar{I}) \quad \text{and} \quad B := a_I + \alpha + \mu - g'(\bar{I})\bar{S}.$$

After some development, we get a transcendental equation with delay dependent coefficients

$$\Delta(\lambda) = -\eta(1-\nu)e^{-a_p \tau}e^{-\lambda \tau}(\lambda - g'(\bar{I})\bar{S} + a_I + \alpha + \mu) + (1-\nu e^{-a_p \tau}e^{-\lambda \tau}) \times \\ [(\lambda + \eta + a_S + g(\bar{I}))(\lambda - g'(\bar{I})\bar{S} + a_I + \alpha + \mu) + g(\bar{I})(g'(\bar{I})\bar{S} - \alpha)].$$
This equation has the form

$$ \Delta(\lambda) = \lambda^2 + p_1(\tau)\lambda + p_0(\tau) + [q_2(\tau)\lambda^2 + q_1(\tau)\lambda + q_0(\tau)]e^{-\lambda\tau} = 0, $$

with $q_2(\tau) = \nu e^{-a_p\tau} > 0$. E. Beretta and Y. Kuang [11], considered the following general transcendental equation

$$ \sum_{k=0}^{n} p_k(\tau)\lambda^k + \left[ \sum_{k=0}^{m} q_k(\tau)\lambda^k \right] e^{-\lambda\tau} = 0, \quad p_n(\tau), q_m(\tau) \neq 0. $$

However, their approach is only applicable where $n > m$. In our case, we have $n = m = 2$.

We first consider the case of a disease-free equilibrium $\bar{E}^0 := (\bar{S}^0, 0, \bar{u}^0)$.

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

The linearized system (23) about the disease-free steady-state $\bar{E}^0 := (\bar{S}^0, 0, \bar{u}^0)$ is,

$$ \left\{ \begin{aligned} s'(t) &= -(a_S + \eta)s(t) - g'(0)\bar{S}^0 i(t) + \alpha i(t) + (1 - \nu)e^{-a_p\tau}v(t - \tau), \\ i'(t) &= g'(0)\bar{S}^0 i(t) - (a_I + \alpha + \mu)i(t), \\ v(t) &= \eta s(t) + \nu e^{-a_p\tau}v(t - \tau). \end{aligned} \right. $$

The characteristic equation $\Delta(\lambda) = 0$ becomes

$$ \begin{aligned} & -\eta(1 - \nu)e^{-a_p\tau}e^{-\lambda\tau}(\lambda - g'(0)\bar{S}^0 + a_I + \alpha + \mu) \\ & + (1 - \nu e^{-a_p\tau}e^{-\lambda\tau})(\lambda + \eta + a_S)(\lambda - g'(0)\bar{S}^0 + a_I + \alpha + \mu) = 0. \end{aligned} \quad (24) $$

We can see that

$$ \lambda_0 = g'(0)\bar{S}^0 - a_I - \alpha - \mu = (a_I + \alpha + \mu)(\mathcal{R}_0 - 1). $$

is a real eigenvalue. The other eigenvalues are solutions of the transcendental equation

$$ \Delta_0(\lambda) := -\eta(1 - \nu)e^{-a_p\tau}e^{-\lambda\tau} + (1 - \nu e^{-a_p\tau}e^{-\lambda\tau})(\lambda + \eta + a_S) = 0. \quad (25) $$

**Theorem 4.** (i) Assume that $\mathcal{R}_0 > 1$. Then, the real root $\lambda_0$ is positive, and the steady-state $\bar{E}^0 := (\bar{S}^0, 0, \bar{u}^0)$ is unstable.

(ii) Assume that $\mathcal{R}_0 < 1$. Then, the real root $\lambda_0$ is negative and all roots of the equation (25) have negative real parts. This means that the steady-state $\bar{E}^0 := (\bar{S}^0, 0, \bar{u}^0)$ is locally asymptotically stable.

*Proof.* We have the real eigenvalue

$$ \lambda_0 = (a_I + \alpha + \mu)(\mathcal{R}_0 - 1). $$
(i) If $\mathcal{R}_0 > 1$, $\lambda_0 > 0$. Then, the steady state $\bar{E}^0 := (\bar{S}^0, 0, \bar{u}^0)$ is unstable.

(ii) Suppose that $\mathcal{R}_0 < 1$. Then, $\lambda_0 < 0$. Now, we show that there is no root with nonnegative real part of the equation (25). Indeed, if we suppose that $\lambda = x + i\omega$, with $x \ge 0$, is a root of the equation (25), we obtain

$$e^{-(x+a_p+i\omega)\tau} = \frac{x+\eta+a_S+i\omega}{\eta+\nu(x+a_S)+i\nu\omega}.$$

By squaring the modulus of this last equality, we obtain

$$\frac{(x + \eta + a_S)^2 + \omega^2}{(\eta + \nu(x + a_S))^2 + \nu^2\omega^2} = |e^{-(x+a_p+i\omega)\tau}|^2 = e^{-2(x+a_p)\tau} < 1.$$

Then, we get

$$(1 - \nu^2)\omega^2 + (1 - \nu)(x + a_S)[(1 + \nu)(x + a_S) + 2\eta] < 0.$$

This is absurd. Hence, no root with positive real part exists. □

We now consider the case of the endemic equilibrium $E^* = (S^*, I^*, u^*)$ of the system (14).

## 4.2 Local asymptotic stability of the endemic equilibrium

The endemic steady-state $E^* = (S^*, I^*, u^*)$ of the system (14) exists if and only if $\mathcal{R}_0 > 1$. The associate characteristic equation is given by

$$-\eta(1 - \nu)e^{-a_p\tau}e^{-\lambda\tau}(\lambda - g'(I^*)S^* + a_I + \alpha + \mu) + (1 - \nu e^{-a_p\tau}e^{-\lambda\tau}) \times \quad (26) \\ [(\lambda + \eta + a_S + g(I^*)) (\lambda - g'(I^*)S^* + a_I + \alpha + \mu) + g(I^*)(g'(I^*)S^* - \alpha)] = 0.$$

The following result establishes the local asymptotic stability of $E^*$.

**Theorem 5.** Suppose that $\mathcal{R}_0 > 1$. Then, all roots of the characteristic equation (26) have negative real parts, and the endemic steady-state $E^*$ is locally asymptotically stable.

*Proof.* The characteristic equation becomes

$$\Delta(\lambda) := [(\lambda + \eta + a_S + g(I^*)) (\lambda - g'(I^*)S^* + a_I + \alpha + \mu) + g(I^*)(g'(I^*)S^* - \alpha)] \\ -\nu e^{-a_p\tau} (\lambda + a_S + g(I^*)) (\lambda - g'(I^*)S^* + a_I + \alpha + \mu) e^{-\lambda\tau} \\ -\nu e^{-a_p\tau} g(I^*)(g'(I^*)S^* - \alpha) e^{-\lambda\tau} - \eta e^{-a_p\tau} (\lambda - g'(I^*)S^* + a_I + \alpha + \mu) e^{-\lambda\tau}, \\ = 0.$$

It can be rewritten as follows,

$$\lambda^2 + (A + \eta)\lambda + E + \eta C = [\nu(\lambda^2 + A\lambda + E) + \eta C + \eta\lambda] e^{-a_p\tau}e^{-\lambda\tau}, \quad (27)$$
with

$$A = -g'(I^*)S^* + a_I + \alpha + \mu + a_S + g(I^*),$$

$$E = (-g'(I^*)S^* + a_I + \alpha + \mu)(a_S + g(I^*)) + g(I^*)(g'(I^*)S^* - \alpha),$$

$$C = -g'(I^*)S^* + a_I + \alpha + \mu.$$

It is not difficult to see that

$$A = C + a_S + g(I^*),$$

$$E = Ca_S + (a_I + \mu)g(I^*).$$

Furthermore, using the second equation of (17), we obtain

$$C := - \left[ g'(I^*) - \frac{g(I^*)}{I^*} \right] S^*.$$

* If $g$ is given by (3), then $C = 0$.

* If $g$ satisfies the conditions (5), then using Lemma 1, we get $C \ge 0$.

In all cases, we have shown that

$$A > 0, \quad E > 0 \quad \text{and} \quad C \ge 0.$$

Let us suppose that there exists a root $\lambda = x + iy$ of the characteristic equation such that $x \ge 0$. In that case, we have

$$|e^{-a_p \tau} e^{-\lambda \tau}| = e^{-a_p \tau} e^{-x \tau} \le 1.$$

Then, by taking the modulus of (27), we obtain

$$| \lambda^2 + (A + \eta)\lambda + E + \eta C | \le | \nu(\lambda^2 + A\lambda + E) + \eta C + \eta \lambda |.$$

By taking square in both sides, we get

$$[\lambda^2 + (A + \eta)\lambda + E + \eta C] [\bar{\lambda}^2 + (A + \eta)\bar{\lambda} + E + \eta C] \\ \le [\nu(\lambda^2 + A\lambda + E) + \eta C + \eta \lambda] [\nu(\bar{\lambda}^2 + A\bar{\lambda} + E) + \eta C + \eta \bar{\lambda}].$$

After arranging this inequality and knowing that $0 < \nu < 1$, we obtain

$$(1 + \nu) | \lambda^2 + A\lambda + E |^2 + [(\lambda^2 + A\lambda + E)(\eta C + \eta \bar{\lambda}) + (\bar{\lambda}^2 + A\bar{\lambda} + E)(\eta C + \eta \lambda)] \le 0.$$

It is clear that the first term of the last inequality is positive. For the second term, we set

$$\Gamma := (\lambda^2 + A\lambda + E)(\eta C + \eta \bar{\lambda}) + (\bar{\lambda}^2 + A\bar{\lambda} + E)(\eta C + \eta \lambda).$$
This can be rewritten as follows

$$ \Gamma = 2\eta C(x^2 - y^2) + 2\eta x(x^2 + y^2) + 2\eta EC + 2E\eta x + 2\eta ACx + 2\eta A(x^2 + y^2). $$

Then, by replacing A by $C + a_S + g(I^*)$ in the last term, we obtain

$$ \Gamma = 4\eta C x^2 + 2\eta(a_S + g(I^*)) (x^2 + y^2) + \eta x |\lambda|^2 + 2\eta EC + 2E\eta x + 2\eta ACx. $$

Since $x \ge 0$ and $C \ge 0$ then $\Gamma \ge 0$. We conclude that

$$ (1 + \nu) |\lambda^2 + A\lambda + E|^2 + \Gamma > 0. $$

This is a contradiction. Hence, when $\mathcal{R}_0 > 1$ only roots with negative real part exist. So, the endemic steady-state is locally asymptotically stable. $\square$

# 5 Global asymptotic stability

In this section, we study the global asymptotic stability of the two steady-states of the hybrid system (14). We consider the SIR and SIS scenarios, as well as different cases of infection force: mass action (3) or a nonlinear function $g$ satisfying the conditions (5). Specifically, we examine the global asymptotic stability of the two steady-states in the two scenarios:

* The force of infection is $g(I) = \beta I$ (mass action) and $\alpha > 0$.

* The force of infection satisfies the conditions (5) and $\alpha = 0$.

## 5.1 Global asymptotic stability of the disease-free steady-state

**Theorem 6.** Suppose that $g(I) = \beta I$, for all $I \ge 0$, $\alpha > 0$ and $\mathcal{R}_0 < 1$. Then, the disease-free steady-state $\bar{E}^0 := (\bar{S}^0, 0, \bar{u}^0)$ of the hybrid system (14) is globally asymptotically stable.

For the proof of this theorem, we use the following version of Barbalat's lemma ([32], Lemma 8.2, page 323).

**Lemma 7.** Let $f: \mathbb{R} \to \mathbb{R}$ be a uniformly continuous function on $[0, +\infty)$. Suppose that $\lim_{t \to +\infty} \int_0^t f(s)ds$ exists and is finite. Then, $\lim_{t \to +\infty} f(t) = 0$.

We assume that $\mathcal{R}_0 < 1$. Then, the disease-free steady-state is the unique equilibrium of (14). Recall that the disease-free steady-state of (14) is given by $\bar{E}^0 := (\bar{S}^0, 0, \bar{u}^0)$ with

$$ \bar{S}^0 = \frac{b}{a_S + \eta \frac{1 - e^{-a_p \tau}}{1 - \nu e^{-a_p \tau}}} \quad \text{and} \quad \bar{u}^0 = \frac{\eta}{1 - \nu e^{-a_p \tau}} \bar{S}^0. $$

Consider the case $g(I) = \beta I$ and $\alpha > 0$. For any arbitrary solution $(S(t), I(t), u(t))$ of System (14), we put, for $t \ge 0$,

$$ \hat{S}(t) = S(t) - \bar{S}^0 \quad \text{and} \quad \hat{u}(t) = u(t) - \bar{u}^0. $$
The first equation of the hybrid system (14) becomes

$$ \hat{S}'(t) = -(\eta + a_S + \beta I(t)) \hat{S}(t) - \beta \bar{S}^0 I(t) + \alpha I(t) + (1 - \nu)e^{-a_p \tau} \hat{u}(t - \tau). \quad (28) $$

The last two equations of (14) are

$$ \left\{ \begin{aligned} I'(t) &= \beta I(t) \hat{S}(t) + \beta \bar{S}^0 I(t) - (a_I + \alpha + \mu) I(t), \\ \hat{u}(t) &= \eta \hat{S}(t) + \nu e^{-a_p \tau} \hat{u}(t - \tau). \end{aligned} \right. \quad (29) $$

We can now prove Theorem 6.

*Proof.* First of all, note that the function $I(t)$ is non-negative, whereas $\hat{S}(t)$ and $\hat{u}(t)$ are not necessarily so. We therefore need to square these latter functions in our construction of a Lyapunov function. Consider the following non-negative function, defined along the solutions of the system (28)-(29),

$$ V(t) = \frac{\delta_1}{2} \hat{S}^2(t) + \int_0^\tau \hat{u}^2(t-s)ds + \delta_2 I(t), $$

with $\delta_1$ and $\delta_2$ two positive parameters to find. Indeed, as $I(t)$, $\hat{S}^2(t)$ and $\hat{u}^2(t)$ are non-negative, then $V(t)$ is non-negative. The function $V(t)$ can be written as follows

$$ V(t) = \frac{\delta_1}{2} \hat{S}(t)^2 + \int_{t-\tau}^{t} \hat{u}^2(s)ds + \delta_2 I(t). $$

So, we can calculate the derivative

$$ V'(t) = \delta_1 \hat{S}'(t) \hat{S}(t) + \hat{u}^2(t) - \hat{u}^2(t - \tau) + \delta_2 I'(t). $$

The first objective is to find $\delta_1 > 0$ and $\delta_2 > 0$ such that $V'(t) < 0$. Using the expression of the solutions of the system (28)-(29), we obtain, for $t > 0$,

$$ \begin{aligned} V'(t) = & -\delta_1(\eta + a_S)\hat{S}^2(t) - \delta_1\beta\hat{S}^2(t)I(t) - \delta_1\beta\bar{S}^0\hat{S}(t)I(t) + \delta_1\alpha\hat{S}(t)I(t) \\ & +\delta_1(1 - \nu)\hat{S}(t)e^{-a_p\tau}\hat{u}(t - \tau) \\ & +\eta^2\hat{S}^2(t) + \nu^2e^{-2a_p\tau}\hat{u}^2(t - \tau) + 2\eta\nu e^{-a_p\tau}\hat{S}(t)\hat{u}(t - \tau) - \hat{u}^2(t - \tau) \\ & +\delta_2\beta\hat{S}(t)I(t) + \delta_2\beta\bar{S}^0I(t) - \delta_2(a_I + \alpha + \mu)I(t). \end{aligned} $$

By grouping some terms, we get

$$ \begin{aligned} V'(t) = & -\delta_1\beta\hat{S}^2(t)I(t) - (\delta_1(\eta + a_S) - \eta^2)\hat{S}^2(t) \\ & + [\delta_1(1 - \nu) + 2\eta\nu]e^{-a_p\tau}\hat{u}(t - \tau)\hat{S}(t) \\ & -(e^{2a_p\tau} - \nu^2)e^{-2a_p\tau}\hat{u}^2(t - \tau) + (\delta_1\alpha + \delta_2\beta - \delta_1\beta\bar{S}^0)\hat{S}(t)I(t) \\ & -\delta_2(a_I + \alpha + \mu)(1 - \mathcal{R}_0)I(t). \end{aligned} $$
We choose

$$a := \delta_1(\eta + a_S) - \eta^2,$$

$$b := \delta_1(1 - \nu) + 2\eta\nu,$$

$$c := e^{2a_p\tau} - \nu^2.$$

As $0 < \nu < 1$, it is clear that parameters $b$ and $c$ are always positive. For the parameter $a$ to be positive, we have to choose

$$\delta_1 > \frac{\eta^2}{\eta + a_S}. \qquad (30)$$

For the rest of the proof, there are two possible scenarios to consider: $\beta\bar{S}^0 - \alpha > 0$ and $\beta\bar{S}^0 - \alpha \le 0$.

Suppose $\beta\bar{S}^0 - \alpha > 0$. Then, we choose $\delta_2 = \delta_1 \left(\bar{S}^0 - \frac{\alpha}{\beta}\right) > 0$. Using the notations $a, b$ and $c$, we get

$$ \begin{aligned} V'(t) &= -a\hat{S}^2(t) + be^{-a_p\tau}\hat{u}(t-\tau)\hat{S}(t) - ce^{-2a_p\tau}\hat{u}^2(t-\tau) \\ &\quad -\delta_1\beta\hat{S}^2(t)I(t) - \delta_2(a_I + \alpha + \mu)(1 - \mathcal{R}_0)I(t), \\ &\le -a\hat{S}^2(t) + be^{-a_p\tau}\hat{u}(t-\tau)\hat{S}(t) - ce^{-2a_p\tau}\hat{u}^2(t-\tau) \\ &\quad -\delta_2(a_I + \alpha + \mu)(1 - \mathcal{R}_0)I(t). \end{aligned} \qquad (31) $$

In other words

$$ V'(t) \le -c \left[ \left( e^{-a_p \tau} \hat{u}(t-\tau) - \frac{b}{2c} \hat{S}(t) \right)^2 - \frac{b^2 - 4ac}{4c^2} \hat{S}^2(t) \right] \\ - \delta_2 (a_I + \alpha + \mu) (1 - \mathcal{R}_0) I(t). \qquad (32) $$

Now we have to determine the sign of the expression

$$\Delta(\delta_1) := b^2 - 4ac.$$

We would like to find $\delta_1 > 0$ such that $\Delta(\delta_1) < 0$. We have,

$$ \begin{aligned} \Delta(\delta_1) &= (\delta_1(1-\nu) + 2\eta\nu)^2 - 4(-\delta_1(\eta + a_S) + \eta^2)(\nu^2 - e^{2a_p\tau}), \\ &= \delta_1^2(1-\nu)^2 + 4\delta_1(1-\nu)\eta\nu + 4\eta^2\nu^2 \\ &\quad -4(-\delta_1(\eta + a_S)(\nu^2 - e^{2a_p\tau}) + \eta^2(\nu^2 - e^{2a_p\tau})), \\ &= (1-\nu)^2\delta_1^2 + 4[(1-\nu)\eta\nu + (\nu^2 - e^{2a_p\tau})(\eta + a_S)]\delta_1 + 4\eta^2e^{2a_p\tau}, \\ &= (1-\nu)^2\delta_1^2 + 4[\eta\nu + a_S\nu^2 - e^{2a_p\tau}(\eta + a_S)]\delta_1 + 4\eta^2e^{2a_p\tau}. \end{aligned} $$

The function $\Delta$ is a second-degree polynomial of $\delta_1$. Its derivative is

$$\Delta'(\delta_1) = 2(1 - \nu)^2 \delta_1 + 4[\eta\nu + a_S\nu^2 - e^{2a_p\tau}(\eta + a_S)].$$
The minimum of the function $\Delta$ is then reached for $\bar{\delta}_1 > 0$ given by

$$ \bar{\delta}_1 := 2 \frac{\eta(e^{2a_p\tau} - \nu) + a_S(e^{2a_p\tau} - \nu^2)}{(1 - \nu)^2} > 0. \quad (33) $$

This value $\bar{\delta}_1$ satisfies the condition (30). Indeed,

$$ \begin{aligned} a &= \bar{\delta}_1(\eta + a_S) - \eta^2, \\ &= 2 \frac{\eta(e^{2a_p\tau} - \nu) + a_S(e^{2a_p\tau} - \nu^2)}{(1 - \nu)^2} (\eta + a_S) - \eta^2, \\ &= 2 \frac{\eta^2(e^{2a_p\tau} - \nu) - \frac{1}{2}(1 - \nu)^2\eta^2 + \eta a_S(e^{2a_p\tau} - \nu) + a_S(\eta + a_S)(e^{2a_p\tau} - \nu^2)}{(1 - \nu)^2}, \\ &= \frac{2}{(1 - \nu)^2} \left[ \eta^2 \left( e^{2a_p\tau} - \frac{1}{2}(1 + \nu^2) \right) + \eta a_S(e^{2a_p\tau} - \nu) + a_S(\eta + a_S)(e^{2a_p\tau} - \nu^2) \right]. \end{aligned} $$

Since $0 < \nu < 1$, we have $0 < \frac{1}{2}(1 + \nu^2) < 1$. Then,

$$ \eta^2 \left( e^{2a_p\tau} - \frac{1}{2}(1 + \nu^2) \right) + \eta a_S(e^{2a_p\tau} - \nu) + a_S(\eta + a_S)(e^{2a_p\tau} - \nu^2) > 0. $$

We still have to show that $\Delta(\bar{\delta}_1) < 0$. Let's return to the second-degree polynomial $\delta_1 \mapsto \Delta(\delta_1)$ and note that its discriminant is given by

$$ \begin{aligned} \bar{\Delta} &= 16 \left[ (\eta\nu + a_S\nu^2 - e^{2a_p\tau}(\eta + a_S))^2 - \eta^2(1 - \nu)^2 e^{2a_p\tau} \right], \\ &= 16 \left[ (\eta\nu + a_S\nu^2 - e^{2a_p\tau}(\eta + a_S)) - \eta(1 - \nu)e^{a_p\tau} \right] \\ & \quad \times \left[ (\eta\nu + a_S\nu^2 - e^{2a_p\tau}(\eta + a_S)) + \eta(1 - \nu)e^{a_p\tau} \right], \\ &= 16 \left[ \eta(1 + e^{a_p\tau})(\nu - e^{a_p\tau}) + a_S(\nu^2 - e^{2a_p\tau}) \right] \\ & \quad \times \left[ \eta(1 - e^{a_p\tau})(\nu + e^{a_p\tau}) + a_S(\nu^2 - e^{2a_p\tau}) \right]. \end{aligned} $$

Since $0 < \nu < 1$, we have

$$ \left\{ \begin{array}{l} \eta(1 + e^{a_p\tau})(\nu - e^{a_p\tau}) + a_S(\nu^2 - e^{2a_p\tau}) < 0, \\ \eta(1 - e^{a_p\tau})(\nu + e^{a_p\tau}) + a_S(\nu^2 - e^{2a_p\tau}) < 0. \end{array} \right. $$

So, $\bar{\Delta} > 0$. Then, for the value $\delta_1 = \bar{\delta}_1$ given by (33), we have

$$ \Delta(\bar{\delta}_1) := b^2 - 4ac < 0. $$

We then conclude from the inequality (32), that

$$ V'(t) < 0. $$
Now, we consider the case $\beta \bar{S}^0 - \alpha \le 0$. Remember, we have

$$ 
\begin{aligned}
V'(t) = & -a\hat{S}^2(t) + be^{-a_p\tau}\hat{u}(t-\tau)\hat{S}(t) - ce^{-2a_p\tau}\hat{u}^2(t-\tau) \quad (34) \\
& -\delta_1\beta\hat{S}^2(t)I(t) + (\delta_1\alpha + \delta_2\beta - \delta_1\beta\bar{S}^0)\hat{S}(t)I(t) - \delta_2(a_I + \alpha + \mu)(1 - \mathcal{R}_0)I(t), \\
= & -a\hat{S}^2(t) + be^{-a_p\tau}\hat{u}(t-\tau)\hat{S}(t) - ce^{-2a_p\tau}\hat{u}^2(t-\tau) \\
& + (-\delta_1\beta\hat{S}^2(t) + [\delta_1(\alpha - \beta\bar{S}^0) + \delta_2\beta]\hat{S}(t) - \delta_2(a_I + \alpha + \mu)(1 - \mathcal{R}_0))I(t).
\end{aligned}
$$

As the parameters $a, b$ and $c$ are independent on $\delta_2$ and with the choice of $\delta_1 = \bar{\delta}_1$ given by (33), we already have

$$ -a\hat{S}^2(t) + be^{-a_p\tau}\hat{u}(t-\tau)\hat{S}(t) - ce^{-2a_p\tau}\hat{u}^2(t-\tau) < 0. $$

It remains to find $\delta_2 > 0$ such that

$$ \hat{S}(t) \mapsto -\delta_1\beta\hat{S}^2(t) + [\delta_1(\alpha - \beta\bar{S}^0) + \delta_2\beta]\hat{S}(t) - \delta_2(a_I + \alpha + \mu)(1 - \mathcal{R}_0) < 0. \quad (35) $$

Clearly, $-\delta_2(a_I + \alpha + \mu)(1 - \mathcal{R}_0) < 0$ and $\delta_1(\alpha - \beta\bar{S}^0) + \delta_2\beta > 0$. We now compute the discriminant of the second-degree polynomial (35),

$$ 
\begin{aligned}
\Delta(\delta_2) &= [\delta_1(\alpha - \beta\bar{S}^0) + \delta_2\beta]^2 - 4\delta_1\delta_2\beta(a_I + \alpha + \mu)(1 - \mathcal{R}_0), \\
&= \beta^2\delta_2^2 + 2\delta_1\beta[\alpha - \beta\bar{S}^0 - 2(a_I + \alpha + \mu)(1 - \mathcal{R}_0)]\delta_2 + \delta_1^2(\alpha - \beta\bar{S}^0)^2, \\
&= \beta^2\delta_2^2 + 2\delta_1\beta[\alpha - \beta\bar{S}^0 - 2(a_I + \alpha + \mu - \beta\bar{S}^0)]\delta_2 + \delta_1^2(\alpha - \beta\bar{S}^0)^2, \\
&= \beta^2\delta_2^2 - 2\delta_1\beta[\alpha - \beta\bar{S}^0 + 2(a_I + \mu)]\delta_2 + \delta_1^2(\alpha - \beta\bar{S}^0)^2.
\end{aligned}
$$

We are looking for $\bar{\delta}_2$ such that $\Delta(\bar{\delta}_2) < 0$, to satisfy (35). We compute the discriminant of $\Delta(\delta_2)$ as a second-degree polynomial of $\delta_2$ and we get

$$ 
\begin{aligned}
\bar{\Delta} &= 4\delta_1^2\beta^2[(\beta\bar{S}^0 - \alpha - 2(a_I + \mu))^2 - (\alpha - \beta\bar{S}^0)^2], \\
&= 16\delta_1^2\beta^2(a_I + \mu)(a_I + \mu + \alpha)(1 - \mathcal{R}_0), \\
&> 0.
\end{aligned}
$$

With any choice of $\delta_1 > 0$ and in particular for $\delta_1 = \bar{\delta}_1$, we have $\bar{\Delta} > 0$. This means that $\Delta(\delta_2) = 0$ has two positive real roots. Indeed, the expression of $\Delta(\delta_2)$ tells us that the product and the sum of these two roots are both positive, so the roots are positive. If we choose any value of $\delta_2$ between these positive real roots, we obtain, for $\delta_2 = \bar{\delta}_2$, that $\Delta(\delta_2) < 0$. Therefore, we can conclude that

$$ -\delta_1\beta\hat{S}^2(t) + [\delta_1(\alpha - \beta\bar{S}^0) + \delta_2\beta]\hat{S}(t) - \delta_2(a_I + \alpha + \mu)(1 - \mathcal{R}_0) < 0. $$

Once again, we proved that $V'(t) < 0$. $V$ is therefore a non-increasing function. Since $V$ is also non-negative, we have

$$ \lim_{t \to +\infty} V(t) = \inf_{s \ge 0} V(s) := V^* \ge 0. $$
Furthermore, according to (32) and (34), we have

$$V'(t) \leq -\alpha_1 \hat{S}^2(t), \qquad (36)$$

with

$$\alpha_1 := \frac{4ac - b^2}{4c} > 0.$$

By integrating both sides of the inequality (36), we obtain

$$\alpha_1 \int_0^t \hat{S}^2(s) ds \leq V(0) - V(t). \qquad (37)$$

All the functions involved in the inequality (37) are non-decreasing and bounded. Then their limits, when $t \to +\infty$, exist and we have

$$\lim_{t \to +\infty} \int_0^t \hat{S}^2(s) ds \leq \frac{1}{\alpha_1} [V(0) - V^*].$$

Using Theorem 2, we know that the function $\hat{S}(t)$ is bounded. So, directly from the equation (28), we see that the derivative is also bounded. Then, this function is uniformly continuous. Applying Lemma 7 to the function $t \mapsto \hat{S}^2(t)$, we deduce that,

$$\lim_{t \to +\infty} \hat{S}^2(t) = 0,$$

and so,

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

We will now prove that

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

We know that $\hat{u}$ satisfies the inhomogeneous difference equation

$$D(\hat{u}_t) = \eta \hat{S}(t),$$

with

$$D(\phi) := \phi(0) - \nu e^{-a_p \tau} \phi(-\tau), \quad \phi \in C([-\tau, 0]) \quad \text{and} \quad \lim_{t \to +\infty} \hat{S}(t) = 0.$$

As the linear operator $D$ is stable, because $0 < \nu e^{-a_p \tau} < 1$, we can apply Lemma 3.5 of [20], to conclude that

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

As $\lim_{t \to +\infty} V'(t) = 0$, we conclude from (31) and (34) that

$$\lim_{t \to +\infty} I(t) = 0.$$
So, we have the global asymptotic stability of the disease-free equilibrium. □

We are now interested in the global asymptotic stability of the disease-free equilib-
rium in the case where the function *g* satisfies the conditions (5), but here we assume
that *α* = 0.

**Theorem 8.** Suppose that the function *g* satisfies the conditions (5), *α* = 0 and *R*₀ < 1. Then, the disease-free steady-state *E*⁰ := (*S*⁰, 0, *u*⁰) of the hybrid system (14) is globally asymptotically stable.

*Proof.* The proof is inspired by that of Theorem 6. Consider again the Lyapunov function

$$
V(t) = \frac{\delta_1}{2} \hat{S}^2(t) + \int_{t-\tau}^{t} \hat{u}^2(s)ds + \delta_2 I(t),
$$

where

$$
\delta_1 = 2 \frac{\eta(e^{2a_p \tau} - \nu) + a_S(e^{2a_p \tau} - \nu^2)}{(1-\nu)^2} > 0. \quad (38)
$$

Suppose for the moment, that $\alpha \ge 0$. Then, the derivative of $V$ along the solutions is

$$
V'(t) = -a\hat{S}^2(t) + be^{-a_p\tau}\hat{u}(t-\tau)\hat{S}(t) - ce^{-2a_p\tau}\hat{u}^2(t-\tau) - \delta_1 g(I(t))\hat{S}^2(t) \\
\qquad + [\delta_1\alpha I(t) + (\delta_2 - \delta_1\bar{S}^0)g(I(t))]\hat{S}(t) + \delta_2 g(I(t))\bar{S}^0 - \delta_2(a_I + \alpha + \mu)I(t),
$$

with the notations

$$
a := \delta_1(\eta + a_S) - \eta^2 > 0, \quad b := \delta_1(1 - \nu) + 2\eta\nu > 0 \quad \text{and} \quad c := e^{2a_p\tau} - \nu^2 > 0.
$$

We know from the proof of Theorem 6 that

$$
-a \hat{S}^2(t) + b e^{-a_p \tau} \hat{u}(t-\tau) \hat{S}(t) - c e^{-2 a_p \tau} \hat{u}^2(t-\tau) < 0.
$$

We have to find $\delta_2 > 0$ such that

$$
-\delta_1 g(I(t)) \hat{S}^2(t) + [\delta_1 \alpha I(t) + (\delta_2 - \delta_1 \bar{S}^0) g(I(t))] \hat{S}(t) + \delta_2 g(I(t)) \bar{S}^0 - \delta_2 (a_I + \alpha + \mu) I(t) < 0.
$$

Thanks to Lemma 1, we have $g(I) < g'(0)I$, for all $I > 0$. So,

$$
\delta_2 g(I(t)) \bar{S}^0 - \delta_2 (a_I + \alpha + \mu) I(t) < - \delta_2 (a_I + \alpha + \mu) (1 - \mathcal{R}_0) I(t) < 0.
$$

For the last expression, it is difficult to use the discriminant to find its sign, because
of the nonlinear function *g*. We consider the particular case *α* = 0. Then, we choose

$$
\delta_2 = \delta_1 \bar{S}^0.
$$

Hence, we obtain

$$
V'(t) = -a\hat{S}^2(t) + be^{-a_p\tau}\hat{u}(t-\tau)\hat{S}(t) - ce^{-2a_p\tau}\hat{u}^2(t-\tau)
$$
$$ -\delta_1 g(I(t)) \hat{S}^2(t) + \delta_2 g(I(t)) \tilde{S}^0 - \delta_2 (a_I + \alpha + \mu) I(t), \\ < -\alpha_1 \hat{S}^2(t) - \delta_2 (a_I + \alpha + \mu) (1 - \mathcal{R}_0) I(t), $$

with

$$ \alpha_1 := \frac{4ac - b^2}{4c} > 0. $$

The latter estimate is similar to the one we obtained in the proof of Theorem 6. Then, using the same approach, based on the integration of both sides of the inequality and applying Barbalat's lemma 7, we conclude that the disease-free equilibrium $\tilde{E}^0 := (\tilde{S}^0, 0, \tilde{u}^0)$ of the hybrid system (14) is globally asymptotically stable. $\square$

## 5.2 Global asymptotic stability of the endemic steady-state

In this section, we assume that $\mathcal{R}_0 > 1$, where $\mathcal{R}_0$ is given by (21). Then, we have the existence and uniqueness of the endemic steady-state $E^* = (S^*, I^*, u^*)$ of the system (14). In the following result, we show the global asymptotic stability of $E^*$ in the case where the nonlinear force of infection $g$ satisfies (5) and $\alpha = 0$.

**Theorem 9.** Assume that $g$ satisfies (5), $\alpha = 0$ and $\mathcal{R}_0 > 1$. Then, the endemic steady-state $E^*$ of the system (14) is globally asymptotically stable.

*Proof.* We proceed to the change of variables

$$ \tilde{S}(t) = S(t) - S^* \quad \text{and} \quad \tilde{u}(t) = u(t) - u^*. $$

We emphasize that the signs of $\tilde{S}(t)$ and $\tilde{u}(t)$ are unknown. We do not change the variable for $I(t)$ in order to preserve its positivity, which is necessary for the proof.

The system (14) becomes

$$ \left\{ \begin{aligned} \tilde{S}'(t) &= -(\eta + a_S + g(I(t)))\tilde{S}(t) - [g(I(t)) - g(I^*)]S^* + (1 - \nu)e^{-a_p\tau}\tilde{u}(t - \tau), \\ I'(t) &= g(I(t))[\tilde{S}(t) + S^*] - (a_I + \mu)I(t), \\ \tilde{u}(t) &= \eta\tilde{S}(t) + \nu e^{-a_p\tau}\tilde{u}(t - \tau). \end{aligned} \right. \qquad (39) $$

It is important to notice that $I$ satisfies

$$ I'(t) = \tilde{S}(t)g(I(t)) + (a_I + \mu)I^* \left[ \frac{g(I(t))}{g(I^*)} - \frac{I(t)}{I^*} \right]. $$

We consider the following Lyapunov function

$$ W(t) = \frac{\delta_1}{2} \tilde{S}^2(t) + \int_0^\tau \tilde{u}^2(t-s)ds + \delta_2 \int_{I^*}^{I(t)} \frac{g(v) - g(I^*)}{g(v)} dv, $$
with

$$ \delta_1 = 2 \frac{\eta(e^{2a_p\tau} - \nu) + a_S(e^{2a_p\tau} - \nu^2)}{(1 - \nu)^2} > 0 \quad \text{and} \quad \delta_2 = \delta_1 S^* > 0. $$

We notice that the function $H: [0, +\infty) \to [0, +\infty)$ defined by

$$ H(I) = \int_{I^*}^{I} \frac{g(v) - g(I^*)}{g(v)} dv, \quad I \ge 0, $$

satisfies

$$ H(I) > 0 \quad \text{if and only if} \quad I \neq I^*. $$

In fact, the function $H$ is convex, decreasing on $[0, I^*)$ and increasing on $(I^*, +\infty)$, with $H(I^*) = H'(I^*) = 0$. The next step is to compute the derivative of $W$ along the solutions. We have

$$ \begin{aligned} W'(t) &= \delta_1 \tilde{S}(t) \tilde{S}'(t) + \tilde{u}^2(t) - \tilde{u}^2(t - \tau) + \delta_2 I'(t) \left(1 - \frac{g(I^*)}{g(I(t))}\right), \\ &= \delta_1 \tilde{S}(t) \left[-(\eta + a_S) \tilde{S}(t) - \tilde{S}(t) g(I(t)) - S^* g(I(t)) + S^* g(I^*)\right] \\ &\quad + (1 - \nu) e^{-a_p \tau} \delta_1 \tilde{S}(t) \tilde{u}(t - \tau) + [\eta \tilde{S}(t) + \nu e^{-a_p \tau} \tilde{u}(t - \tau)]^2 - \tilde{u}^2(t - \tau) \\ &\quad + \delta_2 \tilde{S}(t) g(I(t)) \left(1 - \frac{g(I^*)}{g(I(t))}\right) \\ &\quad + \delta_2 (a_I + \mu) I^* \left(\frac{g(I(t))}{g(I^*)} - \frac{I(t)}{I^*}\right) \left(1 - \frac{g(I^*)}{g(I(t))}\right). \end{aligned} $$

After simplification, we get

$$ \begin{aligned} W'(t) = & -a \tilde{S}^2(t) + be^{-a_p \tau} \tilde{u}(t - \tau) \tilde{S}(t) - ce^{-2a_p \tau} \tilde{u}^2(t - \tau) - \delta_1 g(I(t)) \tilde{S}^2(t) \\ & - (\delta_1 S^* - \delta_2) g(I(t)) \tilde{S}(t) + (\delta_1 S^* - \delta_2) g(I^*) \tilde{S}(t) \\ & + \delta_2 (a_I + \mu) I^* \left(1 - \frac{I(t)g(I^*)}{I^*g(I(t))}\right) \left(\frac{g(I(t))}{g(I^*)} - 1\right), \end{aligned} $$

with

$$ a := \delta_1(\eta + a_S) - \eta^2 > 0, $$

$$ b := \delta_1(1 - \nu) + 2\eta\nu > 0, $$

$$ c := e^{2a_p \tau} - \nu^2 > 0. $$

Using the choice $\delta_2 = \delta_1 S^*$, we obtain,

$$ \begin{aligned} W'(t) = & -a \tilde{S}^2(t) + be^{-a_p \tau} \tilde{u}(t - \tau) \tilde{S}(t) - ce^{-2a_p \tau} \tilde{u}^2(t - \tau) - \delta_1 g(I(t)) \tilde{S}^2(t) \\ & + \delta_1 S^*(a_I + \mu) I^* \left( \frac{I(t)g(I^*)}{I^*g(I(t))} - 1 \right) \left( 1 - \frac{g(I(t))}{g(I^*)} \right). \end{aligned} $$
On one hand, we already know that

$$-a\tilde{S}^2(t) + be^{-a_p\tau}\tilde{u}(t-\tau)\tilde{S}(t) - ce^{-2a_p\tau}\tilde{u}^2(t-\tau) < 0.$$

On the other hand, thanks to Lemma 1, we have

$$\left(\frac{I(t)g(I^*)}{I^*g(I(t))} - 1\right) \left(1 - \frac{g(I(t))}{g(I^*)}\right) < 0, \quad I > 0, \quad I \neq I^*.$$

Once again, we obtain

$$W'(t) < -\alpha_1 \tilde{S}^2(t) - \delta_1 S^*(a_I + \mu) I^* \left(1 - \frac{I(t)g(I^*)}{I^*g(I(t))}\right) \left(1 - \frac{g(I(t))}{g(I^*)}\right) < 0.$$

Then, we have

$$\lim_{t \to +\infty} W(t) = \inf_{s \ge 0} W(s) := W^* \ge 0$$

and

$$\lim_{t \to +\infty} \int_{0}^{t} \tilde{S}^2(s) ds \le \frac{1}{\alpha_1} [W(0) - W^*].$$

By the same reasoning as for the proof of Theorem 6, we obtain that

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

Then, the expression of W implies that

$$\lim_{t \to +\infty} H(I(t)) = \lim_{t \to +\infty} \int_{I^*}^{I(t)} \frac{g(v) - g(I^*)}{g(v)} dv = \frac{W^*}{\delta_2}.$$

Thanks to the fluctuation lemma (see page 154, Lemma A.1 of [43]) there exists a sequence $t_k \to +\infty$ such that $\lim_{k \to +\infty} \tilde{S}'(t_k) = 0$. Then, the first equation of (39) implies that $\lim_{k \to +\infty} [g(I(t_k)) - g(I^*)] = 0$. From the continuity of g, we obtain that $\lim_{k \to +\infty} I(t_k) = I^*$. Consequently,

$$\lim_{k \to +\infty} \int_{I^*}^{I(t_k)} \frac{g(v) - g(I^*)}{g(v)} dv = 0.$$

We conclude that $W^* = 0$ and from the properties of the function H that $\lim_{t \to +\infty} I(t) = I^*$. Then, we obtain the globally asymptotic stability of the endemic steady-state $E^*$ of the system (14). □

Now suppose that the force of infection $g(I) = \beta I$ and $\alpha > 0$. Then, we also have the global asymptotic stability of the endemic steady-state $E^*$.

**Theorem 10.** Assume that $g(I) = \beta I$, $\alpha > 0$ and $R_0 > 1$. Then, the endemic steady-state $E^*$ of the system (14) is globally asymptotically stable.
Proof. In this case, we have

$$ \mathcal{R}_0 = \frac{\beta \bar{S}^0}{a_I + \alpha + \mu} > 1, $$

and

$$ S^* = \frac{a_I + \alpha + \mu}{\beta}, \quad u^* = \frac{\eta}{1 - \nu e^{-a_p \tau}} S^*, \quad I^* = \frac{b(a_I + \alpha + \mu)}{\beta \bar{S}^0 (a_I + \mu)} [\mathcal{R}_0 - 1]. $$

Consider again the Lyapunov function

$$ W(t) = \frac{\delta_1}{2} \tilde{S}^2(t) + \int_0^\tau \tilde{u}^2(t-s)ds + \delta_2 \left( I(t) - I^* - I^* \ln \frac{I(t)}{I^*} \right), $$

and in the previous results, we choose the parameters $\delta_1$ and $\delta_2$ as follows

$$ \delta_1 = 2 \frac{\eta(e^{2a_p \tau} - \nu) + a_S(e^{2a_p \tau} - \nu^2)}{(1 - \nu)^2} > 0 \quad \text{and} \quad \delta_2 = \delta_1 \left( \bar{S}^0 - \frac{\alpha}{\beta} \right). $$

$\delta_2 > 0$ because $\mathcal{R}_0 := \frac{\beta \bar{S}^0}{a_I + \alpha + \mu} > 1$. Furthermore, for $g(I) = \beta I$, we have

$$ \int_{I^*}^{I(t)} \frac{g(v) - g(I^*)}{g(v)} dv = I(t) - I^* - I^* \ln \frac{I(t)}{I^*}. $$

We then follow the same technique as before to show that the endemic steady-state $E^*$ is globally asymptotically stable. □

# 6 Discussion

The unit of time considered in the model depends on the disease in question, and can be, for example, the day, month or year. The units of the other parameters and variables used in the model are summarized in the following table.

<table><caption>**Table 1:** Units of parameters and variables</caption><thead><tr><th>Symbols</th><th>unit</th></tr></thead><tbody><tr><td>S, I, R, P</td><td>indiv</td></tr><tr><td>u</td><td>indiv.(time unit)<sup>-1</sup></td></tr><tr><td>τ</td><td>time unit</td></tr><tr><td>b</td><td>indiv.(time unit)<sup>-1</sup></td></tr><tr><td>a<sub>S</sub>, a<sub>I</sub>, a<sub>R</sub>, a<sub>p</sub>, μ, α, η</td><td>(time unit)<sup>-1</sup></td></tr><tr><td>β</td><td>(indiv.time unit)<sup>-1</sup></td></tr><tr><td>a</td><td>indiv<sup>-1</sup></td></tr><tr><td>ν</td><td>without unit</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Fig. 4:** Numerical simulation of the evolution of susceptible, infected, recovered and new protected individuals (SIR model). Parameters are: $b = 2$, $a_S = 0.1$, $\eta = 0.2$, $\tau = 5$, $\nu = 0.05$, $a_p = 0.1$, $\alpha = 0$, $\mu = 0.15$, $a_I = 0.1$ and $a_R = 0.15$. We have $\mathcal{R}_0 > 1$ and the force of infection $g(I) = \beta I/(1+I)$, with $\beta = 0.1$. The units of the parameters are in Table 1. We consider two different arbitrary initial conditions.</figcaption></figure>

The parameters common to the numerical simulations of Figures 4 and 5 are as follows $b = 2$, $a_S = 0.1$, $\eta = 0.2$, $\tau = 5$, $\nu = 0.05$, $a_p = 0.1$, $a_I = 0.1$ and $a_R = 0.15$. The force of infection chosen is $g(I) = \beta I/(1+I)$ with $\beta = 0.1$. The units of the parameters are summarized in Table 1. For each of the two figures, we consider two different arbitrary initial conditions. Figure 4 is done with $\alpha = 0$. This gives $\mathcal{R}_0 \approx 4.42 > 1$. And the endemic equilibrium is given by $E^* \approx (11.25, 0.07, 0.07, 2.27)$. We take $\mu = 0$ for Figure 5. The basic reproduction is then $\mathcal{R}_0 \approx 5.52 > 1$. We can see that the solutions reaches the endemic equilibrium which is $E^* \approx (11, 0.09, 2.67)$. Figure 6 shows the influence of the main parameters, $\tau$, $\nu$ and $\eta$ on the basic reproduction number $\mathcal{R}_0$. It can be seen that the duration of the protection phase $\tau$ allows $\mathcal{R}_0$ to drop below 1 more readily than the other two.

In this work, we have reduced an epidemiological model with a temporary pro-tection compartment, composed of a system of ordinary differential equations and a partial differential equation structured in time since individuals enter the protection phase, to a hybrid mathematical system involving differential and difference equations. The proposed system can model SIR and SIS epidemic scenarios with a temporary protection phase, and a mass action or a nonlinear force of infection with saturation,
<figure><img src="image_7.png" /><figcaption>**Fig. 5:** Numerical simulation of the evolution of susceptible, infected individuals and new protected individuals (SIS model). Parameters are: $b = 2$, $a_S = 0.1$, $\eta = 0.2$, $\tau = 5$, $\nu = 0.05$, $a_p = 0.1$, $\alpha = 0.1$, $\mu = 0$, $a_I = 0.1$ and $a_R = 0.15$. We have $\mathcal{R}_0 > 1$ and the force of infection $g(I) = \beta I/(1+I)$, with $\beta = 0.1$. The units of the parameters are in Table 1. We consider two different arbitrary initial conditions.</figcaption></figure>

(5). This work is a continuation of recently published articles [1–3, 6]. These models are based on a coupling between differential equations and continuous-time difference equations. The idea of the paper is to promote these equations by proposing a study of some generality allowing to capture several models at once. We have considered the case of the SIS model, i.e., by supposing that there is a return of the infected individuals (expressed by the parameter $\alpha I$ with $\alpha > 0$) to the class of susceptible individuals. We have also studied the case $\alpha = 0$ which means that individuals can not lost their immunity. We could also have chosen $\mu = 0$ to ensure that no individual is immune for life. In addition, we have taken the case of a nonlinear incidence (expressed by the function $g$). For the hybrid system, we determined the disease-free and endemic steady-states and gave the basic reproduction number. We obtained a transcendental characteristic equation that is difficult to handle, but from which we established several results on local asymptotic stability. We then established global asymptotic stability results for both steady-states by constructing appropriate Lyapunov functions. However, we had difficulty in constructing a Lyapunov function for the case where both $\alpha > 0$ and the force of infection $g$ is nonlinear. Despite this difficulty, we have proved global asymptotic stability by Lyapunov functions for the case
<figure><img src="image_8.png" /><figcaption>**Fig. 6:** $\mathcal{R}_0$ as a function of the three parameters $\nu$, $\tau$ and $\eta$, with $b = 2$, $a_I = 0.1$, $\alpha = 0.1$, $\mu = 0.15$, $a_S = 0.1$, $a_p = 0.7$, $g'(0) = \beta = 0.04$, where $g(I) = \beta I$ or $g(I) = \beta I/(1 + aI)$, and (1) $\tau = 0.5$, $\eta = 0.2$, (2) $\nu = 0.05$, $\eta = 0.2$ and (3) $\nu = 0.05$, $\tau = 0.5$. The units of the parameters are in Table 1.</figcaption></figure>

of the SIR model with general incidence and the case of the SIS model with bilinear incidence. The hybrid model introduced in this work provides a range of modeling possibilities, while also presenting many new mathematical challenges.

**Acknowledgments.** A. Chekroun thanks grant, PRFU: C00L03UN29012022002, from DGRSDT of Algeria.

**Declarations**

The authors declare there is no conflict of interest.

**References**

[1] 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(2):1329, 2020.

[2] 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, 1–13, 2022.

[3] M. Adimy, A. Chekroun, T. Kuniya, and H. Meghelli. Global stability of a SEIR discrete delay differential-difference system with protection phase. Mathematical Methods in the Applied Sciences, 1–21, 2023.

[4] 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(9):2765–2791, 2015.
[5] M. Adimy, A. Chekroun, and T.M. Touaoula. Global asymptotic stability for an age-structured model of hematopoietic stem cell dynamics. Applicable Analysis, 96(3):429-440, 2016.

[6] M. Adimy, J. Molina, L. Pujo-Menjouet, G. Ranson and J. Wu. Forecasting the Effect of Pre-Exposure Prophylaxis (PrEP) on HIV Propagation with a System of Differential-Difference Equations with Delay. Mathematics, 10(21): 1-24, 2022.

[7] M. Alexander, C. Bowman, S. Moghadas, R. Summers, A. Gumel, and B. Sahai. A vaccination model for transmission dynamics of influenza. SIAM Journal on Applied Dynamical Systems, 3(4):503-524, 2004.

[8] L.J. Allen and P. van den Driessche. The basic reproduction number in some discrete-time epidemic models. Journal of Difference Equations and Applications, 14(10-11):1127-1147, 2008.

[9] R.M. Anderson. Discussion: The Kermack-McKendrick epidemic threshold theorem. Bulletin of Mathematical Biology, 53(1):3-32, 1991.

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

[11] E. Beretta and Y. Kuang. Geometric Stability Switch Criteria in Delay Differential Systems with Delay Dependent Parameters. SIAM Journal on Mathematical Analysis, 33:1144-1165, 2002.

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

[13] D. Bernoulli. Essai d'une nouvelle analyse de la mortalité causée par la petite vérole, et des avantages de l'inoculation pour la prévenir. Mémoires de Mathématique et de Physique, tirés des registres de l'Académie Royale des Sciences, Paris, 1766.

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

[15] B. Buonomo, A. d'Onofrio, and D. Lacitignola. Global stability of an SIR epidemic model with information dependent vaccination. Mathematical Biosciences, 216(1):9-16, 2008.

[16] V. Capasso and G. Serio. A generalization of the Kermack-McKendrick deterministic epidemic model. Mathematical Biosciences, 42(1):43-61, Nov. 1978.
[17] A. Chekroun, M. N. Frioui, T. Kuniya, and T.M. Touaoula. Global stability of an age-structured epidemic model with general Lyapunov functional. Mathematical Biosciences and Engineering, 16:1525-1553, 2019.

[18] A. Chekroun, M.N. Frioui, T. Kuniya, and T.M. Touaoula. Mathematical analysis of an age structured heroin-cocaine epidemic model. Discrete and Continuous Dynamical Systems - B, 25(11):4449-4477, 2020.

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

[20] M.A. Cruz and J.K. Hale. Stability of functional differential equations of neutral type. Journal of Differential Equations, 7(2):334-355, 1970.

[21] O. Diekmann and J. Heesterbeek. Mathematical Epidemiology of Infective Diseases: Model Building, Analysis and Interpretation. Wiley, New York, 2000.

[22] X. Duan, S. Yuan, and X. Li. Global stability of an SVIR model with age of vaccination. Applied Mathematics and Computation, 226:528-540, 2014.

[23] 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 and Continuous Dynamical Systems - B, 15:61-74, 2011.

[24] 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:2120-2133, 2012.

[25] X.X. Feng, L.L. Liu, and F.Q. Zhang. Dynamical Behavior of SEIR-SVS Epidemic Models with Nonlinear Incidence and Vaccination. Acta Mathematicae Applicatae Sinica, English Series, 38(2):282-303, 2022.

[26] S. Gao, L. Chen, J.J. Nieto, and A. Torres. Analysis of a delayed epidemic model with pulse vaccination and saturation incidence. Vaccine, 24(35):6037-6045, 2006.

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

[28] J.K. Hale, and S.M. Verduyn Lunel, Introduction to Functional Differential Equations. Springer New York, 1993.

[29] F.Z. Hathout, T.M. Touaoula, and S. Djilali. Mathematical analysis of a triple age dependent epidemiological model with including a protection strategy. Discrete and Continuous Dynamical Systems-B, 27(12): 7409-7443, 2022.

[30] M. Iannelli, M. Martcheva, and X.Z. Li. Strain replacement in an epidemic model with super-infection and perfect vaccination. Mathematical Biosciences, 195(1):23-46, 2005.
[31] W.O. Kermack and A.G. McKendrick. A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London, Series A, 115(772):700-721, 1927.

[32] H.K. Khalil. Nonlinear systems; 3rd ed. Prentice-Hall, Upper Saddle River, NJ, 2002.

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

[34] 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:38-47, 2010.

[35] X.-Z. Li, J. Wang, and M. Ghosh. Stability and bifurcation of an SIVS epidemic model with treatment and age of vaccination. Applied Mathematical Modelling, 34(2):437-450, 2010.

[36] 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:29-43, 2018.

[37] X. Liu, Y. Takeuchi, and S. Iwami. SVIR epidemic models with vaccination strategies. Journal of Theoretical Biology, 253(1):1-11, 2008.

[38] P. Magal, C.C. McCluskey, and G.F. Webb. Lyapunov functional and global asymptotic stability for an infection-age model. Applicable Analysis, 89(7):1109-1140, 2010.

[39] M. Martcheva. An Introduction to Mathematical Epidemiology. Springer Publishing Company, Incorporated, 1st edition, 2015.

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

[41] J.D. Murray. Mathematical Biology: I. An Introduction. Springer, New York, 3rd edition, 2002.

[42] R. Ross. The prevention of malaria. John Murray, 1911.

[43] H. Smith. An Introduction to Delay Differential Equations with Applications to the Life Sciences. Springer, New York, 2011.

[44] S. Tipsri and W. Chinviriyasit. The effect of time delay on the dynamics of an SEIR model with nonlinear incidence. Chaos, Solitons & Fractals, 75:153-172, 2015.
[45] Y. Xiao and S. Tang. Dynamics of infection with nonlinear incidence in a simple vaccination model. *Nonlinear Analysis: Real World Applications*, 11(5):4154–4163, 2010.

[46] 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:332–342, 2015.