# Reinfection induced multistability in an epidemic model

Pierre-Alexandre Bliman, Marcel Fang

► To cite this version:

Pierre-Alexandre Bliman, Marcel Fang. Reinfection induced multistability in an epidemic model. Nonlinear Science, 2025. hal-05399128

HAL Id: hal-05399128

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

Submitted on 5 Dec 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Reinfection induced multistability in an epidemic model*

Pierre-Alexandre Bliman¹ and Marcel Fang¹

¹Sorbonne Université, Université Paris Cité, CNRS, Inria, Laboratoire Jacques-Louis Lions, LJLL, EPC MUSCLEES, F-75005 Paris, France

December 4, 2025

## Abstract

We consider in this paper the effects induced on the dynamics of disease transmission by differences between primary and subsequent infections (i.e. reinfections), due e.g. to enhancement or weakening of the susceptibility or infectivity. To this end, an 8-dimensional ‘two-stage’ SEIRS reinfection model is considered, extending the classical 4-dimensional SEIRS model. We characterize the steady states of the model according to the basic reproduction number $R_0$. It is shown that the reinfection induced heterogeneity may cause up to two endemic equilibria when $R_0 \le 1$, and up to three endemic equilibria otherwise. Specifically, this suggests that the model may present backward bifurcation for $R_0 < 1$, a feature quite important from the point of view of disease control; and two successive saddle node bifurcations for $R_0 > 1$. Simulations confirm this situation. The disease persistence of the model is also examined. Finally, we prove, for two specific SIRI and SEIRE models accounting for partial immunity and demography, the asymptotic convergence of every trajectory to a steady state. Notably, these particular cases still admit up to two endemic equilibria (when $R_0 < 1$), which makes this result non trivial. The proof is based on an application of Li & Muldowney theory to epidemiological systems with multiple endemic equilibriums.

**Keywords**— Compartmental models in epidemiology; reinfection models; multistability; backward bifurcations; persistence; asymptotic convergence; Li & Muldowney theory; k-contraction

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Disease-free equilibrium and basic reproduction number R<sub>0</sub></td><td>4</td></tr><tr><td>3</td><td>Endemic equilibria</td><td>6</td></tr><tr><td>3.1</td><td>Equilibria description</td><td>6</td></tr><tr><td>3.2</td><td>Numerical illustrations</td><td>8</td></tr><tr><td>3.2.1</td><td>The case R<sub>0</sub> &lt; 1</td><td>8</td></tr><tr><td>3.2.2</td><td>The case R<sub>0</sub> &gt; 1</td><td>9</td></tr><tr><td>3.2.3</td><td>Asymptotic convergence</td><td>10</td></tr><tr><td>4</td><td>Persistence</td><td>10</td></tr><tr><td>5</td><td>Global analysis of reinfection models with partial immunity</td><td>11</td></tr><tr><td>5.1</td><td>SIRI model with partial immunity</td><td>12</td></tr><tr><td>5.2</td><td>SEIRE model with partial immunity</td><td>13</td></tr><tr><td>6</td><td>Discussion</td><td>14</td></tr><tr><td>7</td><td>Conclusion</td><td>17</td></tr><tr><td></td><td>References</td><td>17</td></tr></table>

<footnote>*The financial support of French National Research Agency (ANR) is acknowledged (ANR-23-CE48-0004)</footnote>
A Elements of Li and Muldowney theory 19

<table><tr><td>B Proofs</td><td>21</td></tr><tr><td>B.1 Proof of Proposition 1</td><td>21</td></tr><tr><td>B.2 Proofs of Theorems 3 to 5</td><td>22</td></tr><tr><td>B.3 Proof of Theorem 2</td><td>26</td></tr><tr><td>B.4 Proof of Theorem 6</td><td>27</td></tr><tr><td>B.5 Proofs of Theorems 7 and 10</td><td>29</td></tr><tr><td>B.6 Proofs of Theorems 9 and 12</td><td>31</td></tr></table>

# 1 Introduction

Since the SIR model of Kermack and McKendrick, the field of mathematical epidemiology has evolved considerably, incorporating increasingly complex factors such as quarantine, vaccine, hosts heterogeneity, variants, spatial dynamics or stochastic effects [33]. Recently, a particular attention has been paid to reinfection models, which are driven by the study of the transmission of several diseases for which reinfection cannot be ignored, for example tuberculosis [13], pertussis [17], or COVID-19 [37]. In most of the available reinfection models, the infection and reinfection processes are assumed to behave essentially the same, and up to our knowledge there has been no specific investigation of the effects of the heterogeneity induced by differences in these processes, on the dynamics of the disease spread. Correspondingly, in the classical reinfection models [21] the population is divided between compartments regardless of their history of infections, an approach not suitable for considering diseases with distinct reinfection behaviors.

For instance, in the SEIRS model (1) (adapted from [3] by normalization of the total population to 1), the individuals in the susceptible state $S$ move to the exposed compartment $E$ upon contact with infectious individuals, then become infectious in compartment $I$, and finally recover with temporary immunity in compartment $R$. Once this immunity wanes, they return to the susceptible state, repeating the same cycle $S \to E \to I \to R \to S$ without keeping memory of the past infections.

$$ 
\begin{aligned}
\dot{S} &= \mu - \beta IS + \omega R - \mu S, \\
\dot{E} &= \beta IS - (\sigma + \mu)E, \\
\dot{I} &= \sigma E - (\gamma + \mu)I, \\
\dot{R} &= \gamma I - (\omega + \mu)R,
\end{aligned}
\qquad (1) $$

The parameters $\mu, \beta, \sigma, \gamma$ and $\omega$ are positive and represent respectively the natural birth and death rate, the transmission rate, the rate at which the exposed become infected, the recovery rate and the rate characterizing the immunity loss. The mathematical properties of this model have been analyzed [24, 31]. Its basic reproduction number is equal to

$$ \mathcal{R}_0 := \frac{\beta_{1,1}}{\gamma_1 + \mu} \frac{\sigma_1}{\sigma_1 + \mu}. \qquad (2) $$

The system always contains a unique disease-free equilibrium, and a unique endemic equilibrium also exists if and only if $\mathcal{R}_0 > 1$. Moreover, the disease-free equilibrium is globally asymptotically stable (GAS) in the biologically feasible set for $\mathcal{R}_0 < 1$; while for $\mathcal{R}_0 > 1$, the endemic equilibrium is GAS in the interior of the biologically feasible set [24, 31].

Notice that $\mathcal{R}_0$ defined in (2) does not depend upon the immunity loss $\omega$. Taking this value equal to zero yields the SEIR model, already analyzed in [25]. For this model, the same results hold concerning the equilibrium points.

Several works have investigated reinfection models with partial immunity, providing in fact a simple example of the type of heterogeneity we aim to study. Notably, the SIRI model (3) with partial immunity has been investigated in [14].

$$ 
\begin{aligned}
\dot{S} &= \mu - \beta IS - \mu S, \\
\dot{I} &= \beta IS + \psi\beta IR - (\gamma + \mu)I, \\
\dot{R} &= \gamma I - \psi\beta IR - \mu R.
\end{aligned}
\qquad (3) $$

It is demonstrated therein (see also Section 5 below) the existence of a threshold on the reinfection parameter $\psi$ by which the system transits from a low-infectivity to a high-infectivity equilibrium. The steady states and the asymptotic convergence of the SIRI model has also been addressed in [38] in the case where the demography effects are negligible ($\mu = 0$). To gain further insight into the effects of heterogeneity between infection and reinfection, we introduce in this paper the two-stage SEIRS model (4), for which an illustration is given by the
<figure><img src="image_2.png" /><figcaption>Figure 1: Flowchart of the SEIRS model (1)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Flowchart of the two-stage SEIRS model (4)</figcaption></figure>

flowchart in Figure 2. As in the standard SEIRS model (1), the classes of compartments S, E, I, R represent the susceptible, exposed, infectious and recovered. The two-stage model (4) explicitly separates the primary infections from the subsequent infections. This is rendered through the subscripts 1 and r (as ‘remaining’).

$$
\begin{align}
\dot{S}_1 &= \mu - (\beta_{1,1}I_1 + \beta_{1,r}I_r + \mu)S_1, \nonumber \\
\dot{E}_1 &= (\beta_{1,1}I_1 + \beta_{1,r}I_r)S_1 - (\sigma_1 + \mu)E_1, \nonumber \\
\dot{I}_1 &= \sigma_1 E_1 - (\gamma_1 + \mu)I_1, \nonumber \\
\dot{R}_1 &= \gamma_1 I_1 - (\omega_1 + \mu)R_1, \nonumber \\
\dot{S}_r &= \omega_1 R_1 + \omega_r R_r - (\beta_{r,1}I_1 + \beta_{r,r}I_r + \mu)S_r, \nonumber \\
\dot{E}_r &= (\beta_{r,1}I_1 + \beta_{r,r}I_r)S_r - (\sigma_r + \mu)E_r, \nonumber \\
\dot{I}_r &= \sigma_r E_r - (\gamma_r + \mu)I_r, \nonumber \\
\dot{R}_r &= \gamma_r I_r - (\omega_r + \mu)R_r, \tag{4}
\end{align}
$$

The notations used here are adapted from the ones used for (1). More precisely, $\beta_{1,1}$, $\beta_{1,r}$ and $\beta_{r,1}$, $\beta_{r,r}$ are transmission rates, respectively for primary infections acquired from primary infected and from reinfected; and for reinfections transmitted by primary infected and by reinfected; the parameters $\sigma_1$, $\sigma_r$ are respectively the rates at which the primary and secondary exposed hosts become infective; $\gamma_1$, $\gamma_r$ are rates of recovery of the primo-infected and reinfected; and $\omega_1$, $\omega_r$ represent the loss of immunity, depending again on whether the infection is primary or secondary.

It is noteworthy to mention that, concurrently with our study and driven by the need for more complex models during the COVID-19 pandemic, two-stage reinfection models have been proposed independently in some papers of mathematical epidemiology. For instance, an age-specific two-stage SEIRS model has been utilized for statistical inference of the true infection level in multiple waves of the SARS-CoV-2 pandemic [4], which reveals major under-reporting and a shifting burden toward children and adolescents during the pandemic. In addition, Kaklamanos et al. [19] investigated the dynamics of a similar but different two-stage SIRS model. The
latter, which does not accommodate demographic effects, features two different levels of immunity, accessible after primary and secondary infections, and different time-scales for their appearance and disappearance. Analysis conducted thanks to singular perturbation theory, allows to exhibit backward bifurcation, in situations which however do not match the ones described by our model. Finally, two-stage models have also found their application for the modelling of animal infections, such as bovine tuberculosis (BTB) or Aujesky's disease (pseudorabies virus) in pigs [1, 15]. These considerations underscore the significance of the two-stage model for diseases that permit reinfections.

In the present paper, our goal is to describe precisely the asymptotic behaviour of some nonlinear systems from mathematical epidemiology, in order to investigate the effects induced by reinfections. Based on the generic model (4) introduced before, our study is less focused on modelling aspects than in describing precisely the dynamical behaviour, especially as regards the convergence to equilibrium in a multistability context. We will first describe the steady states of the model, according to the parameters involved in the reinfection. As we shall see, the heterogeneity induced by the reinfection may lead to a multistable system presenting additional endemic equilibria, even when the basic reproduction number $R_0$ is smaller than 1. This matter is important from a public health perspective, as it may complicate the disease control policies. The study of equilibria also underlines the existence of a reinfection threshold, in line with the one underlined by Gomes et al. [13] for the SIRI reinfection model. In a second time, the persistence of the system is analyzed, in the mathematical framework developed in [12]. Persistence is a well-known topic in mathematical biology [40, 41], and particularly in epidemiology, as it is related to the issue of disease eradication [13, 32]. Finally, the question of the asymptotic behavior of the two-stage model is another focal point of this paper. Indeed, the two-stage model often exhibits multiple endemic equilibria, and this makes the analysis of asymptotic behavior a complicated problem. We will prove, for some reinfection models (SIRI and SEIRE) with partial immunity that are subclasses of the general model, the asymptotic convergence of every trajectory to a steady state. It is noteworthy to mention that these particular models still present up to two endemic equilibria when $R_0$ is subcritical, which usually hinders global asymptotic analysis for this type of model. The convergence result is based on an adaptation of Li & Muldowney theory [27]. Originally, the theory has been successfully applied for the global asymptotic stability problem in several epidemiological models [10, 25, 31], where it is used to demonstrate the convergence to an endemic equilibrium (which is assumed to be the unique equilibrium in the invariant interior of the domain). Moreover, and importantly, it has been noticed that the requirement of uniqueness here can be lifted — a fact that has recently sparked considerable interest in the theory, particularly regarding the global asymptotic convergence of nonlinear systems with multiple equilibria [43]. Therefore, the convergence result in the present paper further highlight the relevance of Li & Muldowney theory in this regard. An explanation of why the method fails in the general setting of (4) will also be given.

The analysis begins in Section 2 with the description of the disease-free equilibrium and the computation of the basic reproduction number $R_0$ of the model. Next, in Section 3, comprehensive results on the existence and number of endemic equilibria are established. It is shown therein that the model may exhibit backward bifurcation when $R_0 < 1$. On the other hand, the model may present saddle node bifurcations and up to three distinct endemic equilibria simultaneously when $R_0 > 1$. The results are moreover illustrated with several numerical simulations. Disease persistence of the model is then considered in Section 4. In Section 5, the findings on the equilibria are applied to several particular cases of the general model, where the heterogeneity consists only of partial immunity, for which we establish also the asymptotic convergence of every trajectory to a steady state. The core of the proof relies on the Li and Muldowney's geometrical approach [27], for which a summary is given in Appendix A. Also, for the reader's ease, the proof of several results are reported in Appendix B. Finally, the paper includes a discussion section, Section 6, in which we explore the implications of our findings for real-life diseases. We also provide in the section a geometrical explanation of why the Li and Muldowney theory does not apply to the general two-stage SEIRS model.

## 2 Disease-free equilibrium and basic reproduction number $R_0$

We begin our journey with an easy computation showing that the system (4) admits a unique disease-free equilibrium, i.e., an equilibrium in which $E_1 = I_1 = E_r = I_r = 0$. Additionally, we may characterize the basic reproduction number $R_0$ by the threshold condition for local stability of the DFE. This value is the same than for the SEIRS model (1) presented below, and remarkably, it depends solely on coefficients intervening in the primary infection.

**Theorem 1.** System (4) admits the unique disease-free equilibrium (1, 0, 0, 0, 0, 0, 0). Moreover, the DFE is locally asymptotically stable for $R_0 < 1$ and unstable when $R_0 > 1$, for $R_0$ defined in (2).

*Proof.* By definition, $E_1 = I_1 = E_r = I_r = 0$ at disease-free equilibrium, which means that any DFE $E_0$ has the
form

$$E_0 = (S_1^0, 0, 0, R_1^0, S_r^0, 0, 0, R_1^0).$$

Moreover, the remaining values $S_1^0, R_1^0, S_r^0, R_1^0$ of the steady state fulfill the equations

$$\mu - \mu S_1^* = 0, \quad -(\omega_1 + \mu)R_1^* = 0, \quad \omega_1 R_1^* + \omega_r R_r^* - \mu S_r^* = 0, \quad -(\omega_r + \mu)R_r^* = 0.$$

Hence, $S_1^0 = 1$ and $R_1^0 = S_r^0 = R_1^0 = 0$. The Jacobian of the system (4) at the DFE is

$$A := \begin{pmatrix} -\mu & 0 & -\beta_{1,1} & 0 & 0 & 0 & -\beta_{1,r} & 0 \\ 0 & -(\sigma_1 + \mu) & \beta_{1,1} & 0 & 0 & 0 & \beta_{1,r} & 0 \\ 0 & \sigma_1 & -(\gamma_1 + \mu) & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \gamma_1 & -(\omega_1 + \mu) & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \omega_1 & -\mu & 0 & 0 & \omega_r \\ 0 & 0 & 0 & 0 & 0 & -(\sigma_r + \mu) & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \sigma_r & -(\gamma_r + \mu) & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & \gamma_r & -(\omega_r + \mu) \end{pmatrix}.$$

Let us compute the characteristic polynomial associated to the Jacobian matrix.

$$\det(A - \lambda I_8) = (-\mu - \lambda)^2 (-\omega_r - \mu - \lambda) (-\gamma_r - \mu - \lambda) (-\sigma_r - \mu - \lambda) (-\omega_1 - \mu - \lambda) \\ \times \begin{vmatrix} -(\sigma_1 + \mu + \lambda) & \beta_{1,1} \\ \sigma_1 & -(\gamma_1 + \mu + \lambda) \end{vmatrix}.$$

The factorization shows that at least 6 eigenvalues counting multiplicity exhibited in the factor above are negative, and the remaining two are the eigenvalues of the matrix

$$\begin{pmatrix} -(\sigma_1 + \mu) & \beta_{1,1} \\ \sigma_1 & -(\gamma_1 + \mu) \end{pmatrix}.$$

As the trace of this matrix is negative, it is necessary and sufficient that the determinant of this matrix is positive to have two negative eigenvalues, assumption which is equivalent to

$$(\sigma_1 + \mu)(\gamma_1 + \mu) - \beta_{1,1}\sigma_1 > 0 \Leftrightarrow \mathcal{R}_0 < 1,$$

where $\mathcal{R}_0$ is defined in (2).

Additionally, one may verify, using the method of next generation matrix [42], that $\mathcal{R}_0$ is indeed the basic reproduction number of the model. The system (4) has 4 infected compartments $E_1, I_1, E_r, I_r$, which will be alternatively denoted $x_1, x_2, x_3, x_4$ respectively for simplicity. On the other hand, the remaining compartments $S_1, R_1, S_r, R_r$ are denoted respectively $x_5, x_6, x_7, x_8$. We write now the dynamics of $x := (x_i)_{1 \le i \le 8}$ in the form

$$\dot{x}_i = F_i(x) - V_i(x).$$

Here, $F_i$ is the rate of appearance of new infections in the compartment $x_i$ and $V_i = V_i^- - V_i^+$, where $V_i^+$ (resp. $V_i^-$) is the rate of transition into by all others means (resp. out of) the compartment $x_i$. Denoting $F(x) := (F_i(x))_{1 \le i \le 8}$ and $V(x) := (V_i)_{1 \le i \le 8}$, and the Jacobians of $F$ and $V$ at the disease-free equilibrium are

$$DF(E_0) = \begin{pmatrix} \mathbb{F} & \mathbf{0}_{4 \times 4} \\ \mathbf{0}_{4 \times 4} & \mathbf{0}_{4 \times 4} \end{pmatrix}, \text{ where } \mathbb{F} = \begin{pmatrix} 0 & \beta_{1,1} & 0 & \beta_{1,r} \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix},$$

$$DV(E_0) = \begin{pmatrix} \mathbb{V} & \mathbf{0}_{4 \times 4} \\ J_1 & J_2 \end{pmatrix}, \text{ where } \mathbb{V} = \begin{pmatrix} \sigma_1 + \mu & 0 & 0 & 0 \\ -\sigma_1 & \gamma_1 + \mu & 0 & 0 \\ 0 & 0 & \sigma_r + \mu & 0 \\ 0 & 0 & -\sigma_r & \gamma_r + \mu \end{pmatrix}.$$

In addition,

$$\mathbb{V}^{-1} = \begin{pmatrix} (\sigma_1 + \mu)^{-1} & 0 & 0 & 0 \\ \sigma_1(\sigma_1 + \mu)^{-1}(\gamma_1 + \mu)^{-1} & (\gamma_1 + \mu)^{-1} & 0 & 0 \\ 0 & 0 & (\sigma_r + \mu)^{-1} & 0 \\ 0 & 0 & \sigma_r(\sigma_r + \mu)^{-1}(\gamma_r + \mu)^{-1} & (\gamma_r + \mu)^{-1} \end{pmatrix}.$$

The basic reproduction number of the model (4) is then given by the spectral value of $\mathbb{F}\mathbb{V}^{-1}$, which is equal to $\frac{\beta_{1,1} \sigma_1}{\gamma_1 + \mu \sigma_1 + \mu}$. $\square$
# 3 Endemic equilibria

We tackle now the trickier issue of characterizing the endemic equilibria of (4). After reducing the number of parameters through an appropriate reparameterization, the number of endemic equilibria is provided in Section 3.1. Numerical simulations presented in Section 3.2 corroborate these findings.

## 3.1 Equilibria description

To start with, let us introduce the dimensionless coefficients

$$A_1 = \frac{\mu}{\gamma_1 + \mu} \frac{\sigma_1}{\sigma_1 + \mu}, \quad B_1 = \frac{\beta_{1,r}}{\beta_{1,1}}, \quad C_1 = \frac{\mu}{\beta_{1,1}}, \\ A_r = \frac{\sigma_r}{\sigma_r + \mu} \frac{\gamma_r}{\gamma_r + \mu} \frac{\omega_r}{\omega_r + \mu}, \quad B_r = \frac{\beta_{r,1}}{\beta_{r,r}}, \quad C_r = \frac{\mu}{\beta_{r,r}}, \quad D_r = \frac{\omega_1}{\omega_r} \frac{\omega_r + \mu}{\omega_1 + \mu} \frac{\gamma_1}{\gamma_r} \quad (5)$$

Clearly, the coefficients introduced are all positive, and $A_1, A_r < 1$. With these definitions, it can be demonstrated (see the proof of Proposition 1 in Appendix B.1) that, at any endemic equilibrium $(S_1, E_1^*, I_1^*, R_1^*, S_r^*, E_r^*, I_r^*, R_r^*)$, $I_1^*$ and $I_r^*$ must fulfill the equations

$$I_1^* = A_1 \frac{I_1^* + B_1 I_r^*}{I_1^* + B_1 I_r^* + C_1}, \quad I_r^* = A_r \frac{B_r I_1^* + I_r^*}{B_r I_1^* + I_r^* + C_r} (D_r I_1^* + I_r^*). \quad (6)$$

Introducing in addition the following coefficients, also dimensionless:

$$a_1 = A_1 - C_1, \quad b_1 = B_1 > 0, \quad c_1 = A_1 B_1 > 0, \\ a_r = \frac{C_r}{1 - A_r} > 0, \quad b_r = \frac{1}{1 - A_r} (A_r (B_r + D_r) - B_r), \quad c_r = \frac{A_r B_r D_r}{1 - A_r} > 0, \quad (7)$$

and denoting for simplicity

$$x := I_1^*, \quad y := I_r^*,$$

the equations (6) can be equivalently written as

$$x^2 + (b_1 y - a_1)x - c_1 y = 0, \quad y^2 + (a_r - b_r x)y - c_r x^2 = 0. \quad (8)$$

Notice that the signs of $a_1$ and $b_r$ are unknown in general. Moreover

$$a_1 > 0 \Leftrightarrow \mathcal{R}_0 > 1, \quad a_1 < 0 \Leftrightarrow \mathcal{R}_0 < 1, \quad (9)$$

and in any case $a_1 b_1 < c_1$ is verified. In addition, as shown in [11, Proposition 2.5.3], the coefficients defined in (7) fulfill the inequality

$$\frac{c_1}{b_1} < \frac{2}{2 + b_r + \sqrt{b_r^2 + 4c_r}}. \quad (10)$$

It is noteworthy that, for any quintuplet $(a_1, b_1, c_1, a_r, b_r, c_r)$ fulfilling the inequalities $a_1 < \frac{c_1}{b_1}$, (10) and $b_1, c_1, a_r, c_r > 0$, there exist (at least) one assortment of the 11 positive parameters $\mu, \beta_{1,1}, \beta_{1,r}, \beta_{r,1}, \beta_{r,r}, \sigma_1, \sigma_r, \gamma_1, \gamma_r, \omega_1, \omega_r$ to which it corresponds by (5) and (7), see [11, Proposition 2.5.3]. Therefore, the reparameterization performed above provides exactly the same dynamical behaviours than the initial one.

Finally, the new coefficients in (7) allow us to establish a one-to-one correspondence of the endemic equilibria with the roots of a cubic polynomial.

**Proposition 1.** Let the polynomial $\mathcal{P}$ be defined as

$$\mathcal{P}(x) := x(x - a_1)^2 - (x - a_1)(c_1 - b_1 x)(b_r x - a_r) - c_r(c_1 - b_1 x)^2 x, \quad (11)$$

with the coefficients $a_1, b_1, c_1, a_r, b_r, c_r$ defined by (5) and (7). The endemic equilibria of system (4) are in one-to-one correspondence with the solutions of the equation $\mathcal{P}(x) = 0$ in the interval $(0, 1) \cap (a_1, \frac{c_1}{b_1})$, which give all possible values of $I_1$ at endemic equilibrium.

The detailed proof of Proposition 1 is reported in Appendix B.1. For the subsequent analysis of the polynomial $\mathcal{P}$, we consider $a_r > 0$ as a bifurcation parameter. In fact, writing separately the terms depending on $a_r$ from the rest of the polynomial $\mathcal{P}$ gives

$$\mathcal{P}(x) = xP(x) + a_rQ(x), \\ P(x) := (x - a_1)^2 - b_r(x - a_1)(c_1 - b_1x) - c_r(c_1 - b_1x)^2, \quad Q(x) := (x - a_1)(c_1 - b_1x).$$

As the roots of $Q$ are well known, this decomposition helps to visualize the roots of $\mathcal{P}$ geometrically and motivates the choice of taking $a_r$ as the bifurcation parameter.
**Remark 1.** Notice that the parameter $a_r$ is independent of the value of $R_0$, which is expressed only as a function of the parameters intervening in the first stage of the infection dynamic (see the formula (2) that provides its value). Keeping the coefficients $\mu, \beta_{1,1}, \beta_{1,r}, \sigma_1, \sigma_r, \gamma_1, \gamma_r, \omega_1, \omega_r$ fixed and considering (5) and (7), we see that a variation in $\beta_{r,1}$ and $\beta_{r,r}$ of the same proportion amounts to modify $a_r$ in an inversely proportional way. Therefore, an increase in the parameter $a_r$ can be interpreted as a decrease in the susceptibility and/or infectivity to reinfections.

In order to study the roots of $\mathcal{P}$, we recall the following definition and results on cubic polynomials, see e.g. [18].

**Definition 1.** For a given cubic polynomial $ax^3 + bx^2 + cx + d$ in the variable $x$, the discriminant is defined by
$$b^2c^2 - 4ac^3 - 4b^3d - 27a^2d^2 + 18abcd.$$

The number of roots of a 3rd-degree polynomial depends upon the sign of its discriminant: if the latter is positive, then the polynomial has 3 distinct real roots; if it is null, then the polynomial has a real root which is double or triple; last, if it is negative, then it has a unique real root, which is simple. For future use, we denote by $\Delta(a_r)$ the discriminant of the cubic polynomial $P(x)$, considered here as a function of the variable $a_r$. It is easy to see that $\Delta(a_r)$ is a polynomial of degree up to 4 in the variable $a_r$.

We proceed now to the description of the endemic equilibria. As we shall see below, system (4) may admit multiple endemic equilibria. It is thus interesting to check whether these equilibria are endowed with a total ordering, and we begin our investigation by giving a positive answer to this question. This is quite interesting, as it means that the endemic equilibria may be strictly ordered in terms of severity. The proof is postponed to Appendix B.3.

**Theorem 2.** Let $(S_1^*, E_1^*, I_1^*, R_1^*, S_r^*, E_r^*, I_r^*, R_r^*)$ and $(S_1^{**}, E_1^{**}, I_1^{**}, R_1^{**}, S_r^{**}, E_r^{**}, I_r^{**}, R_r^{**})$ be two distinct endemic equilibria of system (4). Assume $S_1^*$ is larger (resp. smaller) than $S_1^{**}$. Then every remaining component of $(S_1^*, E_1^*, I_1^*, R_1^*, S_r^*, E_r^*, I_r^*, R_r^*)$ is smaller (resp. larger) than its analogue in $(S_1^{**}, E_1^{**}, I_1^{**}, R_1^{**}, S_r^{**}, E_r^{**}, I_r^{**}, R_r^{**})$.

Next, the existence and number of endemic equilibria for the cases $R_0 < 1$, $R_0 = 1$, $R_0 > 1$ are addressed respectively in Theorems 3 to 5. The proofs of these results are contained in Appendix B.2.

In the case $R_0 < 1$ (notice that this inequality equivalent to $a_1 < 0$, see formula (5)-(7)), the following theorem holds.

**Theorem 3.** Assume $R_0 < 1$ (that is $a_1 < 0$). Then,

1. if $a_1^2 + a_1c_1b_r - c_r c_1^2 \ge 0$, then system (4) admits no endemic equilibrium;

2. if $a_1^2 + a_1c_1b_r - c_r c_1^2 < 0$, then there exists at least a value at which $\Delta$ vanishes, and the roots of $\Delta$ are necessarily positive. Moreover,

(a) if $a_r > 0$ is smaller than the smallest zero of $\Delta$, then system (4) admits two endemic equilibria,

(b) if $a_r$ is equal to the smallest zero of $\Delta$, then system (4) admits a unique endemic equilibrium,

(c) if $a_r$ is larger the smallest zero of $\Delta$, then system (4) admits no endemic equilibrium.

To provide some insights, Theorem 3 states that, for $R_0 < 1$ and $a_1 < 0$, the system contains endemic equilibria if $a_1^2 + a_1c_1b_r - c_r c_1^2 < 0$ and $a_r$ is less than or equal to the smallest zero of $\Delta$. Here, the first condition appears to be difficult to interpret. However, the second condition indicates the endemic equilibria exist only if $a_r$ is in the interval between 0 and the smallest root of $\Delta$. From the Remark 1, this means that endemic equilibria appear when susceptibility to reinfections is sufficiently high.

For the critical case $R_0 = 1$ (that is $a_1 = 0$), it is again possible to have up to 2 endemic equilibria.

**Theorem 4.** Assume $R_0 = 1$. Then, there exists a value at which $\Delta$ vanishes, and the roots of $\Delta$ are necessarily positive. Moreover,

1. if $a_r$ is smaller than the smallest zero of $\Delta$, system (4) admits a unique endemic equilibrium.

2. if $1 + b_r b_1 - c_r b_1^2 \le 0$,

(a) if $a_r$ is equal to or larger than the smallest zero of $\Delta$, system (4) admits no endemic equilibrium.

3. if $1 + b_r b_1 - c_r b_1^2 > 0$,

(a) if $3b_1c_r - b_r \ge 0$, system (4) admits no endemic equilibrium for $a_r$ equal to or larger than this zero of $\Delta$.

(b) if $3b_1c_r - b_r < 0$ system (4) admits at least a second zero of $\Delta$. Furthermore,
i. if $a_r$ is between the smallest and the second-smallest zero of $\Delta$, system (4) admits two endemic equilibria.

ii. if $a_r$ is equal to the second-smallest zero of $\Delta$, system (4) admits a unique endemic equilibrium.

iii. if $a_r$ is larger than the second zero of $\Delta$, system (4) admits no endemic equilibrium.

Finally, we proceed to the description of the endemic equilibria in case where $R_0 > 1$ (that is $a_1 > 0$).

**Theorem 5.** Assume $R_0 > 1$. Then,

1. if $1 + b_r b_1 - c_r b_1^2 \le 0$, then system (4) admits a unique endemic equilibrium;

2. if $1 + b_r b_1 - c_r b_1^2 > 0$, then $\Delta$ vanishes at two, three, or four distinct values, which are necessarily positive. Moreover,

(a) if $\Delta$ has two distinct zeroes, then system (4) admits a unique endemic equilibrium,

(b) if $\Delta$ has three distinct zeroes, then system (4) admits two endemic equilibria at the intermediate zero and otherwise a unique endemic equilibrium,

(c) if $\Delta$ has four distinct zeroes, then system (4) admits two endemic equilibria at the two intermediate zeroes, three endemic equilibria between the two intermediate zeroes and otherwise a unique endemic equilibrium.

## 3.2 Numerical illustrations

We consider now various sets of parameters verifying either $R_0 < 1$ or $R_0 > 1$ and show numerically that the corresponding numbers of endemic equilibriums coincide with the ones predicted by Theorems 3 and 5. One can also refer to [11] for other examples of coefficients and extensive simulations corresponding to all cases in Theorems 3 to 5. The system of ordinary differential equations was numerically integrated using MATLAB's ode45 solver. The roots of $\mathcal{P}$ were computed within the same environment using the function roots().

### 3.2.1 The case $R_0 < 1$

When $R_0 < 1$, Theorem 3 applies. Consider first the setting where the inequality $a_1^2 + a_1 c_1 b_r - c_r c_1^2 \ge 0$ is verified, then Theorem 3 states the non-existence of endemic equilibrium, regardless of the value of $a_r$. The numerical computation of roots of $\mathcal{P}$ in one such instance in Figure 3a confirms the non-existence of admissible roots generating endemic equilibria, which must be located between the dashed lines in the graph. For information, the polynomial coefficients in this example are

$$a_1 = -0.001, b_1 = 0.81, c_1 = 0.004, b_r = -10, c_r = 1.$$

On the other hand, if the previous inequality is not fulfilled, Theorem 3 predicts instead the appearance of a

<figure><img src="image_4.png" /><figcaption>(a) Case 1 of Theorem 3: 1 endemic equilibrium (b) Case 2 of Theorem 3: up to 2 endemic equilibria</figcaption></figure>

<caption>Figure 3: Example of numerical computations of roots of the polynomial $\mathcal{P}$, $R_0 < 1$</caption>

couple of endemic equilibria for values of $a_r$ sufficiently close to zero. This is again confirmed numerically with Figure 3b, in which we specified

$$a_1 = -0.001, b_1 = 0.81, c_1 = 0.0004, b_r = 15, c_r = 15, \qquad (12)$$
as the parameter values for $\mathcal{P}$. Moreover, an example of biologically coefficients $\gamma_1, \sigma_1, \omega_1, \beta_{1,1}, \beta_{1,r}, \gamma_r, \sigma_r, \omega_r$ and ratio $\frac{\beta_{r,1}}{\beta_{r,r}}$ matching the polynomial coefficients in (12) can be computed (see [11, Proposition 2.5.3] for more details), which is

$$ 
\begin{aligned}
\beta_{1,1} &= 6.15 \times 10^{-3} \text{ days}^{-1}, \quad \beta_{1,r} = 4.98 \times 10^{-3} \text{ days}^{-1}, \quad \omega_1 = 4.13 \times 10^{-5} \text{ days}^{-1}, \\
\omega_r &= 3.11 \times 10^{-4} \text{ days}^{-1}, \quad \sigma_1 = 4.16 \times 10^{-5} \text{ days}^{-1}, \quad \sigma_r = 6.96 \times 10^{-4} \text{ days}^{-1}, \\
\gamma_1 &= 3.90 \times 10^{-3} \text{ days}^{-1}, \quad \gamma_r = 5.80 \times 10^{-4} \text{ days}^{-1}, \quad \frac{\beta_{1,r}}{\beta_{r,r}} = 0.94.
\end{aligned} 
\quad (13) $$

For this set of system parameters, the corresponding bifurcation diagram according to $a_r$ (recall from Remark 1 that $a_r$ is in this case inversely proportional to $\beta_{r,1}, \beta_{r,r}$) is shown in Section 3.2.1. The dashed line here corresponds to an unstable equilibrium, while the solid line represents a stable one. In particular, the diagram illustrates the occurrence of a backward bifurcation [16], i.e., the emergence of both a stable and an unstable endemic equilibrium, as the coefficient $a_r$ decreases (which is equivalent to say that $\beta_{r,1}, \beta_{r,r}$ increases).

<figure><img src="image_5.png" /><figcaption>Figure 4: Bifurcation diagram of (4) fulfilling (13)</figcaption></figure>

### 3.2.2 The case $R_0 > 1$

For $R_0 > 1$, the condition $1 + b_r b_1 - c_r b_1^2 \le 0$ in Theorem 5 guarantees the uniqueness of endemic equilibrium. If this condition is not satisfied, additional steady states may arise as the number of roots of $\Delta$ increases. Although these conditions are difficult to be interpreted biologically, it can be observed that multiple endemic equilibria typically arise when $\Delta$ admits four zeroes and for values of $a_r$ between two intermediate roots of $\Delta$. Hence, it represents a transition between a mode with low susceptibility to reinfections to a higher one.

To facilitate the visualization of the above phenomenon, the roots of $\mathcal{P}$ in an example fulfilling $1 + b_r b_1 - c_r b_1^2 > 0$ is illustrated in Figure 5b. The polynomial coefficients of this example are in fact

$$ a_1 = 0.00025, \quad b_1 = 0.81, \quad c_1 = 0.0004, \quad b_r = 150, \quad c_r = 30. \quad (14) $$

Commenting on the Figure 5b, the red stars in the x-axis represent the zeroes of $\Delta$, of which there are four. This matches indeed the case 2c of Theorem 5. Clearly, the smallest and greatest zero of $\Delta$ are accompanied by an appearance or disappearance of roots of $\mathcal{P}$ outside the set of interest $(a_1, \frac{c_1}{b_1})$ marked by the dashed lines. On the other hand, the intermediate zeroes of $\Delta$ represent an appearance or disappearance of roots inside the set of interest. Notably, Figure 5b illustrates the existence of a swift transition from a low-infectivity endemic equilibrium to a high-infectivity equilibrium as the susceptibility to reinfection increases in some settings of the model (4), (recall that the roots of $\mathcal{P}$ in the admissible region correspond to the values $I_1^*$ of endemic equilibria). This suggests the existence of a reinfection threshold, already underscored by Gomes et al. [14] for simpler reinfection SIRI models.

However, the existence of a reinfection threshold is not always guaranteed, as illustrated numerically by the simulation in Figure 5a corresponding to an example of the case 1 of Theorem 5. The polynomial coefficients in this example are

$$ a_1 = 0.0009, \quad b_1 = 1.5, \quad c_1 = 0.004, \quad b_r = 0.01, \quad c_r = 150. \quad (15) $$

As in the previous subsection, we tested numerically, in Figure 6b, the stability of the endemic equilibria with
<figure><img src="image_6.png" /><figcaption>(a) Case 1 of Theorem 5: 1 endemic equilibrium</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) Case 2 of Theorem 5: up to 3 endemic equilibria</figcaption></figure>

Figure 5: Example of numerical computations of roots of the polynomial $P, R_0 > 1$

<figure><img src="image_8.png" /><figcaption>(a) Biological parameters matching (15): 1 endemic equi- (b) Biological parameters matching (14): up to 3 endemic
librium equilibria</figcaption></figure>

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

Figure 6: Bifurcation diagrams of (4), $R_0 > 1$

two set of biological coefficients that match respectively the polynomial coefficients in (14) and (15). Figure 5a
hints that the unique endemic equilibrium given by the case 1 of Theorem 5 is locally stable. On the other hand,
when multiple endemic equibria arise as in the case 2 of Theorem 5, Figure 5b suggests an unstable intermediate
endemic equilibrium between two locally asymptotically stable endemic equilibria. More specifically, the diagram
shows in fact the occurrence of two saddle node bifurcations. At the first saddle-node bifurcation, a stable and an
unstable equilibrium emerge simultaneously. As $a_r$ increases, the unstable equilibrium moves toward the other
stable one, until a second bifurcation causes them to collide and annihilate each other.

3.2.3 Asymptotic convergence

In this brief subsection, several examples of trajectories are computed around endemic equilibria in the settings
of (12) and (14) such that the systems exhibit respectively 3 and 2 endemic equilibria. The trajectories shown in
Figure 7 illustrate furthermore the stability characteristics of the equilibria. Moreover, every trajectory computed
converges to a steady state and this leads us to conjecture that the $\omega$-limit set of any trajectory of system (4) is
simple, i.e., an equilibrium. This issue will be considered in Section 5 below.

4 Persistence

Disease persistence, i.e., whether the disease eventually dies out after a certain time, is an interesting topic in
epidemiological models. We have seen that the two-stage SEIRS model presents quite complicated dynamics, with
possibly up to three endemic equilibria and one disease-free equilibrium. The multistability and the dimension
<figure><img src="image_10.png" /><figcaption>Figure 7: Temporal evolution of several trajectories around equilibria</figcaption></figure>

of the system render the global analysis of the system difficult and motivate first the study of persistence. For this purpose, let us introduce some notions from [12].

**Definition 2.** Let $X$ be a metric space with metric $d$, and $\phi$ be a continuous flow defined on $X$. Let $E$ be a closed subset of $X$ with boundary $\partial E$ and non-empty interior $\tilde{E}$ nonempty. Assume moreover that $E$ is positively invariant under the flow $\phi$. Denoting $\phi_t(x)$ as the flow map, the flow $\phi$ is called uniformly persistent (regarding the set $E$) if there exists $\varepsilon > 0$ such that for all $x \in \tilde{E}$

$$ \liminf_{t \to \infty} d(\phi_t(x), \partial E) > \varepsilon. $$

Consider then the space $X := \mathbb{R}_+^8$, and the invariant biologically feasible set $\Sigma$ defined by

$$ \Sigma := \{(S_1, E_1, I_1, R_1, S_r, E_r, I_r, R_r) \in X | S_1 + E_1 + I_1 + R_1 + S_r + E_r + I_r + R_r = 1\}. \quad (16) $$

We denote $\partial\Sigma$ the boundary of the set $\Sigma$,

$$ \partial\Sigma_{\max} := \{(S_1, E_1, I_1, R_1, S_r, E_r, I_r, R_r) \in \Sigma | E_1 = 0, I_1 = 0, E_r = 0, I_r = 0\}, \quad (17) $$

one deduce easily with these definitions that $\partial\Sigma_{\max}$ is invariant from the equations (4). Indeed, $\partial\Sigma_{\max}$ is the maximal invariant set in the boundary $\partial\Sigma$.

The following result establishes the uniform persistence of system (4) in the case where $\mathcal{R}_0 > 1$.

**Theorem 6.** System (4) is uniformly persistent in the interior $\tilde{\Sigma}$ of $\Sigma$ defined in (16) if $\mathcal{R}_0 > 1$ and non-persistent if $\mathcal{R}_0 < 1$.

The proof of the theorem is included in Appendix B.4. In particular, the theorem indicates that the disease persists for the SEIRS system (4) when $\mathcal{R}_0 > 1$. Moreover, the proportion of the population bearing the disease is asymptotically bounded from below by a positive constant, uniformly regardless of the initial conditions.

# 5 Global analysis of reinfection models with partial immunity

The exposé above shows that the SEIRS model (4) exhibits multistability. In fact, the system presents up to 2 endemic equilibria for $\mathcal{R}_0 < 1$, 3 endemic equilibria for $\mathcal{R}_0 > 1$ alongside the disease-free equilibrium which always exists. The analysis of asymptotic behavior for systems exhibiting multistability is often tedious. In this section, we derive global asymptotic behavior results for the SIRI and SEIRE models, which include reinfection-induced heterogeneity solely via partial immunity. These models can be considered in fact as subclasses of the general two-stage SEIRS model. The technique employed is inspired from Li and Muldowney's theory, a summary of which can be found in Appendix A.
## 5.1 SIRI model with partial immunity

We examine here the endemic equilibria and the asymptotic behavior of the SIRI model, already mentioned in (3), which takes the assumption that the primary infection leads to a partial immunity to subsequent reinfections. We recall that the model was initially studied in [14]. The asymptotic behavior and stability of the model without demography ($\mu = 0$) has been analyzed in [20, 38]. Taking a step further, we investigate the asymptotic behavior and stability of the SIRI model including demography ($\mu > 0$). This addition is essential in terms of modelling, given that the study of reinfections requires a framework with longer time-scale, and therefore demography effects are not negligible in shaping the disease dynamics.

To get started, let us consider first the two-stage SIS model

$$ 
\begin{aligned}
\dot{S}_1 &= \mu - (\beta_{1,1}I_1 + \beta_{1,r}I_r + \mu)S_1, \\
\dot{I}_1 &= (\beta_{1,1}I_1 + \beta_{1,r}I_r)S_1 - (\gamma_1 + \mu)I_1, \\
\dot{S}_r &= \gamma_1I_1 + \gamma_rI_r - (\beta_{r,1}I_1 + \beta_{r,r}I_r + \mu)S_r, \\
\dot{I}_r &= (\beta_{r,1}I_1 + \beta_{r,r}I_r)S_r - (\gamma_r + \mu)I_r.
\end{aligned} 
\quad (18) $$

This system is in fact a simplification of the two-stage SEIRS model and can be formally obtained by approaching $\sigma_1, \sigma_r, \omega_1, \omega_r$ to $+\infty$ in (4). One may check that the proofs of Theorems 2 to 6 can be directly adapted to the SIS model (18) (see [11, Section 2.6]), albeit with a simpler definition of $a_1, b_1, c_1, a_r, b_r$ and $\mathcal{R}_0$.

$$ 
\begin{aligned}
a_1 &= \frac{\mu}{\beta_{1,1}}(\mathcal{R}_0 - 1), & b_1 &= \frac{\beta_{1,r}}{\beta_{1,1}}, & c_1 &= \frac{\beta_{1,r}}{\beta_{1,1}}\frac{\mu}{\gamma_1 + \mu}, \\
a_r &= \frac{\gamma_r + \mu}{\beta_{r,r}}, & b_r &= -\frac{\beta_{r,1}}{\beta_{r,r}} + \frac{\gamma_1}{\mu}, & c_r &= \frac{\gamma_1}{\mu}\frac{\beta_{r,1}}{\beta_{r,r}}, & \mathcal{R}_0 &= \frac{\beta_{1,1}}{\gamma_1 + \mu}.
\end{aligned} 
\quad (19) $$

The SIRI model with demography introduced in [14] can then be viewed as a special case of system (18). To see this, we set $\beta_{1,1} = \beta_{1,r} = \beta, \beta_{r,1} = \beta_{r,r} = \psi\beta, \gamma_1 = \gamma_r = \gamma$. This results in the following two-stage SIS system:

$$ 
\begin{aligned}
\dot{S}_1 &= \mu - \beta(I_1 + I_r)S_1 - \mu S_1, \\
\dot{I}_1 &= \beta(I_1 + I_r)S_1 - (\gamma + \mu)I_1, \\
\dot{S}_r &= \gamma(I_1 + I_r) - \psi\beta(I_1 + I_r)S_r - \mu S_r, \\
\dot{I}_r &= \psi\beta(I_1 + I_r)S_r - (\gamma + \mu)I_r.
\end{aligned} 
\quad (20) $$

Introducing in addition the state variables

$$ I := I_1 + I_r, \quad S := S_1, \quad R := S_r, $$

we obtain the same SIRI system with partial immunity and demography as in [13].

$$ 
\begin{aligned}
\dot{S} &= \mu - \beta IS - \mu S, \\
\dot{I} &= \beta IS + \psi\beta IR - (\gamma + \mu)I, \\
\dot{R} &= \gamma I - \psi\beta IR - \mu R.
\end{aligned} 
\quad (21) $$

We describe in the following theorem the endemic equilibria of the systems (20) and (21), we apply Theorem 5 and Theorem 3 respectively to the case $\mathcal{R}_0 > 1$ and $\mathcal{R}_0 < 1$. For $\mathcal{R}_0 > 1$, this implies the uniqueness of endemic equilibrium. On the other hand, the matter is more intricate for $\mathcal{R}_0 < 1$. The proof of the result is postponed to Appendix B.5.

**Theorem 7.** The systems (20) and (21) admit

1. a unique endemic equilibrium if $\mathcal{R}_0 > 1$, that is $\beta > \gamma_1 + \mu$,

2. no endemic equilibrium if $\mathcal{R}_0 < 1$ and $\beta \le \mu$.

Moreover, if $\mathcal{R}_0 < 1$ and $\beta > \mu$, that is $\mu < \beta < \gamma_1 + \mu$, then there exists a threshold such that the systems (20) and (21) admit

1. no endemic equilibrium if $\psi < \psi^*$,

2. a unique endemic equilibrium if $\psi = \psi^*$,

3. two endemic equilibrium if $\psi > \psi^*$.

**Remark 2.** For a disease ensuring partial immunity after infection, one has $\psi < 1$. Consequently, in such case, systems (20) and (21) do not exhibit endemic equilibrium when $\mathcal{R}_0 < 1$.
On the disease persistence, the uniform persistence still holds when $\mathcal{R}_0 > 1$ for both models. The proof of this statement can be adapted straightforwardly from the one of Theorem 6 or found in [11, Subsection 2.6.2].

**Theorem 8.** System (20) (resp. system (21)) is uniformly persistent in the interior of $\{x \in \mathbb{R}_+^4 | S_1 + I_1 + S_r + I_r = 1\}$ (resp. $\{x \in \mathbb{R}_+^3 | S + I + R = 1\}$) if $\mathcal{R}_0 > 1$ and non-persistent if $\mathcal{R}_0 < 1$.

Subsequently, we turn our attention to the asymptotic convergence of the solutions. In the final result of this subsection, we prove that any trajectory converges to an equilibrium for the reduced system (21), assuming that $\mathcal{R}_0 \neq 1$. In particular, this implies the same for the two-stage SIS system (20). The details of the proof can be found in Appendix B.6.

**Theorem 9.** Assume $\mathcal{R}_0 \neq 1$, any solution $x(t)$ of system (20) (resp. system (21)) with initial condition $x(0) \in \{x \in \mathbb{R}_+^4 | S_1 + I_1 + S_r + I_r = 1\}$ (resp. $x(0) \in \{x \in \mathbb{R}_+^3 | S + I + R = 1\}$) converges to an equilibrium.

In particular, thanks to the uniform persistence, the unique endemic equilibrium of both systems are GAS in the interior of the respective invariant sets when $\mathcal{R}_0 > 1$.

## 5.2 SEIRE model with partial immunity

The previous techniques can also be applied to a more sophisticated class of reinfection model compared to (20) and (21), by adding another class E modelling the disease incubation period. Again, the model we intend to study can also be considered as a special case of the two-stage SEIRS model. Indeed, by setting in (4)

$$ \beta_{1,1} = \beta_{1,r} = \beta, \quad \beta_{r,1} = \beta_{r,r} = \psi\beta, \quad \sigma_1 = \sigma_r = \sigma, \quad \gamma_1 = \gamma_r = \gamma, \quad \omega_1 = \omega_r = +\infty, $$

we obtain the following two-stage SEIS system:

$$ \begin{aligned} \dot{S}_1 &= \mu - (\beta(I_1 + I_r) + \mu)S_1, \\ \dot{E}_1 &= \beta(I_1 + I_r)S_1 - (\sigma + \mu)E_1, \\ \dot{I}_1 &= \sigma E_1 - (\gamma + \mu)I_1, \\ \dot{S}_r &= \gamma(I_1 + I_r) - (\psi\beta(I_1 + I_r) + \mu)S_r, \\ \dot{E}_r &= \psi\beta(I_1 + I_r)S_r - (\sigma + \mu)E_r, \\ \dot{I}_r &= \sigma E_r - (\gamma + \mu)I_r. \end{aligned} \tag{22} $$

Now, let us introduce the state variables

$$ S := S_1, \quad E := E_1 + E_R, \quad I := I_1 + I_r, \quad R := S_r. $$

This leads to a reduced SEIRE system with partial immunity and demography:

$$ \begin{aligned} \dot{S} &= \mu - \beta SI - \mu S, \\ \dot{E} &= \beta SI + \psi\beta IR - (\sigma + \mu)E, \\ \dot{I} &= \sigma E - (\gamma + \mu)I, \\ \dot{R} &= \gamma I - \psi\beta IR - \mu R. \end{aligned} \tag{23} $$

It can be checked straightforwardly that the results of Theorems 2 to 6 still hold for system (22), albeit with $A_r = \frac{\sigma}{\sigma+\mu} \frac{\gamma}{\gamma+\mu}$ and $D_r = 1$ instead of the definition in (5). Moreover in this case,

$$ 1 - A_r = 1 - \frac{\sigma\gamma}{\sigma\gamma + \mu(\sigma + \gamma) + \mu^2} = \frac{\mu(\sigma + \gamma + \mu)}{\sigma\gamma + \mu(\sigma + \gamma + \mu)}, $$

the polynomial coefficients $a_1, b_1, c_1, a_r, b_r, c_r$ are

$$ \begin{aligned} a_1 &= \mu \frac{\beta\sigma - (\gamma + \mu)(\sigma + \mu)}{\beta(\gamma + \mu)(\sigma + \mu)}, \quad b_1 = 1, \quad c_1 = \frac{\mu}{\gamma + \mu} \frac{\sigma}{\sigma + \mu}, \\ a_r &= \frac{\sigma\gamma + \mu(\sigma + \gamma + \mu)}{\psi\beta(\sigma + \gamma + \mu)}, \quad b_r = \frac{\sigma\gamma - (\sigma + \gamma + \mu)\mu}{\mu(\sigma + \gamma + \mu)}, \quad c_r = \frac{\sigma\gamma}{\mu(\sigma + \gamma + \mu)}, \end{aligned} $$

and the basic reproduction number is expressed by

$$ \mathcal{R}_0 = \frac{\beta}{\gamma + \mu} \frac{\sigma}{\sigma + \mu}. $$

Similarly to systems (20) and (21), the endemic equilibrium is unique for $\mathcal{R}_0 > 1$ and the description is more complicated for $\mathcal{R}_0 < 1$. The proof of the next result can be found in Appendix B.5.
**Theorem 10.** The systems (22) and (23) admit

1. a unique endemic equilibrium if $\mathcal{R}_0 > 1$, that is $\beta > \frac{(\gamma + \mu)(\sigma + \mu)}{\sigma}$;

2. no endemic equilibrium if $\mathcal{R}_0 < 1$ and

$$ \beta \le \frac{\mu(\sigma + \gamma + \mu)}{\sigma}. \qquad (24) $$

Moreover, if $\mathcal{R}_0 < 1$ and (24) is not fulfilled, that is

$$ \frac{\mu(\sigma + \gamma + \mu)}{\sigma} < \beta < \frac{(\gamma + \mu)(\sigma + \mu)}{\sigma}, $$

then there exists a threshold $\psi^* > 1$ and systems (22) and (23) admit

1. no endemic equilibrium if $\psi < \psi^*$,

2. a unique endemic equilibrium if $\psi = \psi^*$,

3. two endemic equilibria if $\psi > \psi^*$.

**Remark 3.** Again, for disease ensuring partial immunity after infection, i.e., $\psi < 1$, systems (22) and (23) cannot exhibit endemic equilibria when $\mathcal{R}_0 < 1$ (see Remark 2).

We also enunciate a persistence result, which will be needed for demonstrating the asymptotic convergence, and whose proof can be again easily adapted from the one of Theorem 6.

**Theorem 11.** System (22) (resp. system (23)) is uniformly persistent in the interior of $\{x \in \mathbb{R}_+^6 | S_1 + E_1 + I_1 + S_r + E_r + I_r = 1\}$ (resp. $\{x \in \mathbb{R}_+^4 | S + E + I + R = 1\}$) if $\mathcal{R}_0 > 1$ and non-persistent if $\mathcal{R}_0 < 1$.

We have shown that systems (22) and (23) present again one unique endemic equilibrium when $\mathcal{R}_0 > 1$ and up to two endemic equilibria when $\mathcal{R}_0 < 1$. Excluding the critical case $\mathcal{R}_0 = 1$, the next result shows that any solution of (23) converges either to the disease-free equilibrium or to an endemic equilibrium.

**Theorem 12.** Assume $\mathcal{R}_0 \ne 1$, any solution $x(t)$ of (22) (resp. of (23)) with initial condition $x(0) \in \{x \in \mathbb{R}_+^6 | S_1 + E_1 + I_1 + S_r + E_r + I_r = 1\}$ (resp. $x(0) \in \{x \in \mathbb{R}_+^4 | S + E + I + R = 1\}$) converges to an equilibrium.

The proof of Theorem 12 is contained in Appendix B.6. Notice also that, combined with the persistence result of Theorem 11, the unique endemic equilibrium is again GAS in the interior of the invariant set when $\mathcal{R}_0 > 1$.

# 6 Discussion

We gather in this section some comments on the results, classified by themes.

**Multistability in concrete diseases** We first examine whether multistability is a plausible phenomenon in real-world scenarios. We found no data-based estimates of the parameters of the second stage of model (4) for real life diseases. Nevertheless, we attempted to apply with incomplete data the model previously studied to a real disease example. For this application, we draw inspiration from the following parameter estimates for COVID-19, borrowed from [9] and fitted to the SEIRS model (1):

$$ \beta = 1.29 \text{ days}^{-1}, \quad \sigma = \frac{1}{4.5} \text{ days}^{-1}, \quad \gamma = \frac{1}{3.1} \text{ days}^{-1}, \quad \omega = \frac{1}{365} \text{ days}^{-1}, \quad \mu = \frac{1}{75} \text{ years} = \frac{1}{27,375} \text{ days}^{-1}, $$

yielding $\mathcal{R}_0 = 3.998$. Consequently, we put

$$ \sigma_1 = \sigma_r = \frac{1}{4.5} \text{ days}^{-1}, \quad \gamma_1 = \gamma_r = \frac{1}{3.1} \text{ days}^{-1}, \quad \omega_1 = \omega_r = \frac{1}{365} \text{ days}^{-1}, \quad \mu = \frac{1}{27,375} \text{ days}^{-1}. \quad (25) $$

In addition, in every test below, we set

$$ \beta_{1,1} = \beta_{1,r}, \quad \beta_{r,1} = \rho\beta_{r,r}, $$

for different values of $\rho > 0$. Therefore, in the following simulations, the only remaining parameters are $\beta_{1,1} = \beta_{1,r}$, $\rho := \frac{\beta_{r,1}}{\beta_{r,r}}$ and $\beta_{r,1}$, the latter being the bifurcation parameter of the diagrams.

• We first display in Figure 8a the bifurcation diagram for model (4) corresponding to $\rho = 1$, i.e., $\beta_{r,1} = \beta_{r,r}$, with the parameters mentioned above in (25) and

$$ \beta_{1,1} = \beta_{1,r} = 1.29 \text{ days}^{-1}. $$
<figure><img src="image_11.png" /><figcaption>(a) Case $R_0 > 1$, $\beta_{1,1} = \beta_{1,r} = 1.29 \text{ days}^{-1}$, $\rho = 1$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Case $R_0 > 1$, $\beta_{1,1} = \beta_{1,r} = 1.29 \text{ days}^{-1}$, $\rho = 10$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) Case $R_0 > 1$, $\beta_{1,1} = \beta_{1,r} = 0.33 \text{ days}^{-1}$, $\rho = 0.8$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) Case $R_0 < 1$, $\beta_{1,1} = \beta_{1,r} = 0.29 \text{ days}^{-1}$, $\rho = 1$</figcaption></figure>

<caption>Figure 8: Bifurcation diagrams of (4)</caption>

In this case, $R_0 = 3.998$. Clearly, the diagram shows a configuration with unique endemic equilibrium and the existence of a reinfection threshold with a rapid increase in endemicity as $\beta_{r,1} = \beta_{r,r}$ increase to $\approx 0.5 \text{ days}^{-1}$. This indicates that it is essential to study the reinfection threshold for diseases such as COVID-19, for example in order to reduce the endemicity of the disease.

• We now consider the case where $\rho = 10$. The value of the basic-offspring number is unchanged. The bifurcation diagram in Figure 8b shows again an unique equilibrium. However the transition pace for the reinfection threshold here is much slower compared to the previous case in Figure 8a.

• As a next example, set

$$ \beta_{1,1} = \beta_{1,r} = 0.33 \text{ days}^{-1}, \quad \rho = 0.8. $$

In this case, the basic reproduction number $R_0 = 1.023$ is still greater than 1 and the corresponding bifurcation diagram in Figure 8c exhibits now a transition phase with coexistence of three endemic equilibria while the system transits from low to a higher level of endemicity.

• Last, we take the parameters

$$ \beta_{1,1} = \beta_{1,r} = 0.29 \text{ days}^{-1} $$

in order to derive a case where $R_0 = 0.899 < 1$. The corresponding diagram is shown in Figure 8d, which underlines the presence of a backward bifurcation, a phenomenon which has been predicted in Theorem 3. For reference, the bifurcation happens for $\beta_{r,1} = \beta_{r,r} \approx 0.4 \text{ days}^{-1}$.

To conclude, we were able to reproduce the complex multistable behaviours exhibited theoretically in the paper with parameters for the two-stage SEIRS model (4) chosen realistically close from values measured for an example of real-life disease fitted to the classical SEIRS model (1). However, this was only achievable with a significant restriction on the choice of the contact rates $\beta_{1,1}, \beta_{1,r}$: the latter must be chosen so that $R_0$ is close to one and $\beta_{1,r} \le \beta_{1,1}$. This possibility demonstrates the potential interest of this study for real-life epidemic.
**Attempt at a comprehensive interpretation** As the previous paragraph revealed the possibility of occurrence of multiple equilibrium points for realistic parameter values, we attempt here to understand the conditions under which this phenomenon appears.

* As a general rule, from the numerical explorations presented above, it seems that multiple endemic equilibria (either a backward bifurcation when $R_0 < 1$ or a transitional regime with 3 endemic equilibria when $R_0 > 1$) arise when the basic reproduction number is close to the critical value $R_0 = 1$.

* Moreover, in the simulations, we notice that in the configurations where multiple endemic equilibria arise, at least one of the contact rates $\beta_{r,1}, \beta_{r,r}$, of the second stage is larger than those of the first stage, $\beta_{1,1}, \beta_{1,r}$. This is clear and has been mathematically established for the simplified systems (21), (23), for which it is necessary to have $\psi > 1$ in order to obtain multiple endemic equilibria, specifically a backward bifurcation. Whether this necessity remains true for the general system is an open and interesting question. While numerical simulations suggest that the tendency persists, a rigorous mathematical characterization of the parameter regimes leading to multiple endemic equilibria in the general model remains elusive.

* To provide biological plausibility for scenarios for which the contact rates related to reinfection could exceed those of the primary infection, notice that several diseases may lead to weaker immunity after first infection, in particular for those presenting antibody-dependent enhancement [22] as HIV or dengue fever. Apart from such biological reasons, reinfection-induced higher susceptibility and/or infectiousness may be caused by behavioural changes, e.g. loosening of the distancing rules following recovery from primary infection.

* Commenting now more specifically on the phenomenon of backward bifurcation arising when $R_0 < 1$, Theorem 3 states that the inequality $a_1^2 + a_1c_1b_r - c_r c_1^2 \ge 0$ precludes the existence of backward bifurcation and ensures that no endemic equilibrium may arise. This condition can therefore serve as a basis for control in order to eliminate endemicity. Even though the inequality itself is difficult to interpret in the general case, in the simplified SIRI and SEIRE models (21) and (23), it leads to an upper bound for the contact rate $\beta$ (see Theorem 7 and Theorem 10), which strongly suggests intervention in the control of $\beta$ in order to achieve the objective. Traditional methods of disease control (vaccination campaign, quarantine...) usually aim to reduce the effective basic reproduction number to $R_0 < 1$. However in the presence of backward bifurcation in our models, this objective may be insufficient to eliminate the disease: an additional reduction of $\beta$ may be required to achieve this goal. A detailed analysis of the consequences on vaccination strategies of the backward bifurcation is out of the scope of this paper.

* On the other hand, for $R_0 > 1$, the existence of reinfection threshold, i.e., the rapid transition from a low-endemicity equilibria to a higher one while the susceptibility of the reinfected increases, is also confirmed in the numerical simulations with realistic parameters. Furthermore, the reinfection threshold may be accompanied, in some configurations such as in Figure 8c, by saddle node bifurcations and the co-existence of three endemic equilibria simultaneously. The existence of reinfection threshold highlights the importance of intervening on the control of the transmission parameters of reinfected individuals, for example with a vaccination campaign targeting also the seropositive, as an effective way to reduce endemicity.

* In Figure 8c, we notice that the parameter configuration leading to the emergence of 3 endemic equilibria is characterized by $R_0 \ge 1, \beta_{1,1} \ge \beta_{1,r}$ and the inequality

$$ \frac{\beta_{r,1}}{\beta_{r,r}} < 1. $$

The fact that the inequality above may yield multiple endemic equilibria has an intuitive interpretation. Specifically, this assumption suggests that reinfections are predominantly caused by contact with other reinfected individuals, rather than with primary infectives. In such a scenario, if the initial number of reinfected individuals is negligible, reinfections cannot be sustained, and the dynamics are largely governed by the first stage of the model. As a result, the system effectively behaves like a SEIR model without reinfection. Conversely, if the population contains a sufficiently large proportion of reinfected individuals, the second stage becomes more influential in shaping the dynamics. This shift gives rise to a qualitatively distinct behavior, and thus to the possibility of multiple endemic equilibria.

* The notion of a reinfection threshold introduced by Gomes et al. [14] has sparked a notable controversy with Breban et al. [5], primarily due to the absence of any qualitative change in the dynamics of the original SIRI model at the proposed threshold. In contrast, our model can exhibit genuine qualitative bifurcations characterized by the emergence of three endemic equilibria, a phenomenon typically associated with the transition from a low-endemicity to a high-endemicity regime. Interestingly, the intuition behind the reinfection threshold in [14] lies in the observation that, as the reinfection parameter increases, the intermediate SIRI model shifts from approximating the classical SIR model (without reinfection) to resembling the SIS reinfection model. This perspective shares similarity with our interpretation of the heuristic condition $\frac{\beta_{r,1}}{\beta_{r,r}} < 1$ as a means of inducing saddle node bifurcations with the appearance of three coexisting endemic equilibria in our framework.
**Multistability and asymptotic convergence in the general case** As a final comment, we would like to revisit some mathematical aspects related to multistability. This paper has established several new results concerning the asymptotic convergence of the SIRI (21) and SEIRE (23) models, both of which are specific instances of the general two-stage system (4), by leveraging the theory developed by Li and Muldowney. A natural next step would be to extend this analysis to the full system (4) using the same theoretical framework. However, two major obstacles prevent such a generalization:

* First, application of Li & Muldowney involves the study of the second additive compound matrix of the Jacobian of the non-linear system (see Appendix A). For a $n$-dimensional system, the associated second additive compound matrix is of size $\binom{n}{2}$, which can rapidly become intractable for the analysis as $n$ increases. For instance, the compound matrix to be analyzed for the general two-stage model (4) is of size $28 \times 28$.

* Even if we set aside the practical difficulties involved in analyzing the associated compound matrix, the characteristics of the general model (4) theoretically rules out the applicability of the Li & Muldowney approach. Indeed, as shown in Theorem 5, the system may admit up to three endemic equilibria in the interior $\bar{\Sigma}$ of the domain $\Sigma$ defined in (16). Since $\bar{\Sigma}$ is positively invariant (by persistence arguments from Section 4), the analysis is restricted to this subset. We affirm then that the presence of three distinct equilibria within the invariant set directly obstructs the application of Li & Muldowney's geometric criterion. This becomes clear when considering recent generalizations of the theory via the concept of $k$-contraction [43, 44], which provide an useful geometrical insight on the theory. In fact, at its core, the Li & Muldowney's approach aim to establish that the system is 2-contractive, from which global convergence can be inferred. A system is said to be $k$-contractive on an invariant set $D \subset \mathbb{R}^n$ if all $k$-dimensional parallelo topes in $D$ contract exponentially under the flow [43]. In the case $k = 2$, this implies that the flow must contract the area of any parallelogram. Now, as the invariant set $\bar{\Sigma}$ contains at least three equilibria, these can be used to form the vertices of a parallelogram. Since equilibria are fixed points and remain stationary under the flow, therefore the area of such a parallelogram cannot decrease. This simple geometrical argument demonstrates that 2-contraction cannot hold in the general case, thereby precluding the use of Li & Muldowney theory.

Nevertheless, the successful application in this paper to the SIRI and SEIRE models demonstrates that the approach provided by Li & Muldowney theory can be particularly relevant for classes of models exhibiting backward bifurcations.

# 7 Conclusion

In this paper, we proposed and analyzed an eight dimensional two-stage SEIRS model describing the behavior of diseases with reinfections. The second stage of this model considers a reinfection process that behaves differently from the infection process in the first stage. The model possesses a unique disease-free equilibrium, which is shown to be locally asymptotically stable when the threshold condition $\mathcal{R}_0 < 1$ is satisfied, where $\mathcal{R}_0$ is the basic reproduction number of the model. The endemic equilibria of the system have been comprehensively described, revealing that it is possible to have up to 2 (resp. 3) endemic equilibria when $\mathcal{R}_0 \le 1$ (resp. $\mathcal{R}_0 > 1$). Numerical simulations confirmed the possibility for the models to exhibit backward bifurcations for $\mathcal{R}_0 < 1$ and saddle node bifurcation for $\mathcal{R}_0 > 1$. Then, we established the uniform persistence for the two-stage model, within the framework of [12]. The results were applied afterwards to two specific SIRI and SEIRE models intervening in the modelling of diseases with partial immunity, which are subclasses of the general two-stage model. We found that the endemic equilibrium is unique for these models when $\mathcal{R}_0 > 1$, and simple necessary and sufficient conditions for the appearance of two endemic equilibria have been derived for the case $\mathcal{R}_0 < 1$. We also established global asymptotic convergence to one of their steady states for the SIRI and SEIRE models, even in settings where two endemic equilibria coexist with the disease-free equilibrium. The core of the proof is based on a generalized convergence criterion from [30], within the framework of Li and Muldowney theory [27]. Finally, we offered several numerical examples with realistic parameters, showing that the complex equilibrium behaviors exhibited may appear in real-world scenarios, and the concrete epidemiological implications of our findings were also discussed.

# Acknowledgements

The financial support of French National Research Agency (ANR) is acknowledged (Project NOCIME, ANR-23-CE48-0004).

# References

[1] R. Anguelov, S. M. Garba, and S. Usaini. Backward bifurcation analysis of epidemiological model with partial immunity. *Computers & Mathematics with Applications*, 68(9):931–940, 2014. BIOMATH 2013.
[2] M. M. Ballyk, C. C. McCluskey, and G. S. K. Wolkowicz. Global analysis of competition for perfectly substitutable resources with linear response. *Journal of Mathematical Biology*, 51(4):458–490, 2005.

[3] O. N. Björnstad, K. Shea, M. Krzywinski, and N. Altman. The SEIRS model for infectious disease dynamics. *Nature Methods*, 17, 6 2020.

[4] O. Boldea, A. Alipoor, S. Pei, J. Shaman, and G. Rozhnova. Age-specific transmission dynamics of SARS-CoV-2 during the first 2 years of the pandemic. *PNAS Nexus*, 3(2), 01 2024.

[5] R. Breban and S. Blower. The reinfection threshold does not exist. *Journal of Theoretical Biology*, 235(2):151–152, 2005.

[6] G. Butler and P. Waltman. Persistence in dynamical systems. *Journal of Differential Equations*, 63(2):255–263, 1986.

[7] H.-D. Chiang, M. Hirsch, and F. Wu. Stability regions of nonlinear autonomous dynamical systems. *IEEE Transactions on Automatic Control*, 33(1):16–27, 1988.

[8] W. Coppel. *Stability and Asymptotic Behavior of Differential Equations*. Heath mathematical monographs. Heath, 1965.

[9] T. Crellen, L. Pi, E. Davis, T. Pollington, T. Lucas, D. Ayabina, A. Borlase, J. Turner, K. Prem, G. Medley, P. Klepac, and T. Hollingsworth. Dynamics of sars-cov-2 with waning immunity in the uk population. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 376:20200274, 05 2021.

[10] M. Fan, M. Y. Li, and K. Wang. Global stability of an SEIS epidemic model with recruitment and a varying total population size. *Mathematical Biosciences*, 170(2):199–208, 2001.

[11] M. Fang. *Mathematical modelling, observation and identification of epidemiological models with reinfection*. PhD thesis, Sorbonne Université, 2024.

[12] H. I. Freedman, S. Ruan, and M. Tang. Uniform persistence and flows near a closed positively invariant set. *Journal of Dynamics and Differential Equations*, 6:583–600, 1994.

[13] M. C. Gomes, A. Margheri, and C. Rebelo. Stability and persistence in a compartment model of pulmonary tuberculosis. *Nonlinear Analysis: Theory, Methods & Applications*, 48(4):617–636, 2002.

[14] M. G. M. Gomes, L. J. White, and G. F. Medley. Infection, reinfection, and vaccination under suboptimal immune protection: epidemiological perspectives. *Journal of Theoretical Biology*, 228, 6 2004.

[15] D. Greenhalgh, O. Diekmann, and M. C. de Jong. Subcritical endemic steady states in mathematical models for animal infections with incomplete immunity. *Mathematical Biosciences*, 165(1):1–25, 2000.

[16] A. Gumel. Causes of backward bifurcations in some epidemiological models. *Journal of Mathematical Analysis and Applications*, 395(1):355–365, 2012.

[17] H. W. Hethcote. Simulations of pertussis epidemiology in the United States: effects of adult booster vaccinations. *Mathematical Biosciences*, 158(1):47–73, 1999.

[18] R. Irving. *Integers, Polynomials, and Rings: A Course in Algebra*. Undergraduate Texts in Mathematics. Springer New York, 2003.

[19] P. Kaklamanos, A. Pugliese, M. Sensi, and S. Sottile. A geometric analysis of the SIRS model with secondary infections. *SIAM Journal on Applied Mathematics*, 84(2):661–686, 2024.

[20] G. Katriel. Epidemics with partial immunity to reinfection. *Mathematical Biosciences*, 228, 12 2010.

[21] M. J. Keeling and P. Rohani. *Modeling Infectious Diseases in Humans and Animals*. Princeton University Press, 2008.

[22] R. Kulkarni. *Antibody-Dependent Enhancement of Viral Infections*, pages 9–41. Springer Singapore, Singapore, 2020.

[23] M. Li and J. Muldowney. On Bendixson's Criterion. *Journal of Differential Equations*, 106(1):27–39, 1993.

[24] M. Li, J. Muldowney, and P. Driessche. Global stability of SEIRS models in epidemiology. *The Canadian Applied Mathematics Quarterly*, 7, 01 1999.

[25] M. Y. Li and J. S. Muldowney. Global stability for the SEIR model in epidemiology. *Mathematical Biosciences*, 125, 2 1995.

[26] M. Y. Li and J. S. Muldowney. On R.A. Smith's Autonomous Convergence Theorem. *Rocky Mountain Journal of Mathematics*, 25(1):365–378, 1995.

[27] M. Y. Li and J. S. Muldowney. A geometric approach to global-stability problems. *SIAM Journal on Mathematical Analysis*, 27, 7 1996.
[28] M. Y. Li and J. S. Muldowney. Dynamics of differential equations on invariant manifolds. *Journal of Differential Equations*, 168(2):295-320, 2000.

[29] W. Lohmiller and J.-J. E. Slotine. On contraction analysis for non-linear systems. *Automatica*, 34(6):683-696, 1998.

[30] G. Lu and Z. Lu. Geometric approach to global asymptotic stability for the SEIRS models in epidemiology. *Nonlinear Analysis: Real World Applications*, 36:20-43, 2017.

[31] G. Lu and Z. Lu. Global asymptotic stability for the SEIRS models with varying total population size. *Mathematical Biosciences*, 296, 2 2018.

[32] A. Margheri and C. Rebelo. Some examples of persistence in epidemiological models. *Journal of mathematical biology*, 46(6):564, 2003.

[33] M. Martcheva. *An Introduction to Mathematical Epidemiology*. Texts in Applied Mathematics. Springer US, 2015.

[34] C. C. McCluskey. A strategy for constructing Lyapunov functions for non-autonomous linear differential equations. *Linear Algebra and its Applications*, 409:100-110, 2005. Special Issue in honor of Pauline van den Driessche.

[35] J. S. Muldowney. Compound matrices and ordinary differential equations. *The Rocky Mountain Journal of Mathematics*, pages 857-872, 1990.

[36] J. S. Muldowney. Compound matrices and applications. *Universidad de los Andes, Bogotá, Colombia*, 1998.

[37] A. V. Okhuse et al. Estimation of the probability of reinfection with COVID-19 by the susceptible-exposed-infectious-removed-undetectable-susceptible model. *JMIR Public Health and Surveillance*, 6(2):e19097, 2020.

[38] R. Pagliara, B. Dey, and N. E. Leonard. Bistability and resurgent epidemics in reinfection models. *IEEE Control Systems Letters*, 2:290-295, 2018.

[39] C. C. Pugh. The Closing Lemma. *American Journal of Mathematics*, 89(4):956-1009, 1967.

[40] H. L. Smith and H. R. Thieme. *Dynamical Systems and Population Persistence*. Graduate studies in mathematics. American Mathematical Society, 2011.

[41] H. R. Thieme. Uniform persistence and permanence for non-autonomous semiflows in population biology. *Mathematical Biosciences*, 166(2):173-201, 2000.

[42] P. van den Driessche and J. Watmough. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. *Mathematical Biosciences*, 180, 11 2002.

[43] C. Wu, I. Kanevskiy, and M. Margaliot. *k*-contraction: Theory and applications. *Automatica*, 136:110048, 2022.

[44] S. Zoboli, A. Cecilia, U. Serres, D. Astolfi, and V. Andrieu. LMI conditions for *k*-contraction analysis: A step towards design. In *2023 62nd IEEE Conference on Decision and Control (CDC)*, pages 1903-1910, 2023.

## A Elements of Li and Muldowney theory

For non-planar systems exhibiting multistability, the asymptotic behavior analysis is often difficult. One possible solution is to adapt the technique developed by Li and Muldowney [27], well-known for demonstrating global stability of a unique endemic equilibrium in an absorbing set. Indeed, the same theory may be used to prove more generally that any trajectory in an invariant set of a given system converges to an equilibrium, which is not necessarily unique.

We present a quick review of Li and Muldowney autonomous convergence result. Let $f$ be a $C^1$ function on a domain $D$ and consider the non-linear system

$$ \dot{x} = f(x), \qquad (A.1) $$

let us denote $\delta x_1(t), \delta x_2(t), \dots, \delta x_k(t) \in \mathbb{R}^n$ as $k$ time-varying vectors representing virtual displacements which are solutions of the variational equation of system (A.1):

$$ \dot{\delta x} = \frac{df(x)}{dx} \delta x. \qquad (A.2) $$

We may represent the $k$-dimensional parallelotope generated by the vectors $\delta x_1(t), \delta x_2(t), \dots, \delta x_k(t)$ as an exterior product $z(t) := \wedge_{1 \le i \le k} \delta x_i(t)$ evolving in $\mathbb{R}^{(n)}$, where $(n)$ is the binomial coefficient. With this definition, $z(t)$ is a solution of the $k$-th compound equation below associated to (A.2) (see e.g. [28, Section 2]).

$$ \dot{z} = \frac{df}{dx}^{[k]}(x)z, \qquad (A.3) $$
where, for any matrix A, we denote $A^{[k]}$ as the k-th additive compound matrix of A. For an introduction to exterior product and compound matrices, one refers to [28] and [35, 36]. With these definitions, we introduce the convergence result proved by Li and Muldowney in [27] and reformulated by Ballyk et al. [2] who relaxed the need for a unique equilibrium. The statement of the following theorem is slightly adapted from [2, Theorem 6.2].

**Theorem A.1.** Let $k = 2$, $D$ be a simply connected open subset of $\mathbb{R}^n$ such that solutions of equation (A.1) with $x(0) = x_0 \in D$ remain in $D$ for all finite time. Let $\|\cdot\|$ be a norm on $\mathbb{R}^{(n_2)}$ and let $Q$ be an $(n_2) \times (n_2)$ matrix-valued function on $D$ such that $\|Q^{-1}\|$ is bounded on $D$. If there exist real numbers $T, g > 0$ such that for all solutions of (A.1) in $D$ and all $t \ge T$, any solution $z(t)$ of (A.3) satisfies

$$ \|Q(x(t))z(t)\| \le \|Q(x(0))z(0)\|e^{-gt}, \quad (A.4) $$

then every solution of (A.1) in $D$ bounded away from the boundary converges to an equilibrium.

In other words, if the origin of the 2-nd compound equation (A.3) is exponentially asymptotically stable uniformly w.r.t. a state dependent norm $\|Q(x)z\|$ for any trajectory of (A.1) in the set $D$, then the $\omega$-limit set of any solution of (A.1) in $D$ not approaching its boundary is simple, i.e., is reduced to an equilibrium point.

To give a brief idea of the proof of Li and Muldowney's result in [27], the condition (A.4) is first proved to be a Bendixson criterion [23] that precludes the existence of periodic solution for system (A.1). Indeed, let $D = \mathbb{R}^n$ for simplicity, the quantity $\|Q(x)z\|$ gives a measure of the volume of a $k$-parallelotope $z$ in $\mathbb{R}^n$. Hence, the condition (A.4) for $k = 2$ translates the fact that the area of any surface in $\mathbb{R}^n$ shrinks exponentially under the flow of (A.1), and this rules out periodic solutions. In order to obtain the convergence to an equilibrium, the criterion is moreover demonstrated to be robust under $C^1$ perturbations of the vector field $f$. The Pugh's closing Lemma [39] is finally invoked to deduce the non-existence of non-constant non-wandering points, hence any non-wandering point is an equilibrium.

Practically, the Bendixson criterion (A.4) is demonstrated through the use of Lyapunov-like function. In fact, a carefully selected vector norm $|\cdot|$ on $\mathbb{R}^{(n_2)}$ is employed as Lyapunov function to prove the property (A.4). In this case, it is interesting to introduce the Lozinskii measure $\mu(E)$ of a matrix $E$ with respect to the norm $|\cdot|$, defined by

$$ \mu(E) = \lim_{h \to 0^+} \frac{|\mathbf{I} + hE| - 1}{h}, $$

which is used to estimate the rate of expansion $-g$ in (A.4). Indeed, it is possible to verify that (see [34])

$$ D_+|z(t)| \le \mu \left( \frac{df}{dx}^{[2]}(x(t)) \right) |z(t)|, $$

where $D_+$ is the right-hand derivative. More generally, for $V(x, z) := |Q(x)z|$, where $Q(\cdot)$ is a $C^1$ non-singular $(n_2) \times (n_2)$ matrix-valued function

$$ D_+|Q(x(t))z(t)| \le \mu \left( Q_f(x(t))Q^{-1}(x(t)) + Q(x(t)) \frac{df}{dx}^{[2]}(x(t))Q^{-1}(x(t)) \right) |Q(x(t))z(t)|, $$

where $Q_f$ is the directional derivative of the matrix $Q$ in the direction of the vector field $f$. For more details about Lozinskii measure, one may refer to [8].

$$ \text{Denoting } B(x) := Q_f(x)Q^{-1}(x) + Q(x) \frac{df}{dx}^{[2]}(x)Q^{-1}(x), \text{ if} \\ \sup_{x \in D} \mu(B(x)) < 0, \quad (A.5) $$

or more generally,

$$ \limsup_{t \to \infty} \sup_{x_0 \in D} \frac{1}{t} \int_0^t \mu(B(x(s, x_0))ds < 0, \quad (A.6) $$

then the property (A.4) is fulfilled. In fact, the conditions (A.5) and (A.6) are the original Bendixson criterions proposed by Li and Muldowney [26, 27]. Insofar as they are used in the proof of results in Section 5, it is also worth noting that the Bendixson criterion (A.6) can be further relaxed, see for example [30, Theorem 2.4] or [31, Theorems 2.2, 2.3 and 2.4].

As epidemiological models are often considered in invariant and biologically feasible sets of co-dimension 1 (generated by the fact that the total population is constant along the time), a useful extension of the convergence result is available for differential equations evolving on invariant manifolds [28].

Indeed, if the system (A.1) possesses a simply connected invariant manifold $\Gamma := \{x \in \mathbb{R}^n_+ | g(x) = 0\}$, where $g: \mathbb{R}^n \to \mathbb{R}^m$ is $C^2$ and $\dim(\frac{\partial g}{\partial x}) = m$, then there is a function $N: \mathbb{R}^n \to \mathbb{R}^m$ [28] such that

$$ \frac{\partial g}{\partial x}(x)f(x) = N(x)g(x), $$
and the convergence result of Li and Muldowney holds as well if we replace the 2-nd compound equation (A.3) by the $(m+2)$-th compound equation (A.7) [28]:

$$ \dot{z} = \left( \frac{df}{dx} \right)^{[m+2]} (x) - \nu(x) \mathbf{I} z, \quad (A.7) $$

where

$$ \nu(x) := \mathrm{Tr}(N(x)), \quad (A.8) $$

is the trace of $N(x)$. By the same arguments then as in [2, Theorem 6.2], we deduce a result similar to Theorem A.1 in the context of systems evolving in an invariant manifold.

**Theorem A.2.** D be a simply connected open subset of $\Gamma$ such that solutions of equation (A.1) with $x(0) = x_0 \in D$ remain in D for all finite time. Let $\|\cdot\|$ be a norm on $\mathbb{R}^{(m+2)}$ and let $Q$ be a $(\binom{n}{m+2}) \times (\binom{n}{m+2})$ matrix-valued function on D such that $\|Q^{-1}\|$ is bounded on D. If there exist reals $T, g > 0$ such that for all solutions of (A.1) in D and all $t \ge T$, any solution $z(t)$ of (A.3) satisfies

$$ \|Q(x(t))z(t)\| \le \|Q(x(0))z(0)\| e^{-gt}, \quad (A.9) $$

then every solution of (A.1) in D bounded away from the boundary converges to an equilibrium.

The interest of studying the associated compound (A.7) instead of (A.3) lies in the fact that the $(m+2)$-th additive compound matrix may be in some cases smaller in size than the 2-nd additive compound matrix. This is essential because compound matrices explode rapidly in size as the dimension $n$ increases.

Finally, if the same assumptions in Theorems A.1 and A.2 hold in a compact invariant set $\Sigma$ instead, then [28, Theorem 6.1] states that any trajectory in $\Sigma$ converges to an equilibrium.

It is noteworthy that the convergence criterion (A.4) has been introduced recently in the control community as a property called *k-contraction*, (see e.g. [44, Definition 1]), with $k = 2$ in the context of Li and Muldowney theory. Indeed, the property of *k-contraction* can be viewed as a generalization of the incremental stability property of *contraction* [29] and constitutes a novel field of active research [43].

# B Proofs

## B.1 Proof of Proposition 1

Let $(S_1^*, E_1^*, I_1^*, R_1^*, S_r^*, E_r^*, I_r^*, R_r^*)$ be an equilibrium of system (4), then the following equations hold

$$
\begin{aligned}
0 &= \mu - (\beta_{1,1}I_1^* + \beta_{1,r}I_r^* + \mu)S_1^*, \\
0 &= (\beta_{1,1}I_1^* + \beta_{1,r}I_r^*)S_1^* - (\sigma_1 + \mu)E_1^*, \\
0 &= \sigma_1 E_1^* - (\gamma_1 + \mu)I_1^*, \\
0 &= \gamma_1 I_1^* - (\omega_1 + \mu)R_1^*, \\
0 &= \omega_1 R_1^* + \omega_r R_r^* - (\beta_{r,1}I_1^* + \beta_{r,r}I_r^* + \mu)S_r^*, \\
0 &= (\beta_{r,1}I_1^* + \beta_{r,r}I_r^*)S_r^* - (\sigma_r + \mu)E_r^*, \\
0 &= \sigma_r E_r^* - (\gamma_r + \mu)I_r^*, \\
0 &= \gamma_r I_r^* - (\omega_r + \mu)R_r^*.
\end{aligned}
\quad (B.1)
$$

In other words,

$$
\begin{aligned}
R_r^* &= \frac{\gamma_r}{\omega_r + \mu} I_r^*, & E_r^* &= \frac{\gamma_r + \mu}{\sigma_r} I_r^*, \\
E_r^* &= \frac{\beta_{r,1}I_1^* + \beta_{r,r}I_r^*}{\sigma_r + \mu} S_r^*, & S_r^* &= \frac{\omega_1 R_1^* + \omega_r R_r^*}{\beta_{r,1}I_1^* + \beta_{r,r}I_r^* + \mu}, \\
R_1^* &= \frac{\gamma_1}{\omega_1 + \mu} I_1^*, & E_1^* &= \frac{\gamma_1 + \mu}{\sigma_1} I_1^*, \\
E_1^* &= \frac{\beta_{1,1}I_1^* + \beta_{1,r}I_r^*}{\sigma_1 + \mu} S_1^*, & S_1^* &= \frac{\mu}{\beta_{1,1}I_1^* + \beta_{1,r}I_r^* + \mu},
\end{aligned}
\quad (B.2)
$$

After elimination of $S_1^*, E_1^*, R_1^*, S_r^*, E_r^*, R_r^*$, the remaining quantities $I_1^*$ and $I_r^*$ are related by

$$
\begin{aligned}
I_1^* &= \frac{\sigma_1}{\sigma_1 + \mu} \frac{\mu}{\gamma_1 + \mu} \frac{\beta_{1,1}I_1^* + \beta_{1,r}I_r^*}{\beta_{1,1}I_1^* + \beta_{1,r}I_r^* + \mu}, & (B.3) \\
I_r^* &= \frac{\beta_{r,1}I_1^* + \beta_{r,r}I_r^*}{\beta_{r,1}I_1^* + \beta_{r,r}I_r^* + \mu} \frac{\sigma_r}{\sigma_r + \mu} \frac{\gamma_r}{\gamma_r + \mu} \frac{\omega_r}{\omega_r + \mu} \left( \frac{\omega_1 \omega_r + \mu \gamma_1}{\omega_r \omega_1 + \mu \gamma_r} I_1^* + I_r^* \right).
\end{aligned}
$$
All the coefficients being positive, any nonnegative nonzero solution $(I_1^*, I_r^*)$ of (B.3) yields a nonnegative nonzero equilibrium value $(S_1^*, E_1^*, I_1^*, R_1^*, S_r^*, E_r^*, I_r^*, R_r^*)$ such that $S_1^*+E_1^*+I_1^*+R_1^*+S_r^*+E_r^*+I_r^*+R_r^* = 1$. The endemic equilibria are therefore in one-to-one correspondence with the nonnegative nonzero solutions of system (B.3). With the coefficients $A_1, B_1, C_1, A_r, B_r, C_r, D_r$ defined in (5), it is easy to check that equations (B.3) can be rewritten as (6).

Denoting now $x := I_1^*, y := I_r^*$, we retrieve the equations (8) using the coefficients $(a_1, b_1, c_1, a_r, b_r, c_r)$ defined in (7). Moreover, for any $y > 0$,

$$y = \Psi_1(x) := \frac{1}{2} \left( b_r x - a_r + \sqrt{(b_r x - a_r)^2 + 4c_r x^2} \right), \qquad (B.4)$$

is the unique solution of the second equation (the other root is negative). On the other hand, the expression of $y$ given by the first identity in (8) must be

$$y = \Psi_r(x) := x \frac{x - a_1}{c_1 - b_1 x}. \qquad (B.5)$$

Notice that,

$$\max\{0, a_1\} < \frac{c_1}{b_1} < 1.$$

Therefore, $\Psi_r(x) > 0$ with $x > 0$ if and only if

$$x \in (0, 1) \cap (a_1, \frac{c_1}{b_1}), \qquad (B.6)$$

and the endemic equilibria of system (4) are in one-to-one correspondence with the roots of the equation

$$\Psi_1(x) = \Psi_r(x), \quad x \in (0, 1) \cap (a_1, \frac{c_1}{b_1}),$$

for $\Psi_1, \Psi_r$ given in (B.4)-(B.5). The problem may be further transformed by isolating the square root in $\Psi_1(x)$ and squaring both sides of the equality, hence the equality above yields

$$2x \frac{x - a_1}{c_1 - b_1 x} - (b_r x - a_r) = \sqrt{(b_r x - a_r)^2 + 4c_r x^2}, \qquad (B.7)$$

thus

$$\left(2x \frac{x - a_1}{c_1 - b_1 x} - (b_r x - a_r)\right)^2 = (b_r x - a_r)^2 + 4c_r x^2,$$

that is

$$4x^2 \left( \frac{x - a_1}{c_1 - b_1 x} \right)^2 - 4x \frac{x - a_1}{c_1 - b_1 x} (b_r x - a_r) = 4c_r x^2.$$

The following relation is therefore verified for any solution of (B.8) fulfilling (B.6),

$$x \frac{x - a_1}{c_1 - b_1 x} \left( x \frac{x - a_1}{c_1 - b_1 x} - (b_r x - a_r) \right) = c_r x^2 > 0. \qquad (B.8)$$

As $x \frac{x-a_1}{c_1-b_1x} > 0$ in the interval given in (B.6), the previous inequality implies

$$b_r x - a_r < x \frac{x - a_1}{c_1 - b_1 x} < 2x \frac{x - a_1}{c_1 - b_1 x},$$

and the left-hand side of (B.7) is indeed positive for any solution of (B.8) in the interval (B.6). Therefore, the problem is reduced to finding the solutions of polynomial equation (B.8) instead of (B.7).

Last, we multiply by $\frac{(c_1-b_1x)^2}{x}$ on both sides of (B.8) to obtain the cubic polynomial equation.

$$(x - a_1)(x(x - a_1) - (c_1 - b_1x)(b_r x - a_r)) = c_r x(c_1 - b_1x)^2,$$

which is equivalent to $P(x) = 0$, with $P$ as defined in the Proposition 1. $\square$

## B.2 Proofs of Theorems 3 to 5

The proofs of Theorems 3 to 5 are based on the following observations, whose proofs are straightforward.

**Lemma B.1.** The polynomial $P$ fullfills

$$P(a_1) = -c_r a_1 (c_1 - b_1 a_1)^2, \quad P(0) = -a_1 c_1 a_r, \quad P\left(\frac{c_1}{b_1}\right) = \frac{c_1}{b_1} \left(\frac{c_1}{b_1} - a_1\right)^2 > 0,$$

and the values $P(a_1), P(0)$ are thus both positive iff $a_1 < 0$, that is $R_0 < 1$; and both negative iff $a_1 > 0$, that is $R_0 > 1$.
**Lemma B.2.** The polynomials P and Q fulfill

$$P(a_1) = -c_r(c_1 - b_1a_1)^2 < 0, \quad P\left(\frac{c_1}{b_1}\right) = \left(\frac{c_1}{b_1} - a_1\right)^2 > 0, \qquad (B.9a)$$

$$Q(a_1) = Q\left(\frac{c_1}{b_1}\right) = 0, \quad Q > 0 \text{ on } \left(a_1, \frac{c_1}{b_1}\right), \qquad (B.9b)$$

$$P(0) = a_1^2 + a_1 c_1 b_r - c_r c_1^2, \quad \lim_{x \to +\infty} \frac{P(x)}{x^2} = 1 + b_r b_1 - c_r b_1^2. \qquad (B.9c)$$

**Proof of Theorem 3.** For $\mathcal{R}_0 < 1$, notice that $a_1 < 0$, and the endemic equilibria are in one-to-one correspondence with the zeros of $\mathcal{P}$ in $(0, \frac{c_1}{b_1})$.

First, let us consider the zeroes of the polynomial P. From (B.9a), we deduce that P has a zero in $(a_1, \frac{c_1}{b_1})$ and there can be no other zero in this interval, as P is of degree 1 or 2. If it exists, the second zero of P is then either in $(-\infty, a_1)$ (and necessarily P is convex) or in $(\frac{c_1}{b_1}, +\infty)$ (and necessarily P is concave). Moreover, the zero of P located in $(a_1, \frac{c_1}{b_1})$ may be either in $(a_1, 0)$ (in which case necessarily $P(0) = a_1^2 + a_1 c_1 b_r - c_r c_1^2 > 0$), in $(0, \frac{c_1}{b_1})$ (in which case $P(0) < 0$) or at 0 (in which case $P(0) = 0$).

Consider now these three cases separately. Notice that the derivative of $xP(x)$ at $x = 0$ is equal to $P(0)$. Therefore,

* if $P(0) = a_1^2 + a_1 c_1 b_r - c_r c_1^2 > 0$, then $xP(x)$ is strictly increasing at $x = 0$ and in addition $P(x) > 0$ in $(0, \frac{c_1}{b_1})$. As $Q(x) > 0$ in this interval, one concludes that, for any $a_r > 0$, $\mathcal{P}(x) = xP(x) + a_r Q(x) > 0$ in $(0, \frac{c_1}{b_1})$, and there is no endemic equilibrium. This demonstrates point 1 in the case of a strict inequality.

* if $P(0) = a_1^2 + a_1 c_1 b_r - c_r c_1^2 = 0$, 0 is a root of the polynomial P. Thus, due to (B.9a), $P(x) > 0$ in $(0, \frac{c_1}{b_1})$. The same argument shows that, for any $a_r > 0$, $\mathcal{P}(x) > 0$ in $(0, \frac{c_1}{b_1})$: no endemic equilibrium exists. This, together with the previous point, fully demonstrates case 1.

* if $P(0) = a_1^2 + a_1 c_1 b_r - c_r c_1^2 < 0$ (case 2), consider first the situation where $1 + b_r b_1 - c_r b_1^2 \neq 0$ and P is a 2nd-degree polynomial. In this case, $xP(x)$, $\mathcal{P}(x)$ are 3rd-degree polynomials and, for some $x^* \in (0, \frac{c_1}{b_1})$, one has

$$P(x) < 0 \text{ in } (0, x^*), \quad P(x^*) = 0, \quad P(x) > 0 \text{ in } \left(x^*, \frac{c_1}{b_1}\right).$$

Therefore, the 3rd-order polynomial $xP(x)$ vanishes at $x = 0$ and $x = x^*$, is negative in $(0, x^*)$ and positive in $(x^*, \frac{c_1}{b_1})$. As $Q > 0$ in $[0, x^*] \subset (a_1, \frac{c_1}{b_1})$, then, for $a_r$ close to zero, increasing the value of $a_r$ leads to the displacement of the first root in $(0, \frac{c_1}{b_1})$ of $\mathcal{P}(x) = xP(x) + a_r Q(x)$ to the right and the second to the left, and $\mathcal{P}$ still admits two distinct roots in $(0, x^*)$.

Moreover, $xP(x) > 0$ at $x = a_1 < 0$ and at $x = \frac{c_1}{b_1}$ thanks to (B.9a). The polynomial $xP(x)$ being of 3rd degree, it goes to infinity when $|x| \to +\infty$, with different signs when $x \to +\infty$ and $x \to -\infty$. Therefore, one of the two limits is negative, and the third (real) zero of $xP(x)$ is either in $(-\infty, a_1)$, or in $(\frac{c_1}{b_1}, +\infty)$. Let us study now these two cases in more details.

* If $xP(x)$ possesses a zero in $(-\infty, a_1)$, then the latter is unique and $xP(x)$ diverges to $-\infty$ when $x \to -\infty$, while taking positive values in $(\frac{c_1}{b_1}, +\infty)$. Recall in addition that the polynomial $Q$ is of degree 2, vanishes at $x = a_1$ and is negative in $(-\infty, a_1)$, so that $\mathcal{P}(a_1) > 0$, $\lim_{x \to -\infty} \mathcal{P}(x) = -\infty$ and, for any positive value $a_r$, $\mathcal{P}(x) = xP(x) + a_r Q(x)$ has at least a zero in $(-\infty, a_1)$. For $a_r > 0$ close to zero, it is already described that $\mathcal{P}$ possesses two roots in $(0, x^*) \subset (0, \frac{c_1}{b_1})$. When the parameter $a_r$ increases, the same situation holds, until the two roots located in $(0, x^*)$ merge together and disappear, as $a_r Q(x)$ (which is positive in $((0, \frac{c_1}{b_1}))$) became not cancellable by the term $xP(x)$ in this interval. This happens when $a_r$ reaches a critical value $a_r^*$, where $a_r^*$ is such that the associated polynomial $\mathcal{P}(x)$ admits a double root in the interval $(0, \frac{c_1}{b_1})$, that is $\Delta = 0$ for $a_r = a_r^*$.

Notice that a pair of zeros of $\mathcal{P}$ appears in $(\frac{c_1}{b_1}, +\infty)$ for large enough values of $a_r$. One has to wonder whether this happens for a value of $a_r$ smaller than $a_r^*$. For such a value, $\mathcal{P}$ would have one root in $(-\infty, a_1)$, two distinct roots in $(0, \frac{c_1}{b_1})$ and (at least) one in $(\frac{c_1}{b_1}, +\infty)$. Such a configuration with 4 roots is impossible for a 3rd-degree polynomial, and $a_r^*$ is thus the smallest $a_r > 0$ for which $\Delta$ vanishes.
○ On the contrary, if $xP(x)$ possesses a zero on $(\frac{c_1}{b_1}, +\infty)$, then $xP(x)$ diverges to $-\infty$ when $x \to +\infty$ and a similar argument applies to deduce that $\mathcal{P}$ has at least one zero in $(\frac{c_1}{b_1}, +\infty)$. Therefore, $\mathcal{P}$ cannot have a third zero in $(0, \frac{c_1}{b_1}) \subset (a_1, \frac{c_1}{b_1})$. Again, a critical value $a_r^*$ exists at which the corresponding $\Delta$ vanishes. For $a_r$ greater than $a_r^*$, $\mathcal{P}$ has no root in $(a_1, \frac{c_1}{b_1})$. In parallel, a pair of zeros of $\mathcal{P}$ appears in $(-\infty, a_1)$ for large enough values of $a_r$, but one shows with the same argument on the polynomial degree, that this occurs only after the disappearance of the two roots in $(0, \frac{c_1}{b_1})$. Again, the critical value $a_r^*$ is the smallest $a_r > 0$ for which $\Delta = 0$.

○ Last, let us consider the degenerate case where $P(0) = a_1^2 + a_1c_1b_r - c_r c_1^2 < 0$, but $1 + b_r b_1 - c_r b_1^2 = 0$ and the degree of $P$ is 1; The polynomial $P$ is an increasing linear function that cuts the x-axis at a point $x^* \in (0, \frac{c_1}{b_1})$. The function $xP(x)$ is positive in $(-\infty, 0)$ and $(x^*, +\infty)$, negative in $(0, x^*)$. Thus, for $a_r = 0$, $\mathcal{P}$ admits a root at $x = 0$ and $x = x^* \in (0, \frac{c_1}{b_1})$ and no roots elsewhere as its degree 2 precludes the appearance of new roots. Moreover, the polynomial $Q$ is positive on $(a_1, \frac{c_1}{b_1})$. Therefore, for $a_r > 0$ close to 0, an increase in $a_r$ leads to the displacement of the first root of $\mathcal{P}$ to the right and the second to the left. This goes on until the roots merge and disappear at a critical value $a_r^*$, beyond which no root exists in this interval. Again, these are precisely the situations depicted in the cases 2a, 2b and 2c.

Finally, let us prove that any root of $\Delta$ must be positive in this case. Assuming that $a_r < 0$, as $xP(x)$ cancels out at least once in $[0, \frac{c_1}{b_1}) \subset (a_1, \frac{c_1}{b_1})$, is positive at the extremities of $(a_1, \frac{c_1}{b_1})$ and moreover $a_r Q(x)$ is negative in $(a_1, \frac{c_1}{b_1})$ with its roots situated at both extremities of the same interval, this implies that $xP(x) + a_r Q(x)$ is positive at $x = a_1$ and $x = \frac{c_1}{b_1}$, and negative at $x = 0$, thus it vanishes twice in $(a_1, \frac{c_1}{b_1})$. Additionally, $\mathcal{P}$ has a third root in $(-\infty, a_1) \cup (\frac{c_1}{b_1}, +\infty)$. This means that $\mathcal{P}$ has always three real roots and $\Delta$ does not cancel out when $a_r < 0$ and $a_1^2 + a_1 c_1 b_r - c_r c_1^2 < 0$.

The previous considerations thus prove that when $P(0) < 0$, the behavior always obeys the analysis provided by the three cases 2a to 2c.

This completes the proof of Theorem 3. $\square$

**Proof of Theorem 4.** When $\mathcal{R}_0 = 1$, $a_1 = 0$, and the endemic equilibria are in one-to-one correspondence with the zeros of $\mathcal{P}$ in $(0, \frac{c_1}{b_1})$.

The polynomial $P$ has a zero in $(0, \frac{c_1}{b_1})$ due to (B.9a). Additionally, as $P$ is of degree 1 or 2, there cannot be another zero in this interval. If it exists, the second zero of $P$ is thus either in $(-\infty, 0)$ (and necessarily $P$ is convex) or in $(\frac{c_1}{b_1}, +\infty)$ (and necessarily $P$ is concave). Hence, $xP(x)$ has a root at $x = 0$, another one in $(0, \frac{c_1}{b_1})$ and the last one is either in $(-\infty, 0)$ or $(\frac{c_1}{b_1}, +\infty)$. Anyway, because $Q(0) = Q(a_1) = 0$ the polynomial $\mathcal{P}$ has always a zero at $x = 0$ when $\mathcal{R}_0 = 1$. The analysis further divides into two cases:

* If $P$ is strictly concave (which is equivalent to $\frac{P(x)}{x^2} = 1 + b_r b_1 - c_r b_1^2 < 0$), the second root of $P$ is in $(\frac{c_1}{b_1}, +\infty)$, which leads $xP(x)$ and thus $\mathcal{P}$ to diverge to $-\infty$ when $x \to +\infty$. As $\mathcal{P}(\frac{c_1}{b_1}) > 0$, there is always a root in $(\frac{c_1}{b_1}, +\infty)$ and consequently the polynomial $\mathcal{P}$ has always three real roots for every value of $a_r \ge 0$. Now, assume that $a_r = 0$, the third root of $\mathcal{P}$ is in $(0, \frac{c_1}{b_1})$ and is equal to the root of $P$ in the same interval. Recalling again that $Q$ is positive in $(0, \frac{c_1}{b_1})$ and negative outside of it, then, as $a_r$ is increasing, $\mathcal{P}$ increases in $(0, \frac{c_1}{b_1})$ and thus the third root of $\mathcal{P}$ moves to the left and leaves the interval $(0, \frac{c_1}{b_1})$ when reaching $x = 0$. It happens when $a_r$ reaches the smallest positive root of $\Delta$.

* If $P$ is a linear function (i.e. $1 + b_r b_1 - c_r b_1^2 = 0$), it is thus an increasing function as its root is in $(0, \frac{c_1}{b_1})$ and that $P(\frac{c_1}{b_1}) > 0$. $\mathcal{P}$ is then a convex function with roots at $x = 0$ and in $(0, \frac{c_1}{b_1})$ when $a_r = 0$. As $a_r$ is increasing, the later root decreases and leaves the interval $(0, \frac{c_1}{b_1})$ by forming a double root at $x = 0$ when $a_r$ reaches the first positive root of $\Delta$. Thus, endemic equilibrium of the system disappears when $a_r \ge a_r^*$.

* If $P$ is strictly convex (i.e. $1 + b_r b_1 - c_r b_1^2 > 0$) the second root is in $(-\infty, 0)$, $\mathcal{P}(x)$ vanishes now once respectively at $x = 0$, in $(-\infty, 0)$ and in $(0, \frac{c_1}{b_1})$ when $a_r = 0$. The function $xP(x)$ is positive between the first zero and $x = 0$, negative between $x = 0$ and the third zero, and again positive between the third zero and $+\infty$. Moreover, an increase of $a_r > 0$ leads the function $\mathcal{P}$ to decrease in $(-\infty, a_1)$ and increase in $(a_1, \frac{c_1}{b_1})$. Thus, the first root of $\mathcal{P}$ is displacing to the right and the third is decreasing as $a_r$ increases before merging and vanish. As the two roots have opposite sign for $a_r = 0$, this implies that one of them must reach the third root $x = 0$ in order to merge and form a double root with each other. This happens when $a_r$ reaches the smallest positive root $a_r^*$ of $\Delta$. Therefore, we have three different scenarios:

- If it is the positive root that reaches $x = 0$ before becoming negative, then the system loses its unique endemic equilibrium when $a_r > a_r^*$.
- It the two roots converge exactly at the point $x = 0$, then they also vanish at the same point and, as previously, the system loses its endemic equilibrium when $a_r > a_r^*$.

- If it is the negative root that passes through $x = 0$ and becomes positive, then the system gains a second endemic equilibrium before that the two come together and then vanish when $a_r$ reaches the second positive root of $\Delta$.

We may further determine which root is reaching zero at the first critical value $a_r^*$ by computing the value of $P''(0)$ when $a_r = a_r^*$.

- If $P''(0) > 0$, $P$ has a double root at $x = 0$ and $P$ is positive in a neighborhood of $x = 0$, then thanks to $\lim_{x \to -\infty} P(x) = -\infty$ (due to the convexity of $P$), the third root is in $(-\infty, 0)$ and it is the initially positive root that is passing through zero.

- If $P''(0) < 0$, $P$ has a double root at $x = 0$ and is negative in a neighborhood of $x = 0$, by the same argument we have $\lim_{x \to +\infty} P(x) = +\infty$ and then the third root is positive. Hence, it is the initially negative root that is passing through zero.

- If $P''(0) = 0$, $P$ has a triple root at $x = 0$, then both the positive and negative roots are reaching zero together.

Let us compute $P''(0)$:

$$ 
\begin{aligned}
P' &= P(x) + xP'(x) + a_r^*Q'(x), \\
P'' &= 2P'(x) + xP''(x) + a_r^*Q''(x), \\
P''(0) &= 2P'(0) + a_r^*Q''(0) \\
&= 2(-2a_1 - b_r c_1 - a_1 b_1 b_r + 2c_r b_1 c_1) - 2b_1 a_r^* = 2(-b_r c_1 + 2c_r b_1 c_1) - 2b_1 a_r^*.
\end{aligned}
$$

Thanks to the fact that $0 = P'(0) = P(0) + a_r^*Q'(0)$, one has

$$ a_r^* = -\frac{P(0)}{Q'(0)} = \frac{-c_r c_1^2}{c_1 + a_1 b_1} = -c_r c_1. $$

Thus,

$$ P''(0) = 2(3c_r b_1 c_1 - b_r c_1) = 2c_1(3c_r b_1 - b_r), $$

and which scenario occurs ultimately depends on the sign of $3c_r b_1 - b_r$.

Finally, let us show that any root of $\Delta$ must be positive. Notice in any case that, for a negative value of $a_r$, $a_r Q(x)$ is negative in $(0, \frac{c_1}{b_1})$ and its roots are located at the extremities of the interval. On the other hand, as $P(0) < 0$, $xP(x)$ is negative on the right side of $x = 0$ and positive on the left side of $x = 0$ and at $x = \frac{c_1}{b_1}$ (see (B.9a)). Hence, $P$ is negative on the right side of $x = 0$ and positive at $x = \frac{c_1}{b_1}$, and $P$ has at least two distinct roots, one at the point $x = 0$ and the other one in $(0, \frac{c_1}{b_1})$. Moreover, if $P$ is of degree 3, there must be a third root in the set $(-\infty, 0) \cup (\frac{c_1}{b_1}, +\infty)$, as $P > 0$ at $x = \frac{c_1}{b_1}$ and on the right side of the point $x = 0$. Therefore, the corresponding $\Delta$ of $P$ cannot vanish for a non-positive value of $a_r$.

All cases being examined, this concludes the proof of Theorem 4. $\square$

**Proof of Theorem 5.** When $R_0 > 1$, then $a_1 > 0$, and the endemic equilibria are in one-to-one correspondence with the zeros of $P$ in $(a_1, \frac{c_1}{b_1})$.

As in the proof of Theorem 3, the polynomial $P$ has a zero in $(a_1, \frac{c_1}{b_1})$, and due to (B.9a) and the fact that it is of degree 1 or 2, there can be no other zero in this interval. If the degree of $P$ is 2, then the second zero of $P$ is either in $(-\infty, a_1)$ (and necessarily $P$ is convex) or in $(\frac{c_1}{b_1}, +\infty)$ (and necessarily $P$ is concave).

Let us study the behavior of $xP(x)$, and afterward of $P(x) = xP(x) + a_rQ(x)$.

* If $\lim_{x \to +\infty} \frac{P(x)}{x^2} = 1 + b_r b_1 - c_r b_1^2 < 0$, then $P$ is concave, and its second zero is in $(\frac{c_1}{b_1}, +\infty)$. The 3rd-degree polynomial $xP(x)$ possesses one zero at $x = 0 \in (-\infty, a_1)$, one in $(a_1, \frac{c_1}{b_1})$, and one in $(\frac{c_1}{b_1}, +\infty)$. Moreover, it converges to $+\infty$ (resp. $-\infty$) when $x \to -\infty$ (resp. $+\infty$). The polynomial $Q$ being negative in each interval $(-\infty, a_1)$ and $(\frac{c_1}{b_1}, +\infty)$, positive in $(a_1, \frac{c_1}{b_1})$ and null at $a_1$ and $\frac{c_1}{b_1}$, the polynomial $P(x)$ thus possesses, for any $a_r > 0$, (at least) one zero in each of these three intervals. Applying Proposition 1, one deduces existence of a unique endemic equilibrium of the system (4).
* If $\lim_{x \to +\infty} \frac{P(x)}{x^2} = 1 + b_r b_1 - c_r b_1^2 = 0$. In this case, $P$ is a polynomial of degree 1, which is increasing due to (B.9a). The 2nd-degree polynomial $xP(x)$ then possesses two zeroes, one at $x = 0 \in (-\infty, a_1)$ and one in $(a_1, \frac{c_1}{b_1})$, and is negative between the two roots, and positive outside. Arguing as in the previous case, one deduces that $P(x)$ possesses, for any $a_r > 0$, (at least) one zero in $(a_1, \frac{c_1}{b_1})$. This zero is indeed unique, as otherwise the 2nd-degree polynomial $P(x)$ would have three zeroes. This, together with the previous considerations, establishes the case 1 of the statement.

* If $\lim_{x \to +\infty} \frac{P(x)}{x^2} = 1 + b_r b_1 - c_r b_1^2 > 0$ (case 2), then $P$ is convex, and its second zero is inside $(-\infty, a_1)$: for some $x^* < a_1$ and $x^{**} \in (a_1, \frac{c_1}{b_1})$, one has

$$P(x) < 0 \text{ in } (x^*, x^{**}), \quad P(x^*) = P(x^{**}) = 0, \quad P(x) > 0 \text{ in } (x^{**}, \frac{c_1}{b_1}),$$

and similarly for $xP(x)$ except at $x = 0$ where $xP(x) = 0$. Recall in addition that $Q(x) > 0$ in $(a_1, \frac{c_1}{b_1})$ and $Q(x) < 0$ in $\mathbb{R} \setminus [a_1, \frac{c_1}{b_1}]$. Therefore, for $a_r = 0$ or a sufficiently small positive value of $a_r > 0$, $P(x) = xP(x) + a_r Q(x)$ possesses, similarly to $xP(x)$, two zeroes in $(-\infty, a_1)$, and exactly one in $(a_1, \frac{c_1}{b_1})$ (recall $P$ is a 3rd-degree polynomial). This situation goes on while $a_r$ increases, until the two roots located in $(-\infty, a_1)$ collapse and disappear, at a point where $\Delta$ necessarily vanishes.

On the other hand, for sufficiently large values of $a_r > 0$, we see that $P(x)$ possesses this time two zeroes in $(\frac{c_1}{b_1}, +\infty)$, as $P(\frac{c_1}{b_1}) > 0$, $\lim_{x \to +\infty} P(x) > 0$ and $Q(x) < 0$ for $x > \frac{c_1}{b_1}$. Invoking again the fact that $P$ cannot have more than 3 roots, one sees that this occurs through the appearance of a double root of $P$ as the parameter $a_r$ increases and reaches a new zero of the equation $\Delta = 0$, and afterward splits into two distinct zeroes that subsequently persist for arbitrary large values of $a_r > 0$ greater than this zero of $\Delta$. Moreover, the polynomial $\Delta$ being of degree 4, it may still have, on top of the two previous real zeroes exhibit, up to two distinct zeroes.

* First, if there is no supplementary real zero of $\Delta$, then from the observation above the polynomial $P$ has always a unique zero in $(a_1, \frac{c_1}{b_1})$ for any $a_r > 0$, and by Proposition 1 there exists a unique endemic equilibrium. This situation is covered by the case 2a in the statement.

* Second, assume two distinct real zeroes exist, these zeroes cannot be larger than the second critical value exhibited. As they are also larger than the first one, they are thus sandwiched between the two critical values, indicating the apparition of two supplementary roots of $P$ for any value of $a_r$ between these two zeroes of $\Delta$ (case 2c). Fixing $x$ and varying $a_r$, the quantity $P(x)$ is an increasing function of $a_r$ when $x \in (-\infty, a_1)$ and decreasing one when $x \in (\frac{c_1}{b_1}, +\infty)$. Thus $P(x)$ is positive for $x \in (-\infty, a_1)$ and $a_r$ bigger than the smallest positive root of $\Delta$ and negative for $x \in (\frac{c_1}{b_1}, +\infty)$ and $a_r$ smaller than the greatest root of $\Delta$. Consequently, the two supplementary roots of $P$ may only appear in the interval $(a_1, \frac{c_1}{b_1})$, which implies the apparition of two supplementary endemic equilibria. Moreover, in the degenerate cases where $a_r$ is equal to one of these two roots (case 2b), only a single supplementary endemic equilibrium appears.

* Third, if a unique real zero exists, it is the degenerate case where the two roots of the previous case are equal. In this case, two endemic equilibria exist at this critical point (case 2b), and one otherwise (case 2a).

To conclude, let us show that the roots of $\Delta$ must be positive in the present case. For negative values of $a_r$, due to $Q(x) < 0$ for $x < a_1$, $Q(a_1) = 0$, $\lim_{x \to -\infty} xP(x) < 0$, $a_1 P(a_1) < 0$, then $P(x) = xP(x) + a_r Q(x)$ is negative at $x = a_1$, positive at $x = 0$ and $\lim_{x \to -\infty} P = -\infty$. Thus $P$ has always three real roots when, $a_r < 0$ and $\Delta$ does not vanish in this setting.

All cases are shown to enter one of the categories depicted in the statement. This completes the proof of Theorem 5. $\square$

## B.3 Proof of Theorem 2

Assume that $S_1^*$ is smaller than $S_1^{**}$, let us check first that $I_1^*$ is larger than $I_1^{**}$. Differentiating first $\Psi_r$ in (B.5) give:

$$ \begin{aligned} \Psi'_r(x) &= \frac{x - a_1}{c_1 - b_1 x} + x \frac{c_1 - b_1 x + b_1(x - a_1)}{(c_1 - b_1 x)^2}, \\ &= \frac{(x - a_1)(c_1 - b_1 x) + x(c_1 - b_1 x + b_1(x - a_1))}{(c_1 - b_1 x)^2} > 0, \quad \forall x \in (a_1, \frac{c_1}{b_1}). \end{aligned} \quad (B.10) $$
Thus the endemic value $I_r^*$ increases along with $I_1^*$. We deduce that

$$ \frac{\mu}{\beta_{1,1}x + \beta_{1,r}\Psi_r(x) + \mu} $$

is a decreasing function in $(a_1, \frac{c_1}{b_1})$. Therefore, thanks to (B.2) and $S_1^* \le S_1^{**}$, we conclude that $I_1^{**} \le I_1^*$ and (also) $I_r^{**} \le I_r^*$. The inequalities for $E_1^*, R_1^*, E_r^*, R_r^*$ are easily deduced from (B.2) using $I_r^{**} \le I_r^*$.

It remains to check $S_r^* \ge S_r^{**}$. One has the following equality, derived from (B.2)

$$ S_r^* = \frac{\sigma_r + \mu}{\beta_{r,1}I_1^* + \beta_{r,r}I_r^*} E_r^* = \frac{\sigma_r + \mu}{\beta_{r,1}I_1^* + \beta_{r,r}I_r^*} \frac{\gamma_r + \mu}{\sigma_r} I_r^* = \frac{\sigma_r + \mu}{\sigma_r} \frac{\gamma_r + \mu}{\beta_{r,r}} \frac{I_r^*}{B_rI_1^* + I_r^*} $$

with $B_r$ defined in (5). Let us define then the function $\xi$ as

$$ \xi(x) := \frac{\sigma_r + \mu \gamma_r + \mu}{\sigma_r \beta_{r,r}} \frac{\Psi_r(x)}{B_r x + \Psi_r(x)}. $$

Its derivative is

$$ \begin{aligned} \xi'(x) &= \frac{\sigma_r + \mu \gamma_r + \mu \Psi'_r(x)(B_r x + \Psi_r(x)) - \Psi_r(x)(B_r + \Psi'_r(x))}{\sigma_r \beta_{r,r} (B_r x + \Psi_r(x))^2} \\ &= \frac{\sigma_r + \mu \gamma_r + \mu \Psi'_r(x)B_r x - B_r \Psi_r(x)}{\sigma_r \beta_{r,r} (B_r x + \Psi_r(x))^2} = \frac{\sigma_r + \mu \gamma_r + \mu B_r (\Psi'_r(x)x - \Psi_r(x))}{\sigma_r \beta_{r,r} (B_r x + \Psi_r(x))^2}. \end{aligned} $$

Notice that, thanks to (B.10)

$$ x\Psi'_{r}(x) = \Psi_{r}(x) + x^{2} \frac{c_{1} - a_{1}b_{1}}{(c_{1} - b_{1}x)^{2}} > \Psi_{r}(x) $$

We conclude that $\xi'(x) > 0$, $\forall x \in (a_1, \frac{c_1}{b_1})$, and as $I_1^* \ge I_1^{**}$, then $S_r^* \ge S_r^{**}$. $\square$

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

## B.4 Proof of Theorem 6

The proof requires the following definitions (from [12]) and lemmas.

**Definition B.1.** A nonempty subset $M \subset X$ is called an isolated set if there exists $\varepsilon > 0$ such that for any invariant set $N$ contained entirely in $\{x : x \in X, d(x, M) \le \varepsilon\}$, we have $N \subset M$.

**Definition B.2.** The flow $\phi$ is point dissipative over a nonempty set $M \subset X$ if there exists a compact set $N \subset X$ such that for any $y \in M$, there exists $t(y) > 0$ such that for any $t \ge t(y)$, $\phi_t(y) \in \dot{N}$.

**Lemma B.3.** If the initial condition of system (4) fulfills $E_1(0) + I_1(0) + E_r(0) + I_r(0) > 0$, then $E_1(t) > 0, I_1(t) > 0, E_r(t) > 0, I_r(t) > 0$ for any $t > 0$.

*Proof.* Let us begin by observing that every state $g \in \{S_1, E_1, I_1, R_1, S_r, E_r, I_r, R_r\}$ of the system (4) fulfills the inequality $\dot{g} \ge -a_g g$ where $a_g$ is a non-negative constant. Hence, if any state of the system (4) is positive at a time $T \ge 0$, the positivity is preserved for any $t > T$ thanks to a straightforward application of Grönwall Lemma. Therefore, the Lemma is verified if the initial state vector is positive element-wise, and it remains to check that the Lemma is still valid if the vector of initial condition has some zero components.

If the initial condition fulfills $S_1(0) = 0$, $\dot{S}_1$ is strictly positive and $S_1(t)$ for $t > 0$. Thus, $S_1$ is positive for any $t > 0$ in any case. Adding $\dot{I}_1$ to $\dot{I}_r$, we deduce from (4) that

$$ \dot{I}_1 + \dot{I}_r > \min(\sigma_1, \sigma_r)(E_1 + E_r) - (\max(\gamma_1, \gamma_r) + \mu)(I_1 + I_r). $$

Thus, the condition $E_1(0) + I_1(0) + E_r(0) + I_r(0) > 0$ implies either $I_1(0) + I_r(0) > 0$ or $I_1(0) + I_r(0) = 0$ and $E_1(0) + E_r(0) > 0$. In the second case, thanks to the inequality above, we see that $\dot{I}_1 + \dot{I}_r > 0$ on the right side of $t = 0$. In any case, one deduces $I_1(t) + I_r(t) > 0$ for $t > 0$. This and the fact that $S_1(t) > 0$ for $t > 0$ force $\dot{E}_1$ to be positive on the right side of $t = 0$ in the hypothetical situation where $E_1(0) = 0$. Therefore, $E_1(t) > 0$ for any $t > 0$ also. Finally, the repeated application of the same reasoning leads successively to $I_1(t) > 0, R_1(t) > 0, S_r(t) > 0, E_r(t) > 0, I_r(t) > 0, R_r(t) > 0$ for $t > 0$ which completes the proof. $\square$

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

**Lemma B.4.** Assume that

$$ \lim_{t \to \infty} E_1 = 0, \quad \lim_{t \to \infty} I_1 = 0, \quad \lim_{t \to \infty} E_r = 0, \quad \lim_{t \to \infty} I_r = 0, $$

then,

$$ \lim_{t \to \infty} S_1 = 1, \quad \lim_{t \to \infty} R_r = 0, \quad \lim_{t \to \infty} S_r = 0, \quad \lim_{t \to \infty} R_r = 0. $$
Proof. For $\varepsilon > 0$, there is $T > 0$ such that for $t > T$

$$E_1(t) < \varepsilon, \quad I_1(t) < \varepsilon, \quad E_r(t) < \varepsilon, \quad I_r(t) < \varepsilon.$$

Then, we have the following differential inequalities for $t > T$:

$$\begin{aligned} & \mu - \mu S_1 - (\beta_{1,1} + \beta_{1,r})\varepsilon \le \dot{S}_1 \le \mu - \mu S_1, \\ & -(\omega_1 + \mu)R_1 \le \dot{R}_1 \le \gamma_1\varepsilon - (\omega_1 + \mu)R_1, \\ & -\mu S_r \le \dot{S}_r \le (\gamma_1 + \gamma_r)\varepsilon - \mu S_r, \\ & -(\omega_r + \mu)R_r \le \dot{R}_r \le \gamma_r\varepsilon - (\omega_r + \mu)R_r, \end{aligned}$$

which imply, for $t > T$,

$$\begin{aligned} & S_1(T)e^{-\mu(t-T)} + \int_T^t [\mu - \varepsilon(\beta_{1,1} + \beta_{1,r})]e^{-\mu(t-s)}ds \le S_1(t) \le S_1(T)e^{-\mu(t-T)} + \int_T^t \mu e^{-\mu(t-s)}ds, \\ & E_1(T)e^{-(\omega_1+\mu)(t-T)} \le E_1(t) \le E_1(T)e^{-(\omega_1+\mu)(t-T)} + \int_T^t \gamma_1\varepsilon e^{-(\omega_1+\mu)(t-s)}ds, \\ & S_r(T)e^{-\mu(t-T)} \le S_r(t) \le S_r(T)e^{-\mu(t-T)} + \int_T^t (\gamma_1 + \gamma_r)\varepsilon e^{-\mu(t-s)}ds, \\ & E_r(T)e^{-(\omega_r+\mu)(t-T)} \le E_r(t) \le E_r(T)e^{-(\omega_r+\mu)(t-T)} + \int_T^t \gamma_r\varepsilon e^{-(\omega_r+\mu)(t-s)}ds, \end{aligned}$$

As $\lim_{t \to +\infty} \int_T^t e^{-\lambda(t-s)}ds = \lambda^{-1}$ for $\lambda \in \mathbb{R}$, the inequalities lead to

$$\begin{aligned} & 1 - \mu^{-1}\varepsilon(\beta_{1,1} + \beta_{1,r}) \le \liminf_{t \to \infty} S_1(t) \le \limsup_{t \to \infty} S_1(t) \le 1, \\ & 0 \le \liminf_{t \to \infty} E_1(t) \le \limsup_{t \to \infty} E_1(t) \le (\omega_1 + \mu)^{-1}\gamma_1\varepsilon, \\ & 0 \le \liminf_{t \to \infty} S_r(t) \le \limsup_{t \to \infty} S_r(t) \le \mu^{-1}(\gamma_1 + \gamma_r)\varepsilon, \\ & 0 \le \liminf_{t \to \infty} E_r(t) \le \limsup_{t \to \infty} E_r(t) \le (\omega_r + \mu)^{-1}\gamma_r\varepsilon. \end{aligned}$$

As the previous inequalities are valid for all $\varepsilon > 0$, taking $\varepsilon \to 0$ we obtain

$$\lim_{t \to \infty} S_1(t) = 1, \lim_{t \to \infty} E_1(t) = 0, \lim_{t \to \infty} S_r(t) = 0, \lim_{t \to \infty} E_r(t) = 0.$$

This completes the proof.

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

**Lemma B.5.** For $\mathcal{R}_0 > 1$, no solution of system (4) converges to the DFE except those starting on the invariant set $\partial\Sigma_{\max}$ (17).

Proof. Assume that the contrary is true, i.e., that for $\mathcal{R}_0 > 1$ there is a solution $x$ converging to the DFE and such that $x(0) \notin \partial\Sigma_{\max}$ (note that Lemma B.3 implies then $x(t) \notin \partial\Sigma_{\max}$ for any $t > 0$). Let $L$ be $L := E_1 + \frac{\sigma_1+\mu}{\sigma_1}I_1$, the derivative $\dot{L}$ along trajectories of system (4) is

$$\dot{L} = (\beta_{1,1}I_1 + \beta_{1,r}I_r)S_1 - \frac{(\sigma_1 + \mu)(\gamma_1 + \mu)}{\sigma_1}I_1.$$

Let $\varepsilon_1 > 0$ such that $(1 - \varepsilon_1)\beta_{1,1} > (\gamma_1 + \mu)\frac{\sigma_1+\mu}{\sigma_1}$ (this is possible because $\mathcal{R}_0 > 1$). As the solution $x$ converges toward the DFE, there is a $T > 0$ such that $1 - \varepsilon_1 < S_1(t) < 1$ for $t > T$. Thus for $t > T$, we have the differential inequality:

$$\dot{L} > [(1 - \varepsilon_1)\beta_{1,1} - \frac{(\sigma_1 + \mu)(\gamma_1 + \mu)}{\sigma_1}]I_1 > 0,$$

and the quantity $L(t) = E_1(t) + \frac{\sigma_1+\mu}{\sigma_1}I_1(t)$ is strictly increasing for $t > T$, which contradicts the convergence to the DFE. Finally, the convergence towards the DFE in the invariant set $\partial\Sigma_{\max}$ is self-evident, as the dynamics (4) are reduced to the following system in the set $\partial\Sigma_{\max}$

$$\begin{aligned} & \dot{S}_1 = \mu - \mu S_1, \\ & \dot{R}_1 = -(\omega_1 + \mu)R_1, \\ & \dot{S}_r = -\mu S_r, \\ & \dot{R}_r = -(\omega_r + \mu)R_r, \end{aligned}$$

whose solutions converge to the point $(1, 0, 0, 0)$.

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

The previous lemmas and results from [12] allow us to derive that the system (4) is uniformly persistent.
**Proof of Theorem 6.** The non persistence is deduced from the fact that the DFE is locally asymptotically stable when $\mathcal{R}_0 < 1$.

The uniform persistence follows from an application of [12, Theorem 4.3]. More precisely, let $X$ be the space $\mathbb{R}_+^8$ and $\Sigma$, $\partial\Sigma_{\max}$ defined as in (16) and (17). We show that $\partial\Sigma_{\max}$ is the maximal invariant set in $\partial\Sigma$. Indeed, taking $x_0 \in \partial\Sigma$ but $x_0 \notin \partial\Sigma_{\max}$, then $E_1(0) + I_1(0) + E_r(0) + I_r(0) > 0$ and Lemma B.3 implies that $\phi_t(x_0) \in \bar{\Sigma}$ for $t > 0$, therefore we have proved that $\partial\Sigma_{\max}$ is the maximal invariant set in $\partial\Sigma$.

Moreover, notice that the hypothesis $(\mathcal{H})$ that ensures in the aforementioned theorem by Freedman the nonexistence of cycle in $\partial\Sigma_{\max}$ is fulfilled, as the maximal invariant set $\partial\Sigma_{\max}$ is closed and connected. Indeed, the $(\mathcal{H})$ hypothesis requires that $\partial\Sigma_{\max}$ is a closed invariant set and there exists a cover $\{N_\alpha\}_{\alpha \in A}$ of $\partial\Sigma_{\max}$, where $A$ is a nonempty index set, $N_\alpha \subset \partial E$, $\partial\Sigma_{\max} \subset \cup_{\alpha \in A} N_\alpha$ and $N_\alpha(\alpha \in A)$ are pairwise disjoint closed invariant sets. Furthermore, it is assumed that

* All $N_\alpha$ are isolated invariant sets of the flow $\phi$.

* $\{N_\alpha\}_{\alpha \in A}$ is acyclic, that is, any finite subset of $\{N_\alpha\}_{\alpha \in A}$ does not form a cycle (see [6]).

* Any compact subset of $\partial\Sigma$ contains, at most, finitely many sets of $\{N_\alpha\}_{\alpha \in A}$.

As the set $\partial\Sigma_{\max}$ is already closed and connected, we take simply the cover $\{N_\alpha\}_{\alpha \in A} = \{\partial\Sigma_{\max}\}$ and the three points above can be verified. Furthermore, the system is naturally point-dissipative over the invariant compact set $\Sigma$, thus [12, Theorem 4.3] can be applied. According to this Theorem and Lemma B.3, uniform persistence of the system is equivalent to the non-existence of trajectories starting in $\bar{\Sigma}$ and approaching asymptotically the set $\partial\Sigma_{\max}$. Now, due to Lemma B.4, the convergence towards the set $\partial\Sigma_{\max}$ implies, in fact, the convergence towards the DFE, which is indeed excluded when $\mathcal{R}_0 > 1$ due to Lemma B.5. $\square$

## B.5 Proofs of Theorems 7 and 10

**Proof of Theorem 7.** As $\beta_{1,1} = \beta_{1,r} = \beta$, $\beta_{r,1} = \beta_{r,r} = \psi\beta$, $\gamma_1 = \gamma_r = \gamma$, by formulas in (19), $b_1 = 1$ and

$$1 + b_r b_1 - c_r b_1^2 = 1 + b_r - c_r = 1 + \frac{\gamma}{\mu} - 1 - \frac{\gamma}{\mu} = 0.$$

We conclude the uniqueness of endemic equilibrium for $\mathcal{R}_0 > 1$ by applying point 1 of Theorem 5. This implies the same for system (21).

Let us demonstrate now the statement of the Theorem for system (20) and the case $\mathcal{R}_0 \le 1$. Clearly, from the Theorem 3, a sufficient condition for the non-existence of endemic equilibrium is

$$c = a_1^2 + a_1 c_1 b_r - c_r c_1^2 \ge 0. \qquad (B.12)$$

Thanks to (19), the quantity $c$ can also be expressed in terms of epidemiological parameters

$$c = \left(\frac{\mu\beta - \gamma - \mu}{\beta\gamma + \mu}\right)^2 + \frac{\mu\beta - \gamma - \mu}{\beta\gamma + \mu}\frac{\mu}{\gamma + \mu}\left(\frac{\gamma}{\mu} - 1\right) - \frac{\gamma}{\gamma + \mu}\frac{\mu}{\gamma + \mu} \\ = \frac{\mu\beta - \gamma - \mu}{\beta(\gamma + \mu)^2}\left(\mu\frac{\beta - \gamma - \mu}{\beta} + \gamma - \mu - \gamma\frac{\beta}{\beta - \gamma - \mu}\right).$$

As $\mathcal{R}_0 < 1$, thus $\beta - \gamma - \mu < 0$ and the condition (B.12) is equivalent to

$$\mu \frac{\beta - \gamma - \mu}{\beta} + \gamma \le \gamma \frac{\beta}{\beta - \gamma - \mu} + \mu, \\ \frac{(\beta - \gamma - \mu)\mu + \gamma\beta}{\beta} \le \frac{\gamma\beta + \mu(\beta - \gamma - \mu)}{\beta - \gamma - \mu}, \\ (\gamma + \mu) \frac{\beta - \mu}{\beta} \le (\gamma + \mu) \frac{\beta - \mu}{\beta - \gamma - \mu}.$$

Now, if $\beta > \mu$, the last inequality above is the same as

$$\frac{\beta - \gamma - \mu}{\beta} \ge 1,$$

which is always false. On the other hand, if $\beta < \mu$ then the same condition is now equivalent to

$$\frac{\beta - \gamma - \mu}{\beta} \le 1,$$

which is always verified. Finally, if $\beta = \mu$, the same condition is trivially verified. To sum up, the condition (B.12) is the same as

$$\beta \le \mu.$$
Therefore, if $\beta \le \mu$ the system does not admit endemic equilibrium.

For the contrary case, Theorem 3 states that the smallest zero of $\Delta$ acts as a threshold on $a_r$ through which the system passes from a two endemic equilibrium regime to a zero endemic equilibrium regime. As $a_r = \frac{\gamma+\mu}{\psi\beta}$ is inversely proportional to $\psi$, this is equivalent to say that there exists a threshold $\psi^*$ on $\psi$ for the appearance of two endemic equilibria.

Finally, let $\psi = 1$ and denote $I := I_1 + I_r$, $R := S_1 + S_r$, then the system (20) can be reduced to the SIS system:

$$ \dot{S} = \mu - \beta IS - \mu S + \gamma I, $$

$$ \dot{I} = \beta IS - (\gamma + \mu)I. $$

It is easy to verify that this SIS system does not admit endemic equilibrium for $\mathcal{R}_0 < 1$. In particular, this forces $\psi^*$ to be greater than 1.

We have shown that the system (20) exhibits one endemic and one disease-free equilibrium for $\mathcal{R}_0 > 1$, and up to 2 endemic equilibria and one disease-free equilibrium for $\mathcal{R}_0 < 1$. In addition, the endemic equilibria of the two-stage SIS system are endowed with a total ordering in the same way as in Theorem 2 for the SEIRS system. Hence, Theorem 7 still holds for the reduced SIRI model (21), as the multiple endemic equilibria that may arise in system (20) cannot be confounded in the reduced model (21) due to the ordering. $\square$

**Proof of Theorem 10.** The computation is mostly the same as in the proof of Theorem 7 for the case $\mathcal{R}_0 > 1$.

On the other hand, for $\mathcal{R}_0 \le 1$, we prove first the statement for the system (22). The quantity $c$ in this setting is equal to

$$ \begin{aligned} c &= a_1^2 + a_1 c_1 b_r - c_r c_1^2 \\ &= \sigma \frac{\beta \sigma \mu - \mu(\gamma + \mu)(\sigma + \mu)}{\beta(\gamma + \mu)^2(\sigma + \mu)^2(\sigma + \gamma + \mu)} \left( \mu(\sigma + \gamma + \mu) \frac{\beta \sigma - (\gamma + \mu)(\sigma + \mu)}{\beta \sigma} + \gamma \sigma - (\sigma + \gamma + \mu)\mu \right. \\ & \qquad \left. - \gamma \sigma \frac{\beta \sigma}{\beta \sigma - (\gamma + \mu)(\sigma + \mu)} \right). \end{aligned} $$

For $\mathcal{R}_0 < 1$, $\beta\sigma < (\gamma + \mu)(\sigma + \mu)$ and the condition $c \ge 0$ is equivalent to

$$ \frac{\mu(\sigma + \gamma + \mu)(\beta\sigma - (\gamma + \mu)(\sigma + \mu)) + \gamma\sigma^2\beta}{\beta\sigma} \le \frac{\mu(\sigma + \gamma + \mu)(\beta\sigma - (\gamma + \mu)(\sigma + \mu)) + \gamma\sigma^2\beta}{\beta\sigma - (\gamma + \mu)(\sigma + \mu)}. $$

If $\mu(\sigma + \gamma + \mu)(\beta\sigma - (\gamma + \mu)(\sigma + \mu)) + \gamma\sigma^2\beta > 0$, the condition $c \ge 0$ is the equivalent to

$$ \frac{\beta\sigma - (\gamma + \mu)(\sigma + \mu)}{\beta\sigma} \ge 1, $$

which is always false. On the contrary, if $\mu(\sigma + \gamma + \mu)(\beta\sigma - (\gamma + \mu)(\sigma + \mu)) + \gamma\sigma^2\beta \le 0$, the same condition is equivalent to

$$ \frac{\beta\sigma - (\gamma + \mu)(\sigma + \mu)}{\beta\sigma} \le 1, $$

which is verified. To sum up, the condition $c \ge 0$ is equivalent to

$$ \mu(\sigma + \gamma + \mu)(\beta\sigma - (\gamma + \mu)(\sigma + \mu)) + \gamma\sigma^2\beta \le 0, $$

that is (24). Similarly to the proof of Theorem 7, the existence of the threshold $\psi^*$ is then deduced using Theorem 3. Moreover, for $\psi = 1$, we denote $S := S_1 + S_r$, $I := I_1 + I_r$, $E := E_1 + E_r$, then the system (22) can be again reduced to the SEIS system below.

$$ \dot{S} = \mu - \beta SI - \mu S + \gamma I, $$

$$ \dot{E} = \beta SI - (\sigma + \mu)E, $$

$$ \dot{I} = \sigma E - (\gamma + \mu)I. $$

It is known that this system does not admit endemic equilibrium for $\mathcal{R}_0 < 1$ [10]. In particular, this implies that $\psi^* > 1$.

Notice finally that the result still holds for the reduced system (23), thanks to an adaptation of Theorem 2 which gives an ordering of the equilibria. $\square$
## B.6 Proofs of Theorems 9 and 12

**Proof of Theorem 9.** For $\mathcal{R}_0 \neq 1$, let us show that any trajectory of (21) converges to an equilibrium using the Li and Muldowney theory summarized in Appendix A. The invariant and biologically feasible set $\Gamma$ is defined in this case by $\Gamma = \{x \in \mathbb{R}_+^3 | g(x) = 0\}$, where $g: \mathbb{R}^3 \to \mathbb{R}^m$, $g(x) = S + I + R - 1$. In the rest of the proof, we may consider only the trajectories evolving in the interior $\bar{\Gamma}$ of $\Gamma$. This is justified as trajectories starting in $\Gamma$ either converge to the disease-free equilibrium if $I(0) = 0$, or fulfill $S(t), I(t), R(t) > 0$ for $t > 0$ thanks to the differential inequalities:

$$ \dot{S} \geq \mu(1-S), \quad \dot{I} \geq -(\gamma + \mu)I, \quad \dot{R} \geq -\mu R, $$

and the fact that $\dot{S}(0) > 0$ (resp. $\dot{R}(0) > 0$) if $S(0) = 0$ (resp. $R(0) = 0$). In addition, notice that

$$ \frac{\partial g}{\partial x}(x)f(x) = \mu(1 - S - I - R) = -\mu g(x), $$

thus the function $\nu$ introduced in (A.8) of the Appendix A above corresponds to the constant $-\mu$ and $m = 1$ for the system (21). The associated 3-rd additive compound matrix of $\frac{df}{dx}(x)$ associated to the equation is (see [36])

$$ \dot{z} = \left( \frac{df}{dx} \right)^{[3]} (x) + \mu I z. $$

As

$$ \frac{\partial f}{\partial x}(x) = \begin{pmatrix} -\beta I - \mu & -\beta S & 0 \\ \beta I & \beta(S + \psi R) - \gamma - \mu & \psi \beta I \\ 0 & \gamma - \psi \beta R & -\psi \beta I - \mu \end{pmatrix}, \qquad (B.13) $$

is a 3 × 3 matrix, its 3-rd additive compound is simply its trace $\frac{df}{dx}^{[3]} = \text{Tr}(\frac{df}{dx})$. Therefore, the compound equation (A.7) is reduced to:

$$ \dot{z} = [\beta(S + \psi R) - \gamma - \beta(1 + \psi)I - 2\mu]z. $$

Moreover,

$$ \beta(S + \psi R) - \gamma - \beta I(1 + \psi) - 2\mu = \frac{\dot{I}}{I} - \beta I(1 + \psi) - \mu. $$

Then, after integration of the compound equation, we obtain

$$ z(t) = z(0) \frac{I(t)}{I(0)} e^{-\int_{0}^{t} (\beta(1+\psi)I(s)+\mu)ds}. $$

As $\lim_{t \to +\infty} \int_0^t (\beta(1 + \psi)I(s) + \mu)ds = \infty$, thus $\lim_{t \to +\infty} z(t) = 0$ exponentially, and the system (21) satisfies a Bendixson condition (see [27] or Appendix A) in $\bar{\Gamma}$. Therefore, Theorem A.2 holds with $D = \bar{\Gamma}$, $Q = I$ and any solution $x(t)$ of (21) either converges to an equilibrium or approaches the boundary of $\bar{\Gamma}$. The later case is precluded for $\mathcal{R}_0 > 1$ due to uniform persistence (Theorem 8). On the other hand, for $\mathcal{R}_0 < 1$, it can be observed that the solutions may only approach the component $\{x \in \mathbb{R}_+^3 | S + R = 1, I = 0\}$ of the boundary. Indeed, we have $\dot{S} > 0$ if $S$ is close to 0 and the same for $\dot{R}$ if $I > 0$ and $R$ close to 0. In addition, from the expression of the Jacobian (B.13), the DFE is locally asymptotically stable, thus its basin of attraction is an open subset of $\Gamma$ [7]. As the DFE attracts $\{x \in \mathbb{R}_+^3 | S + R = 1, I = 0\}$, hence it also attracts a close neighborhood of the latter and the solution converges to the DFE in this case. $\square$

For the proof of Theorem 12, the Lemma below is first required.

**Lemma B.6.** Let $x(t) = (S(t), E(t), I(t), R(t))$ be any solution of (23) in the invariant set $\{x \in \mathbb{R}_+^4 | S+E+I+R = 1\}$. Then ultimately the following inequalities are verified

$$ S \geq \frac{\mu}{\beta + \mu}, \quad \gamma \geq \psi \beta R. $$

*Proof of Lemma B.6.* The first inequality is deduced from the inequality

$$ \dot{S} \geq \mu - (\beta + \mu)S. $$

Notice that the derivative of $R$, solution of (23), is non-positive on the set $\{x \in \mathbb{R}_+^4 | \gamma = \psi \beta R\}$, hence $\{x \in \mathbb{R}_+^4 | \gamma \geq \psi \beta R\}$ is an invariant set of (23).

For trajectories such that $I(0) + E(0) = 0$, it is easy to see that $I(t) = 0$ for $t \geq 0$ and $R(t)$ converges also to 0 in this case, hence $\gamma \geq \psi \beta R$ ultimately.
On the other hand, assume $I(0) + E(0) > 0$, then $I(t), R(t) > 0$ for $t > 0$. This is because, if $I(0) = 0$, $E(0) > 0$ and $\dot{I} > 0$ and $I$ is positive in a positive neighborhood of zero. Thus, we may assume with no loss of generality that $I(0) > 0$. In this case, the differential inequalities

$$ \dot{I}(t) \geq -(\gamma + \mu)I, \quad \dot{R} \geq -(\psi\beta I + \mu)R, $$

and the fact that $\dot{R}(0) > 0$ if $R(0) = 0$ imply that $I(t), R(t) > 0$ for $t > 0$ anyway. Therefore, the fact that $I(t), R(t) > 0$ for all $t > 0$ implies $\dot{R}(t) < 0$ if $\gamma \le \psi\beta R(t)$ and again $\gamma \ge \psi\beta R$ ultimately. $\square$

**Proof of Theorem 12.** In a similar fashion as in the proof of Theorem 9 for the SIR system above, notice first that only trajectories evolving in the interior of $\Gamma := \{x \in \mathbb{R}_+^4 | S+E+I+R = 1\}$ matter, as any solution in $\Gamma$ either fulfills $S(t), E(t), I(t), R(t) > 0$ for $t > 0$ or converges to the disease-free equilibrium when $E(0) + I(0) = 0$. We apply again the geometrical approach of Li and Muldowney [28]. The Jacobian matrix of (23) is

$$ \frac{df}{dx} = \begin{pmatrix} -(\beta I + \mu) & 0 & -\beta S & 0 \\ \beta I & -(\sigma + \mu) & \beta S + \psi\beta R & \psi\beta I \\ 0 & \sigma & -(\gamma + \mu) & 0 \\ 0 & 0 & -\psi\beta R + \gamma & -(\psi\beta I + \mu) \end{pmatrix}. $$

and its 3-rd additive compound is defined by (see the Appendix of [36])

$$ \frac{df}{dx} = \begin{pmatrix} k_1 & 0 & -\psi\beta I & 0 \\ -\psi\beta R + \gamma & k_2 & \beta S + \psi\beta R & \beta S \\ 0 & \sigma & k_3 & 0 \\ 0 & 0 & \beta I & k_4 \end{pmatrix}, $$

where

$$ \begin{aligned} k_1 &= -\beta I - \sigma - \gamma - 3\mu, & k_2 &= -(1 + \psi)\beta I - \sigma - 3\mu, \\ k_3 &= -(1 + \psi)\beta I - \gamma - 3\mu, & k_4 &= -\sigma - \gamma - \psi\beta I - 3\mu. \end{aligned} $$

Hence, the associated compound equation (A.7) is in this case defined as

$$ \dot{z} = \left( \frac{df}{dx}^{[3]} + \mu \mathbf{I} \right) z. $$

Let us introduce now

$$ A := \operatorname{diag}(R, I, E, S), $$

then

$$ A^{-1} = \operatorname{diag}(R^{-1}, I^{-1}, E^{-1}, S^{-1}) \quad \text{and} \quad A_f A^{-1} = \operatorname{diag}(\frac{\dot{R}}{R}, \frac{\dot{I}}{I}, \frac{\dot{E}}{E}, \frac{\dot{S}}{S}), $$

where $A_f$ stands for the directional derivative of the matrix $A$ in the direction of the vector field $f$. One has in addition,

$$ A \frac{df}{dx}^{[3]} A^{-1} = \begin{pmatrix} k_1 & 0 & \psi \beta I \frac{R}{E} & 0 \\ (\gamma - \psi \beta R) \frac{I}{R} & k_2 & (\beta S + \psi \beta R) \frac{I}{E} & \beta I \\ 0 & \sigma \frac{E}{I} & k_3 & 0 \\ 0 & 0 & \beta I \frac{S}{E} & k_4 \end{pmatrix}. $$

Thus,

$$ B := A_f A^{-1} + A \frac{df}{dx}^{[3]} A^{-1} + \mu I = \begin{pmatrix} k_5 & 0 & \psi \beta I \frac{R}{E} & 0 \\ (\gamma - \psi \beta R) \frac{I}{R} & k_6 & (\beta S + \psi \beta R) \frac{I}{E} & \beta I \\ 0 & \sigma \frac{E}{I} & k_7 & 0 \\ 0 & 0 & \beta I \frac{S}{E} & k_8 \end{pmatrix}, $$

where

$$ \begin{aligned} k_5 &= \frac{\dot{R}}{R} - \beta I - \sigma - \gamma - 2\mu, & k_6 &= \frac{\dot{I}}{I} - (1 + \psi)\beta I - \sigma - 2\mu, \\ k_7 &= \frac{\dot{E}}{E} - (1 + \psi)\beta I - \gamma - 2\mu, & k_8 &= \frac{\dot{S}}{S} - \sigma - \gamma - \psi\beta I - 2\mu. \end{aligned} $$

Thanks to Lemma B.6, we may assume directly that $\gamma - \psi\beta R \ge 0$ in the sequel. Defining $y(t) := A(x(t))z(t)$, $y$ is determined by the differential equation

$$ \dot{y} = By. $$
Let us prove (A.4) with $Q(x) = A(x)$ and considering the $L^\infty$ norm on the vector space $\mathbb{R}^4$. For a given matrix $M = (M_{i,j})_{1\le i,j\le 4}$ and denoting

$$l_i := M_{i,i} + \sum_{j \ne i} |M_{i,j}|,$$

the Lozinskii measure $\mu_\infty$ associated to the $L^\infty$ norm is defined by (see for example [35])

$$\mu_\infty(M) = \max_{1 \le i \le 4} \{l_i\}.$$

We provide the following formulas, which are required for the computation of $l_i$.

$$\begin{align}
\frac{\dot{S}}{S} &= \mu S^{-1} - \beta I - \mu, \quad \frac{\dot{E}}{E} = \beta S \frac{I}{E} + \psi \beta \frac{I}{E} R - (\sigma + \mu), \tag{B.14} \\
\frac{\dot{I}}{I} &= \sigma \frac{E}{I} - (\gamma + \mu), \quad \frac{\dot{R}}{R} = \gamma \frac{I}{R} - \psi \beta I - \mu. \nonumber
\end{align}$$

Let us compute now the value of $l_i$ for the matrix $B$.

$$\begin{align*}
l_1 &= \frac{\dot{R}}{R} - \beta I - \sigma - \gamma - 2\mu + \psi \beta I \frac{R}{E} \\
&= \frac{\dot{E}}{E} + \frac{\dot{R}}{R} - \beta I \left(1 + \frac{I}{E}\right) - \gamma - \mu, \\
l_2 &= \frac{\dot{I}}{I} - (1 + \psi)\beta I - \sigma - 2\mu + (\gamma - \psi \beta R) \frac{I}{R} + (\beta S + \psi \beta R) \frac{I}{E} + \beta I \\
&= \frac{\dot{E}}{E} + \frac{\dot{I}}{I} + \frac{\dot{R}}{R} - \psi \beta I, \\
l_3 &= \frac{\dot{E}}{E} - (1 + \psi)\beta I - \gamma - 2\mu + \sigma \frac{E}{I} \\
&= \frac{\dot{E}}{E} - (1 + \psi)\beta I - \mu, \\
l_4 &= \frac{\dot{S}}{S} - \sigma - \gamma - \psi \beta I - 2\mu + \beta I \frac{S}{E} \\
&= \frac{\dot{S}}{S} + \frac{\dot{E}}{E} - \gamma - \psi \beta I \left(1 + \frac{R}{E}\right) - \mu.
\end{align*}$$

Therefore, an easy adaptation of [31, Theorem 2.2] to the norm $L^\infty$ makes it possible to say that it is sufficient, for all $1 \le i \le 4$, to have $\bar{g}_i$ such that, for any trajectory of (23) in a given invariant set $\Sigma$,

$$\limsup_{t \to +\infty} \frac{1}{t} \int_0^t l_i(s) ds \le \bar{g}_i < 0, \qquad (B.15)$$

in order to guarantee (A.9) and for the system (23) to satisfy a Bendixson condition in this invariant set. For $l_1$, $l_3$, $l_4$, the property (B.15) is verified for any trajectory in the interior of $\Gamma$. For example, one has

$$\limsup_{t \to +\infty} \frac{1}{t} \int_0^t l_1(s) ds \le \lim_{t \to +\infty} \frac{1}{t} \int_0^t \left( \frac{\dot{E}}{E} + \frac{\dot{R}}{R} - \gamma - \mu \right) ds = \lim_{t \to +\infty} \frac{1}{t} \left( \log\left(\frac{E(t)}{E(0)}\right) + \log\left(\frac{R(t)}{R(0)}\right) \right) - \gamma - \mu = -(\gamma + \mu).$$

However, for $l_2$, the property is not always verified, as we shall see.

$$\begin{align*}
\frac{1}{t} \int_{0}^{t} l_{2}(s) ds &= \frac{1}{t} \left[ \log\left(\frac{E(t)}{E(0)}\right) + \log\left(\frac{I(t)}{I(0)}\right) + \log\left(\frac{R(t)}{R(0)}\right) \right] - \psi \beta \frac{1}{t} \int_{0}^{t} I(s) ds, \\
\limsup_{t \to +\infty} \frac{1}{t} \int_{0}^{t} l_{2}(s) ds &= -\psi \beta \liminf_{t \to +\infty} \frac{1}{t} \int_{0}^{t} I(s) ds.
\end{align*}$$

For $\mathcal{R}_0 > 1$, the persistence result in Theorem 11 guarantees the existence of $\varepsilon > 0$ such that

$$\liminf_{t \to +\infty} \frac{1}{t} \int_{0}^{t} I(s) ds > \varepsilon,$$

and the property (B.15) is indeed verified for any trajectory in the interior $\dot{\Gamma}$ of $\Gamma$. The uniform persistence gives also a compact absorbing subset (see e.g. [27] for a definition) $\mathcal{K}$ in $\dot{\Gamma}$, where there exists $\varepsilon_\mathcal{K} > 0$ such that $S, E, I, R \ge \varepsilon_\mathcal{K}$ and $\|A^{-1}\|$ is uniformly bounded in $\mathcal{K}$. Therefore, the result in [31, Theorem 2.4] adapted to the
$L^\infty$ norm allows us to deduce that any trajectory starting in the interior of $\Gamma$ converges to the unique endemic equilibrium when $\mathcal{R}_0 > 1$.

On the other hand, for $\mathcal{R}_0 < 1$, no persistence result is available, and we need to first examine under which circumstances the solution of (23) may fulfill

$$ \liminf_{t \to +\infty} \frac{1}{t} \int_0^t I(s)ds = 0. \qquad (B.16) $$

An integration of the equation for $\dot{I}$ in (23) gives

$$ 0 = \limsup_{t \to +\infty} \frac{1}{t} [I(t) - I(0)] = \sigma \limsup_{t \to +\infty} \frac{1}{t} \int_0^t E(s)ds - (\gamma + \mu) \liminf_{t \to +\infty} \frac{1}{t} \int_0^t I(s)ds. $$

Assume $\liminf_{t \to +\infty} \frac{1}{t} \int_0^t I(s)ds = 0$, the equation above implies $\limsup_{t \to +\infty} \frac{1}{t} \int_0^t E(s)ds = 0$. As $E(t) \ge 0$, hence $\lim_{t \to +\infty} \frac{1}{t} \int_0^t E(s)ds$ is then defined and equal to zero. Integrating again the same differential equation for $\dot{I}$, we obtain $\lim_{t \to +\infty} \frac{1}{t} \int_0^t I(s)ds = 0$. Furthermore, the differential inequality

$$ \dot{R} \le \gamma I - \mu R, $$

gives

$$ 0 \le \gamma \lim_{t \to +\infty} \frac{1}{t} \int_0^t I(s)ds - \mu \lim_{t \to +\infty} \frac{1}{t} \int_0^t R(s)ds. $$

As $\lim_{t \to +\infty} \frac{1}{t} \int_0^t I(s)ds = 0$ and $R \ge 0$, thus $\lim_{t \to +\infty} \frac{1}{t} \int_0^t R(s)ds = 0$. Finally, thanks to $S + E + I + R = 1$, $\lim_{t \to +\infty} \frac{1}{t} \int_0^t S(s)ds = 1$. In this scenario, the solution must visit any neighborhood containing points arbitrarily close to the disease-free equilibrium. Otherwise, there would be $\varepsilon > 0$ such that $S(t) < 1 - \varepsilon$ for $t \ge 0$, which contradicts $\lim_{t \to +\infty} \frac{1}{t} \int_0^t S(s)ds = 1$. Hence, for all $\varepsilon > 0$, there must be a $t \ge 0$ such that $S(t) \ge 1 - \varepsilon$. As $E(t), I(t), R(t) \ge 0$ and $S(t) + E(t) + I(t) + R(t) = 1$, then $I(t), E(t), R(t) \le \varepsilon$. As $\varepsilon$ can be made arbitrarily small, indeed this is equivalent to saying that the solution visits any neighborhood of the DFE. Moreover, the Jacobian matrix at the disease-free equilibrium

$$ \begin{pmatrix} -\mu & 0 & -\beta & 0 \\ 0 & -(\sigma + \mu) & \beta & 0 \\ 0 & \sigma & -(\gamma + \mu) & 0 \\ 0 & 0 & \gamma & -\mu \end{pmatrix} $$

is Hurwitz for $\mathcal{R}_0 < 1$. Indeed, $-\mu$ is an eigenvalue of the matrix of algebraic multiplicity 2 and the remaining eigenvalues are the same as those of the matrix

$$ \begin{pmatrix} -(\sigma + \mu) & \beta \\ \sigma & -(\gamma + \mu) \end{pmatrix}. $$

As the $2 \times 2$ matrix has a negative trace and moreover its determinant $(\sigma + \mu)(\sigma + \mu) - \beta\sigma$ is positive when $\mathcal{R}_0 < 1$, this yields the local stability of the DFE. Therefore, if $\liminf_{t \to +\infty} \frac{1}{t} \int_0^t I(s)ds = 0$, then the trajectory must converge to the disease-free equilibrium, as it is locally asymptotically stable for $\mathcal{R}_0 < 1$. Inversely, if the solution converges to the DFE, naturally (B.16) is verified. Hence, the solutions verifying (B.16) are those starting in $\Gamma$ and in the basin of attraction of the DFE, which we refer to as $\Lambda \subset \mathbb{R}^4$.

For $x_0 \in \dot{\Gamma}$, we denote now $x(t, x_0) = (S(t, x_0), E(t, x_0), I(t, x_0), R(t, x_0))$ the solution of (23) with initial condition $x_0$. Notice then that

$$ \liminf_{t \to +\infty} \frac{1}{t} \int_0^t I(s, x)ds = \sup_{t \ge 0} \inf_{s \ge t} \frac{1}{s} \int_0^s I(\tau, x)d\tau, $$

is lower semi-continuous with regard to the variable $x$. Moreover, $\Gamma \setminus \Lambda$ is closed and compact ($\Lambda$ is open as the basin of attraction of the locally asymptotically stable DFE [7]). Therefore, there exists $\varepsilon > 0$ such that

$$ \inf_{x_0 \in \Gamma \setminus \Lambda} \liminf_{t \to +\infty} \frac{1}{t} \int_0^t I(s, x_0)ds = \min_{x_0 \in \Gamma \setminus \Lambda} \liminf_{t \to +\infty} \frac{1}{t} \int_0^t I(s, x_0)ds > \varepsilon. $$

Otherwise there exists $x_0 \in \Gamma \setminus \Lambda$ such that $\liminf_{t \to +\infty} \frac{1}{t} \int_0^t I(s, x_0)ds = 0$ and the solution $x(t, x_0)$ converges to the DFE, which contradicts the definition of $\Gamma \setminus \Lambda$.

Hence, the property (B.15) is satisfied in $\Gamma \setminus \Lambda$. However, it is not guaranteed that $\|A^{-1}\|$ is bounded in $\Gamma \setminus \Lambda$, and we have to show that there is an absorbing subset $\Sigma$ of $\Gamma \setminus \Lambda$ where the inverse is defined and bounded.
Thanks to Lemma B.6, an ultimate lower bound

$$ \liminf_{t \to +\infty} S(t) \ge \varepsilon_S := \frac{\mu}{\beta + \mu} > 0 $$

is already derived for any trajectory in $\Gamma$. It remains to derive similar bounds for $E, I$ and $R$. Notice first that $\partial\Gamma_1 := \{x \in \Gamma | E + I = 0\}$ is a compact subset of $\Lambda$. As $\Lambda$ is open, for each point of $x_\Lambda \in \partial\Gamma_1$, there exists $\varepsilon_{x_\Lambda}$ such that $C_{x_\Lambda} := \{x \in \mathbb{R}^4 \|x - x_\Lambda\| < \varepsilon_{x_\Lambda}\} \subset \Lambda$. Thus, $\{C_x\}_{x \in \partial\Gamma_1}$ constitutes a cover of $\partial\Gamma_1$ by open subsets of $\Lambda$. From the compactness of $\partial\Gamma_1$, we extract a finite subcover from $\{C_x\}_{x \in \partial\Gamma_1}$, and this allows to deduce the existence of $\varepsilon_{E+I} > 0$ such that

$$ \min_{(S,E,I,R) \in \Gamma \setminus \Lambda} E + I \ge \varepsilon_{E+I}. $$

This lower bound and equations (23) imply the differential inequality in the invariant set $\Gamma \setminus \Lambda$:

$$ \dot{I} \ge \sigma \varepsilon_{E+I} - (\gamma + \sigma + \mu)I. $$

Hence, any trajectory in $\Gamma \setminus \Lambda$ fulfills ultimately

$$ \liminf_{t \to +\infty} I(t) \ge \varepsilon_I := \frac{\sigma \varepsilon_{E+I}}{\gamma + \sigma + \mu} > 0. $$

Again from (23), this implies that ultimately, the differential inequalities below are verified for any solution in $\Gamma \setminus \Lambda$.

$$ \liminf_{t \to +\infty} \dot{E}(t) \ge \beta \varepsilon_S \varepsilon_I - (\sigma + \mu)E, \quad \liminf_{t \to +\infty} \dot{R}(t) \ge \gamma \varepsilon_I - \mu R. $$

Thus, there exists also $\varepsilon_E, \varepsilon_R > 0$ such that, for any trajectory in $\Gamma \setminus \Lambda$, the inequalities

$$ \liminf_{t \to +\infty} E(t) \ge \varepsilon_E, \quad \liminf_{t \to +\infty} R(t) \ge \varepsilon_R, $$

are verified. The compact $\Gamma \setminus \Lambda \cap \{(S, E, I, R) \in \Gamma | S \ge \varepsilon_S, E \ge \varepsilon_E, I \ge \varepsilon_I, R \ge \varepsilon_R\}$ is then an absorbing subset of $\Gamma \setminus \Lambda$ where $\|A^{-1}\|$ is uniformly bounded and the Bendixson criterion (B.15) is satisfied. An application of [28, Theorem 6.1] concludes that any trajectory in the invariant set $\Gamma \setminus \Lambda$ converges to an equilibrium. This completes the proof of Theorem 12. $\square$