# An Age of Infection Kernel, an R Formula, and Further Results for Arino-Brauer A, B Matrix Epidemic Models with Varying Populations, Waning Immunity, and Disease and Vaccination Fatalities

Florin Avram, Rim Adenane, Lasko Basnarkov, Gianluca Bianchin, Dan Goreac, Andrei Halanay

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

## To cite this version:

Florin Avram, Rim Adenane, Lasko Basnarkov, Gianluca Bianchin, Dan Goreac, et al.. An Age of Infection Kernel, an R Formula, and Further Results for Arino-Brauer A, B Matrix Epidemic Models with Varying Populations, Waning Immunity, and Disease and Vaccination Fatalities. Mathematics, 2023, 11 (6), pp.1307. 10.3390/math11061307. hal-04263236

HAL Id: hal-04263236

https://univ-eiffel.hal.science/hal-04263236v1

Submitted on 6 Sep 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Article

# An Age of Infection Kernel, an *R* Formula, and Further Results for Arino-Brauer *A*, *B* Matrix Epidemic Models with Varying Populations, Waning Immunity, and Disease and Vaccination Fatalities

Florin Avram <fnref n="1" />,*,† ID, Rim Adenane <fnref n="2" />,†, Lasko Basnarkov <fnref n="3" />,† ID, Gianluca Bianchin <fnref n="4" />,†, Dan Goreac <fnref n="5" />,<fnref n="6" />,<fnref n="7" />,<fnref n="8" />,† ID
and Andrei Halanay <fnref n="9" />,†

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

check for
updates

**Citation:** Avram, F.; Adenane, R.; Basnarkov, L.; Bianchin, G.; Goreac, D.; Halanay, A. An Age of Infection Kernel, an *R* Formula, and Further Results for Arino-Brauer *A*, *B* Matrix Epidemic Models with Varying Populations, Waning Immunity, and Disease and Vaccination Fatalities.
*Mathematics* **2023**, *11*, 1307. https://doi.org/10.3390/math11061307

Academic Editor: Mikhail Kolev

Received: 12 December 2022

Revised: 18 January 2023

Accepted: 22 January 2023

Published: 8 March 2023

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

**Copyright:** © 2023 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

**Abstract:** In this work, we first introduce a class of deterministic epidemic models with varying populations inspired by Arino et al. (2007), the parameterization of two matrices, demography, the waning of immunity, and vaccination parameters. Similar models have been focused on by Julien Arino, Fred Brauer, Odo Diekmann, and their coauthors, but mostly in the case of “closed populations” (models with varying populations have been studied in the past only in particular cases, due to the difficulty of this endeavor). Our Arino-Brauer models contain SIR-PH models of Riano (2020), which are characterized by the phase-type distribution $(\vec{\alpha}, A)$, modeling transitions in “disease/infectious compartments”. The $A$ matrix is simply the Metzler/sub-generator matrix intervening in the linear system obtained by making all new infectious terms 0. The simplest way to define the probability row vector $\vec{\alpha}$ is to restrict it to the case where there is only one susceptible class $s$, and when matrix $B$ (given by the part of the new infection matrix, with respect to $s$) is of rank one, with $B = b\vec{\alpha}$. For this case, the first result we obtained was an explicit formula (12) for the replacement number (not surprisingly, accounting for varying demography, waning immunity and vaccinations led to several nontrivial modifications of the Arino et al. (2007) formula). The analysis of $(A, B)$ Arino-Brauer models is very challenging. As obtaining further general results seems very hard, we propose studying them at three levels: (A) the exact model, where only a few results are available—see Proposition 2; and (B) a “first approximation” (FA) of our model, which is related to the usually closed population model often studied in the literature. Notably, for this approximation, an associated renewal function is obtained in (7); this is related to the previous works of Breda, Diekmann, Graaf, Pugliese, Vermiglio, Champredon, Dushoff, and Earn. (C) Finally, we propose studying a second heuristic “intermediate approximation” (IA). Perhaps our main contribution is to draw attention to the importance of $(A, B)$ Arino-Brauer models and that the FA approximation is not the only way to tackle them. As for the practical importance of our results, this is evident, once we observe that the $(A, B)$ Arino-Brauer models include a large number of epidemic models (COVID, ILI, influenza, illnesses, etc.).

**Keywords:** epidemic models; varying population models; stability; next-generation matrix approach; basic replacement number; vaccination; waning immunity; endemic equilibria

<footnote n="1">¹ Laboratoire de Mathématiques Appliquées, Université de Pau, 64012 Pau, France</footnote>

<footnote n="2">² Département des Mathématiques, Université Ibn-Tofail, Kenitra 14000, Morocco</footnote>

<footnote n="3">³ Faculty of Computer Science and Engineering, Ss. Cyril and Methodius University in Skopje, 1000 Skopje, North Macedonia</footnote>

<footnote n="4">⁴ ICTEAM & Department Mathematical Engineering, University of Louvain, 1348 Louvain-la-Neuve, Belgium</footnote>

<footnote n="5">⁵ School of Mathematics and Statistics, Shandong University, Weihai 264209, China</footnote>

<footnote n="6">⁶ Laboratoire d'Analyse et de Mathématiques Appliquées, Université Gustave Eiffel, 94010 Creteil, France</footnote>

<footnote n="7">⁷ Maître de Conférences HDR, UPEM, 77420 Champs-sur-Marne, France</footnote>

<footnote n="8">⁸ LAMA UMR8050, Université Paris Est Creteil, 94010 Creteil, France</footnote>

<footnote n="9">⁹ Department of Mathematics and Informatics, Polytechnic University of Bucharest, 060042 Bucharest, Romania</footnote>

<footnote>* Correspondence: avramf3@gmail.com</footnote>

<footnote>† These authors contributed equally to this work.</footnote>
MSC: 34A34; 34A30; 37N25; 92D30; 92D25

# 1. Introduction

The initial motivation of this paper was to extend an outstanding formula derived in
in reference ([1], [Thm 2.1]) to epidemic models with varying population sizes (linear birth
rates), waning immunity, and vaccinations, expressing the replacement numbers of large
classes of deterministic epidemic models as quadratic forms with respect to certain input
matrices—see (15) below. Unfortunately, this formula is not well known, and particular
cases of it have recently been reproved in numerous papers.

Besides achieving the initial goal, we also achieved several advances in developing the
theory of the "A, B Arino models" with demography, waning immunity, and vaccinations.
One of our main contributions is introducing (or drawing attention to) this class of models,
which we view as new, as our inspiration paper [1] does not allow for changing populations,
waning immunity, and vaccinations.

**Remark 1.** Some confusion took place due to the fact that this paper, which was first submitted on 6 December 2021 to the Journal of Mathematical Biology (see also [2]), underwent a one-year long review process (starting with minor revisions, and ending up with a rejection based mainly on non-standard notations). Naturally, we continued our work on this fascinating model, producing already-published papers [3,4] that discuss examples of the material presented below, and creating—at a first rapid reading—the wrong impressions that preceded the current paper. On the positive side, we are now able to include updates connecting more recent works that tie the A, B Arino models into important works [5–9] on the renewal model.

We will now proceed with a brief historical overview.

**Deterministic mathematical epidemiology.** Deterministic mathematical models have been widely adopted by epidemiologists to model the spread of diseases worldwide, starting with the Bombay plague in 1905–1906 [10], followed by measles, smallpox, chickenpox, mumps, typhoid fever, influenza, and diphtheria (see for example [11]; see [12–14] for comprehensive textbooks). More recently, they have been used for COVID-19 research—see for example [15–27].

**An explanation of our title.** The first foundational result in mathematical epidemiology is the "*R* alternative". As already encountered in the celebrated paper, "A contribution to the mathematical theory of epidemics" [10], its final version states that the "basic replacement number" *R*, which determines the stability of the "disease-free equilibrium", may be computed under certain conditions (see [28,29] and Section 4.1), using only a subset of the ODE dynamics, sometimes called the "infectious equations". The final formula identifies *R* as the spectral radius of a matrix determined by breaking the infectious equations into two parts, one involving a constant "transfer matrix" *A*, and another "new infections matrix" [30]. The [1] paper that motivated our research focuses on the case of a single susceptible (new infection-producing) compartment. This allows introducing a second constant matrix *B*, which is part of the new infection matrix with respect to *S*. Under the further assumption that *B* has rank one, [1] shows that *R* has an explicit simple expression in terms of the two matrices *A*, *B*.

This intriguing result raises several questions. One that is dealt with here is the extension of the [1] result to the varying population case. Another one, as mentioned in the Conclusions section, involves the extension when *B* has a rank greater than one. Finally, one that we hope to deal with in the future involves proving that other fundamental quantities may be simply expressed in terms of the two matrices (*A*, *B*) as suggested by references [6,8].

References [5,6,8], i.e., "The Kermack-McKendrick epidemic model revisited" and "On the formulation of epidemic models (an appraisal of Kermack and McKendrick)" note
that even experienced experts often believe that [10] is just about the SIR system, when in fact this example is just a very special case of a much more general “renewal model” (RE), also called the DDE (delay differential equation) model or KM (Kermack–McKendrick) model. Subsequently, the second paper shows that the behaviors of all infectious compartments may be dealt with by solving just one scalar renewal equation (in the cases of one-strain diseases). This may be achieved in three cases that are easier to treat separately: (a) homogeneous constant population models, (b) varying population models, and (c) models including waning immunity. The relations between the constant population SIR and SEIR models and their renewal function formulations are spelled out at the bottom of page 3 in [6]; these authors then assigned an extension as an exercise to the readers: “we leave it to our readers to elaborate this”; the general idea that delay equations with kernels defined in terms of matrix exponentials correspond to systems of ODEs, called “Linear Chain trick”. It is our impression that this “exercise” has remained unsolved to date in one direction: that of restricting “Arino *A*, *B* models” of rank one, and showing that the row vectors of infectious classes satisfy a renewal equation, with explicit kernels involving *A*, *B*. We provide the solution in (7), following Example 2.

The rest of our title is rather debatable. Our inspiration, the [1] paper, has five very famous authors: Arino, Brauer, van den Driessche, Watmough, and Wu, so it would have been more correct to call this the ABvdDWW model. Arino himself calls this a simple model [31], but we believe that “a deceivingly simple, but extremely challenging model” might have been more appropriate. So, after ruling out the previous two possibilities, we chose to retain only the name of the first two authors of [1]. Note also that in the previous paper [32] (where demography was not taken into account) we used the name “matrix SIR Arino models”; however, we prefer the new name, where we emphasize that there are two matrices (*A*, *B*) involved, and we added the name of Brauer, who died recently (accessed on 10 October 2022). https://math.wisc.edu/2021/10/18/in-memoriam-fred-brauer/.

**Classification of the deterministic epidemic literature.** The deterministic epidemic literature may be divided into three streams.

1. Closed models, with no demography.

2. Models in which the total population may be kept constant, by balancing births and deaths. These models are easier to study, often admit a single endemic equilibrium point, and typically obey the “*R* alternative”—see [33] for a review. However, as death is an essential factor of epidemics, the assumption of a constant population size (clearly an approximation that holds in the short term or for very-large populations) deserves a comment

3. Finally, we arrive at the object of our paper: models with varying total populations. This introduces several challenges, the first being that varying populations typically lead to multiple endemic equilibrium points. This literature reveals the possibility of bi-stability when *R* < 1 (absent from the previous models), even in simple examples [34–36] (thus, for an initial number of "infectives" being high, the trajectory may lie in the basin of attraction of a stable endemic point instead of being eradicated, a discrepancy from what is expected, which suggests that in this range, the deterministic model may be inappropriate). Despite further remarkable works—see, for example, [37–42], the literature studies on models with varying total populations, unlike the two preceding streams, have not yet reached general results. Our paper seems to be the first exception.

**Contributions.** The contribution of this work is three-fold. First, we introduce a class of compartmental models for epidemics that account for varying population sizes, vaccinations, and immunity losses. Secondly, we recommend studying these models at three levels, which may also be viewed as subcases, i.e., (i) the varying population epidemic models, (ii) the “first order” approximation (FA), which is recovered by ignoring certain quadratic terms from the general model proposed (induced by studying the proportions), and which is well-studied in the literature, and (iii) the “intermediate approximation” (IA), introduced here for the first time, and obtained by neglecting only the terms that
are quadratic with respect to the disease/infectious compartments. Finally, we study the equilibrium points of the proposed models and their stability properties.

**Organization.** This paper is organized as follows. In Section 2, we introduce the Arino-Brauer A, B matrix models with demography, vaccination, and waning immunity. Section 3 computes the ages of the infection kernels of these models, related to references [5–9]. Section 4 presents the stability results for the SIR–PH model—see Propositions 2–4. Section 5 provides conclusions and suggestions for future work. Finally, in Appendix A, we present the background on deterministic epidemic models, such as the basic replacement number and the next-generation matrix method.

## 2. Arino-Brauer A, B Matrix Models with Demography, Waning Immunity, Vaccination, and One Susceptible Class

While $\mathcal{R}$ can be defined as the spectral radius of a certain matrix, provided that the next-generation matrix assumptions apply, it is sometimes possible to employ models where $\mathcal{R}$ may be explicitly expressed in terms of the matrices that define the model [1,43–45].

The idea behind these models is to further divide the non-infected compartments into (susceptible)/input classes, defined by producing “new non-linear infections”, and output R classes (e.g., $D, D_e$ in our example), which are fully determined by the rest, and may, therefore, be omitted from the dynamics. Furthermore, It is convenient to restrict to epidemic models the linear forces of infection, as non-linear forces of infection may lead to very complex dynamics [46–49], which are not always easy to interpret epidemiologically. This is in contrast with the Arino models already studied where one may typically establish the absence of periodic solutions (closed orbits, homoclinic loops, and oriented phase polygons) [40,42].

Below, we study the equilibrium points and the dynamical behaviors of such models, restricting to the case of **one susceptible class**. Our goal is to understand the effects of demography, vaccination, and waning immunity, which are missing in the original paper [1].

**Definition 1.** An Arino-Brauer A, B matrix model of type $(1, n, p)$, with demography parameters $(\Lambda, \mu)$ (scalsars), waning immunity column vector $\gamma_r$ and vaccination row vector $\tilde{\gamma}_s$, is characterized by a set of parameters $(A, B, W, v, v_r)$, where $v, v_r$ (in boldface) denote the column vectors of extra death rates, and $A, B, W$ denote respectively matrices of dimensions $n \times n$, $n \times n$ and $n \times p$, respectively. This model contains one susceptible class S, a p-dimensional vector of removed states $\tilde{R}$ (healthy, dead, vaccinated, etc.), and an n-dimensional vector of “disease” state $\tilde{I}$ (latent/exposed, infective, asymptomatic, etc). The dynamics are:

$$ \tilde{I}'(t) = \tilde{I}(t) \left[ \frac{S(t)}{N} B - V \right] - V = A - \text{Diag}(v + \mu \mathbf{1}), \quad (1) $$

$$ S'(t) = -\frac{S(t)}{N}\tilde{I}(t)\mathbf{b} + \Lambda N - (\tilde{\gamma}_s\mathbf{1} + \mu)S(t) + \tilde{R}(t)\gamma_r, \quad \mathbf{b} = B\mathbf{1}, $$

$$ \tilde{R}'(t) = \tilde{I}(t)W + \tilde{\gamma}_s S(t) - \tilde{R}(t)(\text{Diag}(\gamma_r + v_r + \mu\mathbf{1})), $$

$$ N'(t) = S'(t) + \tilde{I}'(t)\mathbf{1} + \tilde{R}'(t)\mathbf{1} = (\Lambda - \mu)N - \tilde{I}(t)v - \tilde{R}(t)v_r, $$

$$ D'(t) = \mu(S(t) + \tilde{I}(t)\mathbf{1} + \tilde{R}(t)\mathbf{1}), $$

$$ D_e'(t) = \tilde{I}(t)v + \tilde{R}(t)v_r. $$

Here,

1. $\tilde{I}(t) \in \mathbb{R}_+^n$ is a row vector whose components model a set of disease states/classes.

2. $S(t) \in \mathbb{R}_+$ is the set of individuals susceptible to being infected.

3. $\tilde{R}(t) \in \mathbb{R}_+^p$ is a row vector whose components model a set of recovered states (or classes), each accounting for individuals who recover from the infection. In what follows, we will focus on the case of one recovered class.
4. B is an $n \times n$ matrix, where each entry $B_{i,j}$ represents the force of infection of the disease class $i$ onto class $j$. We will denote by $b$ the vector containing the sum of the entries in each row of B, namely, $b = B1$.

5. A is an $n \times n$ Markovian sub-generator matrix (i.e., a Markovian generator matrix for which the sum of at least one row is strictly negative), where each off-diagonal entry $A_{i,j}$, $i \neq j$, satisfies $A_{i,j} \ge 0$ and describes the rate of transition from disease class $i$ to disease class $j$; while each diagonal entry $A_{i,i}$ satisfies $A_{i,i} \le 0$ and describes the rate at which individuals in the disease class $i$ leave toward non-infectious compartments. Alternatively, $-A$ is a non-singular M-matrix [1,50]. (An M-matrix is a real matrix $V$ with $v_{ij} \le 0, \forall i \neq j$, having eigenvalues whose real parts are non-negative [51]).

6. $v \in \mathbb{R}_+^n, v_r \in \mathbb{R}_+^p$ are column vectors, giving the death rates in the disease and recovered compartments caused by the epidemic (and possible vaccinations), respectively.

7. $\gamma_r$ is a column vector describing the rates at which individuals lose immunity (i.e., transition from recovered states to susceptible states).

8. $\tilde{\gamma}_s$ is a vector describing the rates at which individuals are vaccinated (immunized).

9. W is an $n \times p$ matrix whose entries model the rates at which individuals in the disease states transfer to recovered or dead states. In what follows, we assume that the $n \times (n+p)$ matrix $\tilde{A} = (A \ W)$ satisfies $\tilde{A}1 = 0$ (namely, the sum of the entries in each row is equal to 0), which implies mass conservation.

**Remark 2.**

1. Note the factorization of the first equation for the diseased compartments $i$, which must appear in any epidemic model, to ensure the existence of a fixed point where these compartments vanish. Furthermore, this equation has two conceptual parts: the first, sometimes called "new infections", contains all interactions of $i$ with other compartments, and the second, $-V$, contains all of the remaining terms.

2. Note that when $p=1$, then $W$ is a column vector, which completes the matrix with negative row sums $A$ to a matrix with zero row sums. Therefore, $W = (-A)1 := a$, where the last notation is standard in the theory of phase-type distributions.

3. A particular but revealing case is that when $p=1$ and matrix $B$ has rank 1, the form $B = b\bar{\alpha}$, where $\bar{\alpha} \in \mathbb{R}^n$ is a probability row vector whose components $a_j$ represent the fractions of susceptibilities entering into the disease compartment $j$, when infection occurs. We will call this SIR-PH, following Riano [50], who emphasized its probabilistic interpretation, i.e., a SIR model where the exponential infection time is replaced by a PH-type distribution. (Some authors, including [52], similarly replace the exponential latency time in class E of SEIR by a PH-type distribution, but this is conceptually unnecessary, since all of the disease classes, may be grouped together in one group, whose phase-type will be determined by those of the components (via Kronecker product operations).). See also [53] for such models, and see Section 3 below on recent connections to the class of non-Markovian epidemic models.
However, the SIR-PH model in previous works precluded important interactions between $S$ and $R$, such as waning immunity and vaccination, and we amended its definition to include these interactions.

**Remark 3.** We did not find any work in the literature on models with linear birth rates at this level of generality. The majority of the literature studies are dedicated to models that may be formally obtained from (1) by letting $N = 1$ (the idea is that $N$ is approximately constant because it is huge or it is observed only over a short period of time). We will call this model "classic/pedagogical", where we add the last qualifier to emphasize that it is unrelated to the model we studied. This is in contrast with the FA and IA approximations introduced later, which approximate the scaled version of (1) introduced below.
It is convenient to reformulate (1) in terms of the fractions normalized by the total population

$$s = \frac{S}{N'}, \qquad \vec{i} = \frac{1}{N}\vec{I}, \qquad \vec{r} = \frac{1}{N}\vec{R}, \qquad N = s + \vec{i} \cdot \mathbf{1} + \vec{r} \cdot \mathbf{1}. \qquad (2)$$

The reader may check that the following equations hold for the scaled variables:

$$\left\{
\begin{aligned}
\vec{i}'(t) &= \vec{i}(t) \left( \left\{ s(t) B + \left( \vec{i}(t) \mathbf{v} + \vec{r}(t) \mathbf{v}_r \right) I_n \right\} + \{ A - \mathrm{Diag}[\mathbf{v} + \Lambda \mathbf{1}] \} \right) \\
s'(t) &= \Lambda - (\Lambda + \gamma_s)s(t) + \vec{r}(t)\gamma_r - s(t) \vec{i}(t)(\mathbf{b} - \mathbf{v}) + s(t) \vec{r}(t)\mathbf{v}_r \\
\vec{r}'(t) &= \vec{i}(t)W + s(t)\vec{\gamma}_s - \vec{r}(t) \left( \mathrm{Diag}[\gamma_r + \mathbf{v}_r + \Lambda \mathbf{1}] - (\vec{i}(t)\mathbf{v} + \vec{r}(t)\mathbf{v}_r)I_p \right) \\
s(t) + \vec{i}(t)\mathbf{1} + \vec{r}(t)\mathbf{1} &= 1
\end{aligned}
\right.
, \quad (3)$$

where we let $\gamma_s := \vec{\gamma}_s \mathbf{1}$. Moreover, by letting $n := s + \vec{i} \mathbf{1} + \vec{r} \mathbf{1}$, we have

$$n'(t) = (\Lambda - \vec{i}(t)\mathbf{v})(1 - n(t)) = 0.$$

Hence, the above equation guarantees that if $s(t_0) + i(t_0) + r(t_0) = 1$ for some $t_0 \in \mathbb{R}_{\ge 0}$, then $s(t) + i(t) + r(t) = 1$ for all $t \ge t_0$. Accordingly, in what follows, we will always assume that $n(t_0) = 1$, which guarantees that $n(t) = 1, \forall t$.

The following definition places (in a common framework) the dynamics of the scaled process and two interesting approximations.

**Definition 2.** Let $\Phi_s, \Phi_i, \Phi_r \in \{0, 1\}$ and let

$$\begin{aligned}
s'(t) &= \Lambda - (\Lambda + \gamma_s)s(t) + \vec{r}(t)\gamma_r - s(t)\vec{i}(t)\mathbf{b} + \Phi_s s(t)(\vec{i}(t)\mathbf{v} + \vec{r}(t)\mathbf{v}_r), \\
\vec{i}'(t) &= \vec{i}(t)\left(s(t)B + A - \mathrm{Diag}[\mathbf{v} + \Lambda \mathbf{1}]\right) + \Phi_i \vec{i}(t)\left(\vec{i}(t)\mathbf{v} + \vec{r}(t)\mathbf{v}_r\right), \\
\vec{r}'(t) &= s(t)\vec{\gamma}_s + \vec{i}(t)W - \vec{r}(t)\mathrm{Diag}[\gamma_r + \mathbf{v}_r + \Lambda \mathbf{1}] + \Phi_r \vec{r}(t)(\vec{i}(t)\mathbf{v} + \vec{r}(t)\mathbf{v}_r), \\
s(t) + \vec{i}(t)\mathbf{1} + \vec{r}(t)\mathbf{1} &= 1.
\end{aligned} \qquad (4)$$

1. The model (4) with $\Phi_s = \Phi_i = \Phi_r = 1$ will be called a scaled model (SM).

2. The model (4) with $\Phi_s = \Phi_i = \Phi_r = 0$ will be called the first approximation (FA).

3. The model (4) with $\Phi_s = \Phi_r = 1$ and $\Phi_i = 0$ will be called the intermediate approximation (IA).

**Example 1.** The SIR-PH-FA models [32,50] (in which there is only one input class s and one output class r) are defined by:

$$\left\{
\begin{aligned}
\vec{i}'(t) &= \vec{i}(t) \left( s(t) B + A - \mathrm{Diag}[\mathbf{v} + \Lambda \mathbf{1}] \right) \\
s'(t) &= -s(t) \vec{i}(t) \mathbf{b} + \Lambda - (\Lambda + \gamma_s)s(t) + r(t)\gamma_r \\
r'(t) &= \vec{i}(t) \mathbf{a} + s(t)\gamma_s - r(t)(\gamma_r + \mathbf{v}_r + \Lambda).
\end{aligned}
\right.
\qquad (5)$$

**Example 2.** The FA-version of the SEIRS model, with $\Phi_s = \Phi_i = \Phi_r = 0$, is given by:

$$\left\{
\begin{aligned}
(e'(t), i'(t)) &= (e(t), i(t)) \begin{pmatrix} -(\gamma_e + \Lambda) & \gamma_e \\ \beta_s(t) & -(\gamma + \Lambda + \nu) \end{pmatrix} \\
s'(t) &= \Lambda - s(t)(\beta i(t) + \gamma_s + \Lambda) + \gamma_r r(t) \\
r'(t) &= \gamma i(t) + \gamma_s s(t) - r(t)(\gamma_r + \mathbf{v}_r + \Lambda)
\end{aligned}
\right.
\qquad (6)$$

is a particular case of the A, B matrix Arino model obtained when
$$A = \begin{pmatrix} -\gamma_e & \gamma_e \\ 0 & -\gamma \end{pmatrix}, B = \begin{pmatrix} 0 & 0 \\ \beta & 0 \end{pmatrix} = b\vec{\alpha}, b = \begin{pmatrix} 0 \\ \beta \end{pmatrix}, \vec{\alpha} = (1,0), \nu = \begin{pmatrix} 0 \\ \nu \end{pmatrix}, W = a = (-A)\mathbf{1} = \begin{pmatrix} 0 \\ \gamma \end{pmatrix}.$$

See Figure 1.

**Remark 4.** In this paper, we wrote the dynamical systems with the vector state variables pre-multiplying the model matrices, diverging from the dynamical systems tradition. This is because we want $\vec{\alpha}$ defined by the factorization of B (which will play an important role below) to be a row vector, following the convention in the theory of phase-type distributions.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Chart flow of the SEIRS model. The force of infection is $F_o = \beta i$. The red edge corresponds to the entrance of susceptibilities into the disease classes, the dashed green edges correspond to the contacts between the diseases and susceptibilities, the brown edge is the rate of the transition matrix $V$, and the cyan dashed line corresponds to the rate of waning immunity. The remaining black lines correspond to the inputs and outputs of the birth and natural death rates, respectively.</figcaption></figure>

Figure 2 below makes the case where epidemic models can be studied at three levels (the exact model, and its FA and IA approximations). The figure depicts the endemic points of an SIR-type example of the (1) model, as discussed in detail in [3].

**Remark 5.** We end this section by pointing out that our scaled model (3) belongs to a class of models introduced, for different motivations, in [54]. Indeed, after dropping the first equation

$$s'(t) = \Lambda - s(t)\vec{i}(t)b - (\Lambda + \gamma_s)s(t) + \vec{r}(t)\gamma_r + \begin{pmatrix} \vec{i}(t)v + \vec{r}(t)v_r \end{pmatrix}s(t)$$

and rewriting the rest as:

$$\begin{pmatrix} \vec{i}'(t) & \vec{r}'(t) \end{pmatrix} = (0 \quad s(t)\gamma_r) + \begin{pmatrix} \vec{i}(t) & \vec{r}(t) \end{pmatrix} \times \begin{bmatrix} s(t)B + A - (\Lambda + \nu)I_n & W \\ \mathbf{0} & -\text{Diag}(\gamma_r) + \Lambda I_p \end{bmatrix} - \text{Diag}(v_x) + \vec{x}\text{Diag}(v_x)\text{Diag}(\vec{x}),$$
where $\vec{x} = (\vec{i}, \vec{r})$, $\nu_x = (\nu, \nu_r)$, we recognize the above as a particular case of ([54], 2.1).

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Parametric $(s, i)$ plots of the scaled epidemic and its FA and intermediate approximations for a SIR-type model with one infectious class, starting from a starting point SP with $i_0 = 10^{-6}$, with $\mathcal{R}_0 = 3.21, \beta = 5, \gamma = 1/2, \Lambda = \mu = 1/10, \gamma_r = 1/6, \gamma_s = 0.01, \nu_i = 0.9, \nu_r = 0$. EESc, EEIn, and EEF0A are the stable endemic points of the scaled model, IA (intermediate) model, and FA model, respectively. The green vertical line denotes the immunity threshold $1/R = s_{EEF0A} = s_{EEIn}$. Note that the epidemic will at first spend a long time (since births and deaths have slow rates as compared to the disease) in the vicinity of the manifold $\vec{i}(t) = 0$, where the three processes are indistinguishable, before turning toward the endemic equilibrium point(s).</figcaption></figure>

**3. The Semi-Groups and Age of Infection Kernels Associated with SIR-PH-FA Models with B of Rank One**

In the case when B has rank one, SIR-PH-FA models have an associated explicit age of infection kernel, which allows obtaining $R_0$ via integrating.

Our attention to this subject was drawn by formulas from [6] (p. 3) for the “renewal/age of infection kernels” for particular cases of the SIR and SEIR models. These authors (as an exercise) extended their formulas to other models; determining which models to extend to was part of the exercise. Six years later, the exercise was first solved by Champredon-Dushoff-Earn [8] for the Erlang-SEIRS models. We provide below a further extension to the case of SIR-PH-FA models with B of rank one—see also [7,9] for related results.

**Proposition 1.** For a SIR-PH-FA model with $B = b\vec{\alpha}$ of rank one, the formula for the age of infection/renewal kernel (see [5–9] for expositions of this concept), is

$$a(\tau) = \vec{\alpha}e^{-\tau V}b, \qquad (7)$$

where $-V = A - (\text{Diag}[\delta + \Lambda\mathbf{1}])$ (it may be checked that this fits the formula on page 3 of [6] for SEIR when $\Lambda = 0, \delta = 0$).

**Remark 6.** The beautiful decomposition (7) explicitizes the age of infection kernel $a(\tau)$ as a product of the three essential factors of an epidemic: the initial proportions of infectious $\vec{\alpha}$, the “transfer semi-group” $e^{-tV}$, and the infectivity vector $b$. This suggests developing calibration approaches that are not model-dependent, but choose instead themselves the order of V and, subsequently, the three factors.

**Proof.** Let us transform the infectious equations into an integral equation by applying the variation of the constant formula. □
The first step, the solution of the homogeneous part of the infectious equations, yields

$$ \vec{\Gamma}'(t) = -\vec{\Gamma}(t)V \text{ is } \vec{\Gamma}(t) = \vec{\Gamma}(0)e^{t(-V)}. \quad (8) $$

**Remark 7.** When $\vec{\Gamma}(0)$ is a probability vector, this solution has the interesting probabilistic interpretation of the survival probabilities in the various components of the semigroup generated by the Metzler/Markovian sub-generator matrix $-V$ (which inherits this property from the phase-type generator A). Practically, $\vec{\Gamma}(t)$ will give the expected fractions of individuals who are still in each compartment at time t.

The variation of the constant formula implies that $\vec{i}(t)$ satisfies the integral equation:

$$ \vec{i}(t) = \vec{i}(0)e^{-tV} + \int_{0}^{t} s(\tau) \vec{i}(\tau) Be^{-(t-\tau)V} d\tau. \quad (9) $$

We relate now this result to the "constituent equation" of the differential delay/renewal model Theorem 2.9 in reference [8]. In the case $B = b\vec{\alpha}$, (9) becomes

$$ \vec{i}(t) = \vec{i}(0)e^{-tV} + \int_{0}^{t} s(\tau) \vec{i}(\tau) b\vec{\alpha} e^{-(t-\tau)V} d\tau. \quad (10) $$

Denoting further by $\tilde{i}(t) = \vec{i}(t)b$ the total force of infection, and putting

$$ a(t) = \vec{\alpha}e^{-tV}b, $$

our system may also be written as a system of two scalar equations

$$ \begin{cases} s'(t) = \Lambda - (\Lambda + \gamma_s)s(t) - s(t)\tilde{i}(t) \\ \tilde{i}(t) = \vec{i}(0)e^{-tV}b + \int_0^t s(\tau)\tilde{i}(\tau)a(t-\tau)d\tau = \int_{-\infty}^t s(\tau)\tilde{i}(\tau)a(t-\tau)d\tau, \end{cases} \quad (11) $$

where the second equality holds if in (9) it holds that $\vec{i}(0) = k\vec{\alpha}$, and if $\tilde{i}(\tau)$ on the interval $(-\infty, 0]$ is $\frac{k}{s_0}\delta_0(\tau)$, where $\delta_0(\tau)$ denotes the generalized Dirac function—compare to Theorem 2.8 in reference [8]. This second form is related to Theorem 2.7b in reference [8] and Theorem 1 in reference [6]. (In fact, these authors work with the related incidence flux between the s and $\vec{i}$ variables Incid := $s\vec{i}b$, denoted by $i(t)$ in [8], and by $F(t)$ in [6]).

**Remark 8.**

1. Note that (11) is a SI system with $\tilde{i}$ satisfying a distributed delay (DD) equation. Such equations have been a constant preoccupation in mathematical epidemiology since the founding paper [10], see for example [55]. This is quite natural since susceptibilities become infectious only after a time $\tau_0$ after their contact. If $\tau_0$ is assumed to be a known constant, then the second equation of the SI model would involve a Dirac kernel $a(\tau) = k\delta_{\tau_0}(\tau)$; however, since this is not the case, it is natural to replace the Dirac kernel by a continuous one. The fact that DD systems can be approximated by ODE systems has long been exploited in epidemic literature, under the name of "linear chain trick" (which is called Erlangization in queuing theory and mathematical finance)—see [7,9,56] for recent contributions and further references. The opposite direction, i.e., the solution of the exercise in [6] of identifying the kernels associated with ODE models, was not resolved in this generality, prior to our paper.

2. Rewriting the infectious equations using variations of constants transformed the "ODE/Markovian" version (5) into a particular case of a non-local/renewal SI system (11), with the initial condition given by Dirac's delta, and with the kernel of the matrix-exponential type.
3. As well-known, for DD models, the kernel may be factored as the product of $\mathcal{R}$ with the density of the age of infection of the "intrinsic generating interval" [8,9,57]. $\mathcal{R}$ may, thus, be obtained by integrating the kernel

$$ \mathcal{R} = \int_{0}^{\infty} a(\tau) d\tau = \int_{0}^{\infty} \vec{\Gamma}(\tau) \mathbf{b} d\tau = \vec{\alpha} V^{-1} \mathbf{b} \quad (12) $$

with $\Gamma(0) = \vec{\alpha}$—see also (2.3) in reference [8] and (5.9) in reference [7]. A second direct proof of (12), without appealing to the first equality above (note that the kernel $a(t) = \vec{\Gamma}(t)\mathbf{b}$ may be interpreted as an average time-dependent total infectivity; the fact that the disease-free equilibrium stability threshold has an integral representation that was first proved in [30]) is provided in Proposition 2 below.

4. Finally, the "intrinsic generation-interval density" is $g(t) = \frac{a(t)}{\mathcal{R}} = \frac{\vec{\Gamma}(t)\mathbf{b}}{\mathcal{R}}$—see (2.6) in reference [8].

### The Eigenstructure of the Jacobian for the Scaled Model (3)

Using $\nabla \vec{i}(t)\vec{v} = \vec{v} \vec{i}(t) + \vec{i}(t)\vec{v} I_n$, and putting $x = \vec{i}v + \vec{r}v_r$, we find that the transpose of the Jacobian matrix of the scaled model (3) is given by

$$ J = \begin{pmatrix} x - \vec{i} \mathbf{b} - (\Lambda + \gamma_s) & \vec{i} \mathbf{B} & \vec{r}_s \\ \mathbf{s}(v - \mathbf{b}) & \mathbf{s} \mathbf{B} + A - \text{Diag}(v + \Lambda \mathbf{1}) + \vec{v} \vec{i} + (\vec{i} v + \vec{r} v_r) I_n & W + v \vec{r} \\ \gamma_r + s v_r & \vec{v}_r \vec{i} & \vec{v}_r - \text{Diag}[\gamma_r + v_r + \Lambda] + (\vec{i} v + \vec{r} v_r) I_p \end{pmatrix}. $$

At the disease-free equilibrium ($s_{dfe}, \vec{0}, \vec{r}_{dfe}$), we have

$$ J_{dfe} = \begin{pmatrix} -(\Lambda + \gamma_s) + \vec{r}_{dfe} v_r & 0 & \vec{r}_s \\ s_{dfe}(v - \mathbf{b}) & s_{dfe} \mathbf{B} + A - \text{Diag}(v + \Lambda \mathbf{1}) + \vec{r}_{dfe} v_r I_n & W + v \vec{r}_{dfe} \\ \gamma_r + s_{dfe} v_r & 0 & \vec{v}_r \vec{r}_{dfe} - \text{Diag}[\gamma_r + v_r + \Lambda] + (\vec{r}_{dfe} v_r) I_p \end{pmatrix}. $$

**Remark 9.** Note the structure of the second column, which is equivalent to the factorization property, and implies Proposition 2, below.

## 4. Stability Results for the SIR-PH Model

We consider here the particular case with a single recovered class, called the SIR-PH model, where $W = \mathbf{a} = (-A)\mathbf{1}$. Notice that, in this case, we have a reduced set of parameters ($\Lambda = \mu, \vec{\alpha}, A, \mathbf{b}, v, v_r, \gamma_s, \gamma_r$).

### 4.1. The Basic Replacement Number for SIR-PH via the Next-Generation Matrix Method [28,30]

We follow up here on a remark preceding Theorem 2.1 in reference [1], and show in the following proposition that the simplified formula for the basic replacement number still holds when waning immunity and vaccination are allowed, provided that $p = 1$ and $B = b\vec{\alpha}$ has rank one.

**Proposition 2.** Consider a SIR-PH model (i.e., a single-recovered class), with parameters ($\Lambda = \mu, \vec{\alpha}, A, \mathbf{b}, v, v_r, \gamma_s, \gamma_r$), and matrix of recovery rates $W = \mathbf{a} = (-A)\mathbf{1}$.

1. (a) When $v_r = 0$, the unique DFE is $(\frac{\Lambda + \gamma_r}{\Lambda + \gamma_r + \gamma_s}, 0, \frac{\gamma_s}{\Lambda + \gamma_r + \gamma_s})$.

(b) When $v_r > 0$, exclude the case $\Lambda = \gamma_r = 0$. Then, there exists a unique DFE ($s_{dfe}, 0, 1 - s_{dfe}) \in \mathcal{D}$, where $s_{dfe}$ satisfies the second-order equation $\mathbf{s}[v_r \mathbf{s} + \Lambda + \gamma_r + \gamma_s - v_r] - (\Lambda + \gamma_r) = 0$ and is given by

$$ s_{dfe} = \frac{\sqrt{4v_r(\Lambda + \gamma_r) + (\Lambda + \gamma_r + \gamma_s - v_r)^2 - (\Lambda + \gamma_r + \gamma_s - v_r)}}{2v_r} \quad (13) $$

2. The weak $\mathcal{R}$ alternative holds for the threshold parameter,

$$ \mathcal{R}_0 = \lambda_{PF}(FV^{-1}), \quad (14) $$
where F, V are defined in (21), and $\lambda_{PF}$ denotes the (dominant) Perron-Frobenius eigenvalue.

3. For rank 1, $B := \mathbf{b}\bar{\alpha}$ and $v_r = 0$, we further have

(a)

$$ \mathcal{R}_0 = s_{dfe} \mathcal{R}, \text{ where } \mathcal{R} = \bar{\alpha} V^{-1} \mathbf{b}. \qquad (15) $$

(b) If $\mathcal{R}_0 \le 1$, and if the perturbation from linearity defined in (A3) is nonnegative, then the scalar combination

$$ Y = \vec{i} V^{-1} \mathbf{b} \qquad (16) $$

is a Lyapunov function for the DFE.

### Proof.

1. The disease-free system (with $i = 0, r = 1-s$) reduces to

$$ s'(t) = \Lambda - (\gamma_s + \Lambda)s(t) + (\gamma_r + v_r s(t))(1 - s(t)). \qquad (17) $$

For the fixed points, depending on $v_r$, we must solve either a quadratic or a linear equation

$$ \begin{cases} \Lambda + \gamma_r - s[v_r s + \Lambda + \gamma_r + \gamma_s - v_r] = 0 & v_r > 0 \\ s[\Lambda + \gamma_r + \gamma_s] - (\Lambda + \gamma_r) = 0 & v_r = 0 \end{cases} $$

One root

$$ s_{dfe} = \begin{cases} \frac{\Lambda + \gamma_r}{\Lambda + \gamma_r + \gamma_s} & v_r = 0 \\ \frac{\sqrt{\Delta_{dfe} - (\Lambda + \gamma_r + \gamma_s - v_r)}}{2v_r}, & \Delta_{dfe} = 4v_r(\Lambda + \gamma_r) + (\Lambda + \gamma_r + \gamma_s - v_r)^2 \quad v_r > 0 \end{cases} \qquad (18) $$

is always in $[0, 1]$ and will be denoted by $s_{dfe}$.

The other root in the quadratic case $v_r > 0$

$$ \frac{v_r - (\Lambda + \gamma_r + \gamma_s) - \sqrt{4v_r(\Lambda + \gamma_r) + (\Lambda + \gamma_r + \gamma_s - v_r)^2}}{2v_r} \qquad (19) $$

is strictly negative, unless

$$ \begin{cases} \Lambda + \gamma_r = 0 \\ v_r \ge \gamma_s + \Lambda + \gamma_r \end{cases} \Leftrightarrow \begin{cases} \Lambda = \gamma_r = 0 \\ v_r \ge \gamma_s \end{cases}, \qquad (20) $$

in which case it yields a second DFE point with $s = 0 = i$, which we exclude (in order to apply the next-generation matrix method).

2. It is enough to show here the conditions of Theorem 2 in reference [29] hold, with respect to the infectious set $\vec{i}$, and a certain splitting.

The DFE and its local stability for the disease-free system are easy to check.

We provide a splitting for the infectious equations:

$$ \vec{i}'(t) = \vec{i}(t) \left[ \mathbf{s}(t)B + \vec{i}(t)\mathbf{v}I_n + v_r r(t)I_n \right] - \vec{i}(t)[\mathrm{Diag}(\mathbf{v} + \Lambda\mathbf{1}) - A] := \mathcal{F}(s, \vec{i}) - \mathcal{V}(\vec{i}) $$

(where $r = 1 - s - \vec{i} \mathbf{1}$). The corresponding gradients at the DFE $\vec{i} = 0$ are

$$ \begin{cases} F = \left[ \frac{\partial \mathcal{F}(\mathbf{X}^{(\mathrm{DFE})})}{\partial \vec{i}} \right] = sB + v_r r I_n, \\ V = \left[ \frac{\partial \mathcal{V}(\mathbf{X}^{(\mathrm{DFE})})}{\partial \vec{i}} \right] = \mathrm{Diag}(\mathbf{v} + \Lambda \mathbf{1}) - A. \end{cases} \qquad (21) $$
We note that $F$ has non-negative elements and that $V$ is an M-matrix; therefore, $V^{-1}$ exists and has non-negative elements, $\forall \Lambda, \nu$. We may check that the conditions (A2) are satisfied.

For example, the last non-negativity condition in (A2)

$$ \vec{i}(t)[\mathrm{Diag}(\nu + \Lambda\mathbf{1}) - A]\mathbf{1} \geq 0, \forall \vec{i} \in \mathcal{D}, \qquad (22) $$

is a consequence of $-A$ being an M-matrix, which implies $-A\mathbf{1} \geq 0$, component-wise.

3. (a) Now if $n = 1$, or if $B = b\bar{\alpha}$ has rank 1, and $v_r = 0$, the matrix $F$ in (21) is the product of a column vector and a row vector, the dimension of its image is one, and the same holds for $BV^{-1}$. Equivalently, $rank(BV^{-1}) = 1$. The "rank-nullity theorem" $rank(BV^{-1}) + \mathrm{nullity}(BV^{-1}) = n$ [58] implies that $(n - 1)$ of the eigenvalues of $BV^{-1}$ are zero, and the Perron-Frobenius eigenvalue is the remaining one. This latter one must be equal to the trace of $BV^{-1} = b\bar{\alpha} V^{-1}$, which may be checked to equal $\bar{\alpha} V^{-1} b$. Finally, the linearity $\mathcal{R}_0 = s_{dfe} \mathcal{R}$ is obvious.

(b) This is a particular case of [59] since the Perron-Frobenius eigenvector in our rank-one case may be taken as $b$.

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

**Remark 10.** $s_{dfe}$ is continuous in $v_r$, since for $v_r$ small,

$$ s_{dfe} \approx \frac{\Lambda + \gamma_r + \gamma_s - v_r + \frac{2v_r(\Lambda + \gamma_r)}{\Lambda + \gamma_r + \gamma_s - v_r} - (\Lambda + \gamma_r + \gamma_s - v_r)}{2v_r} \rightarrow \frac{\Lambda + \gamma_r}{\Lambda + \gamma_r + \gamma_s} $$

(this approximation may be made rigorous by applying the rule of "l'Hospital").

## 4.2. The Classic/Pedagogical SIRS-FA Model

In this section, we present more explicit results for the disease-free equilibrium and the endemic equilibrium of the following model, referred to as SIRS-FA,

$$ \left\{ \begin{array}{l} s'(t) = \Lambda - (\mu + \gamma_s)s(t) + \vec{r}(t)\gamma_r - s(t)\vec{i}(t)\mathbf{b} \\ \vec{i}'(t) = \vec{i}(t)(s(t)\mathbf{b}\bar{\alpha} - V) \\ \vec{r}'(t) = s(t)\vec{\gamma}_s + \vec{i}(t)W - \vec{r}(t)(\mathrm{Diag}[\gamma_r + \nu_r + \mu\mathbf{1}]) \\ s(t) + \vec{i}(t)\mathbf{1} + \vec{r}(t)\mathbf{1} = 1 \end{array} \right. . \qquad (23) $$

**Proposition 3.**

1. The pedagogical ($\Lambda, \mu, A, \vec{\gamma}_s, \gamma_r$) SIRS system (23) has a unique disease-free equilibrium (DFE) fixed point ($s_{dfe}, \bar{0}, \vec{r}_{dfe}$), where

$$ \left\{ \begin{array}{l} s_{dfe} = \frac{\Lambda}{\mu + \gamma_s - \vec{\gamma}_s (\mathrm{Diag}(\gamma_r + \nu_r + \mu \mathbf{1}))^{-1} \gamma_r}, \\ \vec{r}_{dfe} = s_{dfe} \vec{\gamma}_s \mathrm{Diag}(\gamma_r + \nu_r + \mu \mathbf{1})^{-1} \end{array} \right. . $$

In the SIR-PH case, the DFE simplifies to (this formula has appeared already in many particular cases—see for example (19–20) in reference [60]):

$$ (s_{dfe}, \bar{0}, r_{dfe}) = \left( \frac{\Lambda(\mu + \gamma_r + v_r)}{\mu\gamma_r + (\mu + v_r)(\mu + \gamma_s)}, \bar{0}, s_{dfe}\gamma_s \right). \qquad (24) $$

2. If $\mathcal{R} > 1$, then the pedagogical system (23) has a unique second fixed point within its forward-invariant set. This endemic fixed point is such that $1/s_{ee}$ is an eigenvalue of the matrix $BV^{-1}$.
In the SIR-PH case, it must satisfy

$$1/s_{ee} = \mathcal{R} = \vec{\alpha} V^{-1} \mathbf{b}.$$

The disease component $\vec{i}_{ee}$ satisfy:

$$\vec{i}_{ee}(\frac{1}{\mathcal{R}} B + A - \text{Diag}(\mathbf{v} + \mu \mathbf{1})) := \vec{i}_{ee} M = 0 \quad (25)$$

($\vec{i}_{ee}$ is a Perron-Frobenius eigenvector of the matrix $M$ related to the next-generation matrix).

3. The normalization of $i_{ee}$ is given by (30) below. When $\mu = \Lambda$, this becomes:

$$\vec{i} \left( \mathbf{b} - \mathcal{R} \mathbf{a} \frac{\gamma_r}{\gamma_r + \nu_r + \Lambda} \right) = \Lambda \mathcal{R} - \Lambda + \gamma_s \left( \frac{\gamma_r}{\Lambda + \gamma_r + \nu_r} - 1 \right). \quad (26)$$

4. The disease-free equilibrium is locally asymptotically stable if $\mathcal{R} < 1$ and is unstable if $\mathcal{R} > 1$.

5. When $\nu_r = 0$, the critical vaccination defined by solving $\mathcal{R}_0 = 1$ with respect to $\gamma_s$ is given by

$$\gamma_s^* := (\Lambda + \gamma_r)(\vec{\alpha}V^{-1}\mathbf{b} - 1) = (\Lambda + \gamma_r)(\mathcal{R} - 1). \quad (27)$$

We show in Figure 3 a stream plot of the SIRS-FA model that illustrates the above properties.

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Stream plots of $(s, i)$ for an example of a SIRS-FA model with one infectious class, when $\gamma_s \in \{1/100, 3\}$ is smaller and bigger, respectively, than the critical vaccination $\gamma_s^* = 0.239087$ defined in (27).</figcaption></figure>

**Proof.** The proof starts by examining the two cases that arise from factoring in the disease equations. More precisely, we will search separately in the disease-free set $\mathcal{I} := \{\vec{i} = \vec{0}\}$ and in its complement $\mathcal{I}^c$. Then:

1. Either $\vec{i} \in \mathcal{I}$ and solving

$$\begin{cases} \Lambda = (\mu + \gamma_s)s - \vec{r}\gamma_r \\ \vec{0} = s\vec{\gamma}_s - \vec{r}\text{Diag}(\gamma_r + \mathbf{v}_r + \mu\mathbf{1}) \end{cases} = (s \quad \vec{r}) \begin{pmatrix} \mu + \gamma_s & \vec{\gamma}_s \\ -\gamma_r & -\text{Diag}(\gamma_r + \mathbf{v}_r + \mu\mathbf{1}) \end{pmatrix}$$

for $s, \vec{r}$ yields the unique DFE (it may be shown by induction that the determinant is negative). Or,

2. The determinant of the resulting homogeneous linear system for $\vec{i} \neq \vec{0}$ must be 0, which implies that $s = s_{ee}$ satisfies

$$\det[s_{ee} B + A - \text{Diag}(\mathbf{v} + \mu \mathbf{1})] = 0. \quad (28)$$
Note now that $V = \text{Diag}(\nu + \mu\mathbf{1}) - A$ is an invertible matrix. Using $\det(UU') = \det(U)\det(U')$, (28) may be written as

$$ \det[(s_{ee}B - V)V^{-1}] = 0 \Leftrightarrow \det[s_{ee}BV^{-1} - \mathbf{I}] = 0 $$

Dividing then by $s_{ee}$ yields the characteristic polynomial of a matrix

$$ \det[BV^{-1} - \frac{1}{s_{ee}}\mathbf{I}] = 0. \qquad (29) $$

In the SIR-PH case, noting that the rank one matrix $BV^{-1}$ has $(n-1)$ eigenvalues equal to zero, we conclude that the inverse of the susceptible fraction $1/s_{ee}$ of an endemic state must equal the Perron-Frobenius eigenvalue $\mathcal{R}$. Note that $s_{ee} < 1$ follows from our assumption on $\mathcal{R}$. The other coordinates are determined starting with $\vec{i}$, which must be proportional to a Perron-Frobenius nonnegative eigenvector.

3. Recall the system

$$ \left\{ \begin{array}{l} 0 = \Lambda - (\mu + \gamma_s)s + \vec{r}(t)\gamma_r - s y \\ \vec{0} = \vec{i} (s B + A - \text{Diag}(\nu + \mu\mathbf{1})) \\ \vec{0} = \vec{i} W + s\vec{\gamma}_s - \vec{r}(\text{Diag}(\gamma_r + \nu_r + \mu\mathbf{1})) \end{array} \right. . $$

Since $s_{ee} = \frac{1}{\mathcal{R}}$, and $\vec{i}_{ee}$ is known up to the proportionality constant $y = \vec{i}\mathbf{b}$, it only remains to solve the last equation in (30) below:

$$ \left\{ \begin{array}{l} \vec{r} = [\vec{i} W + s\vec{\gamma}_s] \text{Diag}(\gamma_r + \nu_r + \mu\mathbf{1})^{-1} \\ \vec{i} (s\mathbf{b} - W\text{Diag}(\gamma_r + \nu_r + \mu\mathbf{1})^{-1}\gamma_r) = \Lambda - (\mu + \gamma_s)s + s\vec{\gamma}_s \text{Diag}(\gamma_r + \nu_r + \mu\mathbf{1})^{-1}\gamma_r \end{array} \right. \qquad (30) $$

This equation may be solved numerically. When $p = 1 \Rightarrow W := \mathbf{a} = (-A)\mathbf{1}$, the last formula yields

$$ \vec{i} \left( \mathbf{b} - \mathcal{R} \mathbf{a} \frac{\gamma_r}{\gamma_r + \nu_r + \mu} \right) = \Lambda \mathcal{R} - \mu + \gamma_s \left( \frac{\gamma_r}{\gamma_r + \nu_r + \mu} - 1 \right) \qquad (31) $$

When $\mu = \Lambda$, this yields (26).

4. This follows from Theorem 2.1 in reference [59] (it is a consequence of the fact that a linear function proportional to the associated Perron eigenvector is a Lyapunov function when $\mathcal{R} < 1$).

5. The result is immediate by solving $\mathcal{R} = 1$ with respect to $\gamma_s$, where $\mathcal{R}$ is defined in (15).

□

**Remark 11.** Equivalently, it is easy to check that (29) may be reformulated as saying that $s_{ee}$ must satisfy

$$ \lambda_{PF}(s_{ee} B(A - \text{Diag}(\nu + \mu\mathbf{1})))^{-1} = 1 \Leftrightarrow \lambda_{PF}(s_{ee} B + A - \text{Diag}(\nu + \mu\mathbf{1})) = 0, $$

where $\lambda_{PF}$ denotes the Perron-Frobenius eigenvalue. We note that the matrices in the last formula-intervene also in the next-generation matrix approach.
### 4.3. A Glimpse of the Intermediate Approximation Model for Matrix SIRS, with $\Lambda = \mu$
The intermediate approximation associated with (3) is

$$ egin{cases} 
\begin{aligned}
\mathbf{s}'(t) &= \Lambda - (\Lambda + \gamma_s)\mathbf{s}(t) + \vec{r}(t)\gamma_r - \mathbf{s}(t)\vec{i}(t)(\mathbf{b} - \mathbf{v}) + \mathbf{s}(t)\vec{r}(t)\mathbf{v}_r \\
\vec{i}'(t) &= \vec{i}(t)(\mathbf{s}(t)B + A - \text{Diag}[\mathbf{v}] - \Lambda I_n)
\end{aligned} \\
\vec{r}'(t) = \mathbf{s}(t)\vec{\gamma}_s + \vec{i}(t)W - \vec{r}(t)(\text{Diag}[\gamma_r + \mathbf{v}_r + (\Lambda - (\vec{i}(t)\mathbf{v} + \vec{r}(t)\mathbf{v}_r))\mathbf{1}])
\end{cases} 
; \quad (32) $$

when $i = 0$, we have

$$ s_{dfe} = \frac{\Lambda + \vec{r}_{dfe}\gamma_r}{\Lambda + \gamma_s - \vec{r}_{dfe}\mathbf{v}_r}, \quad (33) $$

and from the last equation in (32), $\vec{r}_{dfe}$ satisfies the following third-order equation

$$ (\Lambda + \gamma_s - \vec{r}_{dfe}\mathbf{v}_r) \vec{r}_{dfe} (\text{Diag}[\gamma_r + \mathbf{v}_r + (\Lambda - \vec{r}_{dfe}\mathbf{v}_r)\mathbf{1}]) - \vec{r}_{dfe}\gamma_r\vec{\gamma}_s = \Lambda\vec{\gamma}_s $$

**Proposition 4.** Assume $\mathbf{v}_r = 0$. Then:

(a) The DFE points of the scaled, intermediate approximation, and FA are equal, with

$$ \vec{r}_{dfe} = \Lambda\vec{\gamma}_s ((\Lambda + \gamma_s)\text{Diag}[\gamma_r + \Lambda\mathbf{1}] - \gamma_r\vec{\gamma}_s)^{-1}. $$

(b) In the SIR-PH case with $\vec{r}$ scalar, they all reduce to $(\frac{\Lambda+\gamma_r}{\Lambda+\gamma_r+\gamma_s}, 0, \frac{\gamma_s}{\Lambda+\gamma_r+\gamma_s})$ (18).

(c) The endemic point is unique. It satisfies $s_{ee} = \frac{1}{\mathbb{R}}$, $\vec{i}$ is an eigenvector of the matrix $\frac{1}{\mathbb{R}}B + A - \text{Diag}[\mathbf{v}] - \Lambda I_n$ for the eigenvalue 0, and

$$ \vec{r}_{ee} = (s_{ee}\vec{\gamma}_s + \vec{i}_{ee}W)(\text{Diag}[\gamma_r + (\Lambda - \vec{i}_{ee}\mathbf{v})\mathbf{1}])^{-1} $$

and it satisfies the normalization

$$ \vec{i}(t)(\mathbf{b} - \mathbf{v}) - \mathbb{R}\vec{i}W(\text{Diag}[\gamma_r + (\Lambda - \vec{i}\mathbf{v})\mathbf{1}])^{-1}\gamma_r = \Lambda(\mathbb{R} - 1) - \gamma_s + \vec{\gamma}_s(\text{Diag}[\gamma_r + (\Lambda - \vec{i}\mathbf{v})\mathbf{1}])^{-1}\gamma_r. \quad (34) $$

**Proof.**

(a) The equations determining the three DFEs coincide.

(b) When $\vec{r}$ is scalar, we find

$$ r_{dfe} = \Lambda\gamma_s \frac{1}{(\Lambda + \gamma_s)[\gamma_r + \Lambda] - \gamma_r\gamma_s} = \frac{\gamma_s}{\Lambda + \gamma_s + \gamma_r}. $$

(c) We have

$$ \left\{ \begin{array}{l} 0 = \Lambda - (\Lambda + \gamma_s)s_{ee} + \vec{r}_{ee}\gamma_r - s_{ee}\vec{i}_{ee}(\mathbf{b} - \mathbf{v}), \\ \vec{r}_{ee} = (s_{ee}\vec{\gamma}_s + \vec{i}_{ee}W)(\text{Diag}[\gamma_r + (\Lambda - \vec{i}_{ee}\mathbf{v})\mathbf{1}])^{-1}, \end{array} \right. $$

by substitution, we have

$$ \Lambda(\mathbb{R} - 1) - \gamma_s + \vec{\gamma}_s (\text{Diag}[\gamma_r + (\Lambda - \vec{i}_{ee}\mathbf{v})\mathbf{1}])^{-1} \gamma_r = \vec{i}_{ee}(\mathbf{b} - \mathbf{v}) - \mathbb{R}\vec{i}_{ee}W(\text{Diag}[\gamma_r + (\Lambda - \vec{i}_{ee}\mathbf{v})\mathbf{1}])^{-1} \gamma_r. $$

□

**Remark 12.** When $p = 1$ and $\vec{i}\mathbf{v} = 0$, (34) reduces to (26) when $\mathbf{v}_r = 0$.
## 5. Conclusions and Further Work

We introduced a new class of matrix-parametrized models with varying populations, waning immunity, and vaccinations, and provided a few general results for the particular case of one susceptible class.

By solving the exercise of [6], we found that the $(A, B)$ Arino-Brauer models with the rank of $B$ have a natural “intrinsic age of infection density”. Thus, these deterministic models also have associated stochastic aspects, which reveal themselves when all infectious equations are grouped into one equation.

This raises the question of whether the age of infection kernel may also be associated with

1. $(A, B)$, for which matrix $B$ has a rank bigger than 1;

2. for matrix models involving two or more susceptible classes.
Other directions that are worthy of further work are:

3. Determining the largest domain of attraction of the DFE, in which there exists some linear Lyapunov function that decreases (this might be approachable via linear programming).

4. Calibrating real epidemic data via distributed delay models with non-monotone intrinsic generation-interval densities.

5. Investigating the applicability of the Laplace transform and Laplace-Adomian decomposition approaches to our model, and extensions to fractional models—see [61,62] for possible extension directions.

We believe that studying matrix models, as opposed to individual examples, is crucial for the further development of mathematical epidemiology.

**Author Contributions:** All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by the AB Nexus seed grant from the University of Colorado.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author, F.A, upon reasonable request.

**Acknowledgments:** We thank N. Bacaer for providing the references [63,64], and the reviewers for their thorough reviews and suggestions.

**Conflicts of Interest:** The authors declare no affiliation or involvement in any organization or entity with financial or non-financial interests in the subject matter or materials discussed in this manuscript.

## Appendix A

### Appendix A.1. Deterministic Homogeneous ODE Epidemic Models, cf. [65]

To put our work into perspective, we recall a definition of deterministic ODE epidemic models lifted from [65].

**Definition A1.** A deterministic epidemic model is a dynamic system with two types of variables $\vec{x}(t) := (\vec{i}(t), \vec{z}(t)) \in \mathbb{R}_+^N$, where

1. $\vec{i}(t)$ model the number (or density) in different compartments of infected individuals (i.e., latent, infectious, hospitalized, etc.), which should (ideally) eventually disappear if the epidemic ever ends;

2. $\vec{z}(t)$ model numbers (or densities) in compartments of individuals who are not infected (i.e. susceptibilities, individuals who are immune, recovered individuals, etc.).
The system must admit an equilibrium called disease-free equilibrium (DFE) and, hence, a "quasi-triangular" linearization in the form

$$ \vec{i}'(t) = \vec{i}(t)A_{i,i}(\vec{x}(t)), \vec{x}(t) \in \mathcal{D} \subset \mathbb{R}_+^N \quad (A1) $$

$$ \vec{z}'(t) = \vec{i}(t)A_{z,i}(\vec{x}(t)) + (\vec{z}(t) - \vec{z}_{dfe})A_{z,z}(\vec{x}(t)), $$

where $\mathcal{D}$ is some forward-invariant subset, where "quasi-triangular" refers to the fact that the functions $A_{i,i}$, $A_{z,i}$, $A_{z,z}$ depend on all variables $\vec{x}(t)$, and where $N$ is the dimension of $\vec{x}(t)$.

As shown in [65], any epidemic model admitting an equilibrium point $(\vec{0}, \vec{z}_{dfe})$ admits the representation (A1), under suitable smoothness assumptions. In what follows, we will call the point $x_{dfe} = (\vec{0}, \vec{z}_{dfe})$ a disease-free equilibrium (DFE).

**Remark A1.** Note that the essential feature of (A1) is the "factorization of the disease equations".

## Appendix A.2. The Basic Reproduction Number $\mathcal{R}_0$

The basic reproduction number or "net reproduction rate" $\mathcal{R}_0$ is a pillar concept in demography, population dynamics, branching processes, and mathematical epidemiology—see the introduction of the book [64]. One central objective in these fields involves the stability of DFE, i.e., the conditions that ensure eradicating the sickness (or a part of the population in the population dynamics). In simple models, it was discovered that this amounts to verifying that a famous threshold parameter (called the basic reproduction number) is less than 1.

1. The notation $\mathcal{R}_0$ was first introduced by the mathematical demography Lotka [63,66]. In epidemiology, the basic reproduction number models the expected number of secondary cases, where one infected case would produce a homogeneous, completely susceptible stochastic population in the next generation. In the simplest setup of the branching process, this parameter—being smaller than 1—makes extinction sure. The relation to epidemiology is that epidemics are well approximated by the branching process at inception, a fact that goes back to Bartlett and Kendall.

2. With more infectious classes, one deals at inception with multi-class branching processes; stability holds when the Perron-Frobenius eigenvalue of the "next-generation matrix" (NGM) means it is smaller than 1.

3. For deterministic epidemic models, it seems at first that the basic reproduction number $\mathcal{R}_0$ is lost due to the generation disappearing in this setup; see Chapter 3 in reference [64], who recalled a method that introduced generations, going back to Lotka, and is reminiscent of the iterative Lotka-Volterra approach in solving integrodifferential equations. At the end of the tunnel, a unified method for defining $\mathcal{R}_0$ emerged only much later, via the "next generation matrix" approach [28–30,67]. The final result is that the local stability of the disease-free equilibrium holds off the spectral radius of a certain matrix called "next-generation matrix", which depends only on a set of "infectious compartments" $\vec{i}$ (which we aim to reduce to 0) being less than one. This approach works provided that certain assumptions listed below hold (and, thus, $\mathcal{R}_0$ is undefined when these assumptions are not satisfied.).

(C1) The foremost assumption is that the disease-free equilibrium $(\vec{0}, \vec{z}_{dfe})$ is unique and locally asymptotically stable within the disease-free space $\vec{i} = 0$, meaning that all solutions of

$$ \vec{z}'(t) = (\vec{z}(t) - \vec{z}_{dfe})A_{z,z}(\vec{0}, \vec{z}(t)), \quad \vec{z}(0) = \vec{z}_0 $$

must approach the point $z_{dfe}$ when $t \to \infty$.

(C2) Other conditions are related to an "admissible splitting" as a difference between two parts $\mathcal{F}, \mathcal{V}$, called respectively "new infections", and "transitions".
**Definition A2.** A splitting

$$ \vec{i}'(t) = \mathcal{F}(\vec{i}(t), \vec{z}(t)) - \mathcal{V}(\vec{i}(t), \vec{z}(t)) $$

will be called admissible if $\mathcal{F}, \mathcal{V}$ satisfy the following conditions [29,59]:

$$ \left\{ \begin{array}{l} \mathcal{F}(\vec{0}, \vec{z}(t)) = \mathcal{V}(\vec{0}, \vec{z}(t)) = 0, \\ \mathcal{F}(\vec{i}(t), \vec{z}(t)) \ge 0, \quad \forall (\vec{i}(t), \vec{z}(t)), \\ \mathcal{V}_j(\vec{i}(t), \vec{z}(t)) \le 0, \quad \text{when } \vec{i}_j = 0, \\ \sum_{j=1}^n \mathcal{V}_j(\vec{i}(t), \vec{z}(t)) \ge 0, \forall (\vec{i}(t), \vec{z}(t)), \end{array} \right. \qquad (A2) $$

where the subscript j refers to the j'th component.

**Remark A2.** The splitting of the infectious equations has its origins in epidemiology. Mathematically, it is related to the "splitting of Metzler matrices"—see, for example, [54]. Note that the splitting conditions may be satisfied for several or no subset of compartments (see for example the SEIT model, discussed in [29], Chapter 3 in reference [12]). Unfortunately, for deterministic epidemic models, there is no clear-cut definition of $\mathcal{R}_0$ [13,68,69]. (One possible explanation is that several stochastic epidemiological models may correspond in the limit to the same deterministic model).

(C3) We turn now to the last conditions, which concern the linearization of the infectious equations around the DFE. Using $L = A_{i,i}(\vec{0}, \vec{z}_{dfe})$, and letting $f$ denote the perturbation from the linearization, we may write:

$$ \vec{i}'(t) = \vec{i}(t)L - f(\vec{i}(t), \vec{z}(t)) = \vec{i}(t)(F - V) - f(\vec{i}(t), \vec{z}(t)), \qquad (A3) $$

$$ F := \left[ \frac{\partial \mathcal{F}}{\partial \vec{i}} \right]_{x_{dfe}} V = \left[ \frac{\partial \mathcal{V}}{\partial \vec{i}} \right]_{x_{dfe}}, \quad L = F - V. $$

The "transmission and transition" linearization matrices $F, V$ must satisfy $F \ge 0$ component-wise and $V$ is a non-singular M-matrix, which ensures that $V^{-1} \ge 0$. (The assumption (B) implies that $L = F - V$ is a "stability (non-singular) M-matrix", which is necessary for the non-negativity and boundedness of the solutions Theorems 1–3 in reference [33].)

Under conditions (C1)–(C3), the next-generation matrix method gives an explicit expression for the basic reproduction number, given by $\mathcal{R}_0 := \lambda_{PF}(FV^{-1})$.

The basic reproduction number is a threshold parameter in the following sense ([29], [Thm 1]):

1(a) When $\mathcal{R}_0 < 1$, the DFE is locally asymptotically stable;

1(b) When $\mathcal{R}_0 > 1$, the DFE is unstable;

2 The DFE is globally asymptotically stable when $\mathcal{R} \le 1$, provided that the "perturbation from linearity" $f = i(F - V) - \mathcal{F} + \mathcal{V}$ is non-negative Theorem 2 in reference [29]. (Note that the [70] strong $\mathcal{R}_0$ alternative was only established for a general n-stage progression, which is a particular case of the model we study below, in which A is an "Erlang" upper-diagonal matrix. It is natural to expect that the result continues to hold for other non-singular Metzler matrices).

In what follows, we will call the alternative 1(a)–1(b) the "weak $\mathcal{R}_0$ alternative". In contrast, result 2 is called the "strong $\mathcal{R}_0$ alternative" in [59,70].
## References

1. Arino, J.; Brauer, F.; van den Driessche, P.; Watmough, J.; Wu, J. A final size relation for epidemic models. *Math. Biosci. Eng.* **2007**, *4*, 159. [PubMed]

2. Avram, F.; Adenane, R.; Halanay, A.; Basnarkov, L.; Bianchin, G.; Goreac, D. On matrix-SIR arino models with linear birth rate, waning immunity, disease and vaccination fatalities, and their approximations. *arXiv* **2021**, arXiv:2112.03436.

3. Avram, F.; Adenane, R.; Bianchin, G.; Halanay, A. Stability analysis of an Eight parameter SIR-type model including loss of immunity, and disease and vaccination fatalities. *Mathematics* **2022**, *10*, 402. [CrossRef]

4. Avram, F.; Adenane, R.; Halanay, A. New results and open questions for SIR-PH epidemic models with linear birth rate, waning immunity, vaccination, and disease and vaccination fatalities. *Symmetry* **2022**, *14*, 995. [CrossRef]

5. Brauer, F. The kermack-mckendrick epidemic model revisited. *Math. Biosci.* **2005**, *198*, 119-131. [CrossRef] [PubMed]

6. Breda, D.; Diekmann, O.; De Graaf, W.F.; Pugliese, A.; Vermiglio, R. On the formulation of epidemic models (an appraisal of kermack and mckendrick). *J. Biol. Dyn.* **2012**, *6* (Suppl. 2), 103–117. [CrossRef]

7. Diekmann, O.; Gyllenberg, M.; Metz, J.A.J. Finite dimensional state representation of linear and nonlinear delay systems. *J. Dyn. Differ. Eq.* **2018**, *30*, 1439–1467. [CrossRef]

8. Champredon, D.; Dushoff, J.; Earn, D.J.D. Equivalence of the erlang-distributed seir epidemic model and the renewal equation. *SIAM J. Appl. Math.* **2018**, *78*, 3258–3278. [CrossRef]

9. Diekmann, O.; Inaba, H. A systematic procedure for incorporating separable static heterogeneity into compartmental epidemic models. *arXiv* **2022**, arXiv:2207.02339.

10. Kermack, W.O.; McKendrick, A.G. A contribution to the mathematical theory of epidemics. *Proc. R. Soc. Lond. Ser. A Contain. Pap. A Math. Phys. Character* **1927**, *115*, 700–721.

11. Earn, D.J.D. A light introduction to modelling recurrent epidemics. In *Mathematical Epidemiology*, Springer: Berlin/Heidelberg, Germany, 2008; pp. 3–17.

12. Martcheva, M. *An Introduction to Mathematical Epidemiology*; Springer: New York, NY, USA, 2015; Volume 61.

13. Thieme, H.R. *Mathematics in Population Biology*; Princeton University Press: Princeton, NJ, USA, 2018.

14. Brauer, F.; Castillo-Chavez, C.; Feng, Z. *Mathematical Models in Epidemiology*; Springer: New York, NY, USA, 2019.

15. Schaback, R. On COVID-19 modelling. *Jahresbericht Deutschen Mathematiker-Vereinigung* **2020**, *122*, 167–205. [CrossRef]

16. Bacaër, N. Un modèle mathématique des débuts de l'épidémie de coronavirus en france. *Math. Model. Nat. Phenom.* **2020**, *15*, 29. [CrossRef]

17. Ketcheson, D.I. Optimal control of an SIR epidemic through finite-time non-pharmaceutical intervention. *arXiv* **2020**, arXiv:2004.08848.

18. Charpentier, A.; Elie, R.; Laurière, M.; Tran, V.C. COVID-19 pandemic control: Balancing detection policy and lockdown intervention under ICU sustainability. *Math. Model. Nat. Phenom.* **2020**, *15*, 57. [CrossRef]

19. Djidjou-Demasse, R.; Michalakis, Y.; Choisy, M.; Sofonea, M.T.; Alizon, S. Optimal COVID-19 epidemic control until vaccine deployment. *medRxiv* **2020**. [CrossRef]

20. Sofonea, M.T.; Reyné, B.; Elie, B.; Djidjou-Demasse, R.; Selinger, C.; Michalakis, Y.; Alizon, S. Epidemiological monitoring and control perspectives: Application of a parsimonious modelling framework to the COVID-19 dynamics in France. *medRxiv* **2020**. [CrossRef]

21. Alvarez, E.F.; Argente, D.; Lippi, F. *A Simple Planning Problem for COVID-19 Lockdown*; Technical Report; National Bureau of Economic Research: Cambridge, MA, USA, 2020.

22. Horstmeyer, L.; Kuehn, C.; Thurner, S. Balancing quarantine and self-distancing measures in adaptive epidemic networks. *arXiv* **2020**, arXiv:2010.10516.

23. Di Lauro, F.; Kiss, I.Z.; Miller, J. Optimal timing of one-shot interventions for epidemic control. *PLoS Comput Biol.* **2021**, *17*, e1008763. [CrossRef]

24. Franco, E. A feedback SIR (fSIR) model highlights advantages and limitations of infection-based social distancing. *arXiv* **2020**, arXiv:2004.13216.

25. Baker, R. Reactive social distancing in a SIR model of epidemics such as COVID-19. *arXiv* **2020**, arXiv:2003.08285.

26. Caulkins, J.; Grass, D.; Feichtinger, G.; Hartl, R.; Kort, P.M.; Prskawetz, A.; Seidl, A.; Wrzaczek, S. How long should the COVID-19 lockdown continue? *PLoS ONE* **2020**, *15*, e0243413. [CrossRef] [PubMed]

27. Caulkins, J.P.; Grass, D.; Feichtinger, G.; Hartl, R.F.; Kort, P.M.; Prskawetz, A.; Seidl, A.; Wrzaczek, S. The optimal lockdown intensity for COVID-19. *J. Math. Econ.* **2021**, *93*, 102489. [CrossRef] [PubMed]

28. van den Driessche, P.; Watmough, J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. *Math. Biosci.* **2002**, *180*, 29–48. [CrossRef]

29. Van den Driessche, P.; Watmough, J. Further notes on the basic reproduction number. In *Mathematical Epidemiology*, Springer: Berlin/Heidelberg, Germany, 2008; pp. 159–178.

30. Diekmann, O.; Heesterbeek, J.A.P.; Metz, J.A.J. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. *J. Math. Biol.* **1990**, *28*, 365–382. [CrossRef]

31. Arino, J.; Portet, S. A simple model for COVID-19. *Infect. Dis. Model.* **2020**, *5*, 309–315. [CrossRef]

32. Avram, F.; Adenane, R.; Ketcheson, D. A review of matrix SIR arino epidemic models. *Mathematics* **2021**, *9*, 1513. [CrossRef]
33. De la Sen, M.; Nistal, R.; Alonso-Quesada, S.; Ibeas, A. Some formal results on positivity, stability, and endemic steady-state attainability based on linear algebraic tools for a class of epidemic models with eventual incommensurate delays. *Discret. Dyn. Nat. Soc.* **2019**, *2019*, 8959681. [CrossRef]

34. Busenberg, S.; van den Driessche, P. Analysis of a disease transmission model in a population with varying size. *J. Math. Biol.* **1990**, *28*, 257–270. [CrossRef]

35. Busenberg, S.; van den Driessche, P. A method for proving the non-existence of limit cycles. *J. Math. Anal. Appl.* **1993**, *172*, 463–479. [CrossRef]

36. Derrick, W.R.; van den Driessche, P. A disease transmission model in a nonconstant population. *J. Math. Biol.* **1993**, *31*, 495–512. [CrossRef]

37. Greenhalgh, D. Hopf bifurcation in epidemic models with a latent period and nonpermanent immunity. *Math. Comp. Model.* **1997**, *25*, 85–107. [CrossRef]

38. Li, M.Y.; Graef, J.R.; Wang, L.; Karsai, J. Global dynamics of a SEIR model with varying total population size. *Math. Biosci.* **1999**, *160*, 191–213. [CrossRef] [PubMed]

39. Sun, C.; Hsieh, Y.-H. Global analysis of an SEIR model with varying population size and vaccination. *Appl. Math. Model.* **2010**, *34*, 2685–2697. [CrossRef]

40. Razvan, M.R. Multiple equilibria for an SIRS epidemiological system. *arXiv 2001*, arXiv:0101051.

41. Li, J.; Ma, Z. Qualitative analyses of SIS epidemic model with vaccination and varying total population size. *Math. Comp. Model.* **2002**, *35*, 1235–1243. [CrossRef]

42. Yang, W.; Sun, C.; Arino, J. Global analysis for a general epidemiological model with vaccination and varying population. *J. Math. Anal. Appl.* **2010**, *372*, 208–223. [CrossRef]

43. Ma, J.; Earn, D.J.D. Generality of the final size formula for an epidemic of a newly invading infectious disease. *Bull. Math. Biol.* **2006**, *68*, 679–702. [CrossRef]

44. Feng, Z. Final and peak epidemic sizes for SEIR models with quarantine and isolation. *Math. Biosci.* **Eng. 2007**, *4*, 675.

45. Andreasen, V. The final size of an epidemic and its relation to the basic reproduction number. *Bull. Math. Biol.* **2011**, *73*, 2305–2321. [CrossRef]

46. Liu, W.-M.; Levin, S.A.; Iwasa, Y. Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models. *J. Math. Biol.* **1986**, *23*, 187–204. [CrossRef]

47. Liu, W.-M.; Hethcote, H.W.; Levin, S.A. Dynamical behavior of epidemiological models with nonlinear incidence rates. *J. Math. Biol.* **1987**, *25*, 359–380. [CrossRef] [PubMed]

48. Georgescu, P.; Hsieh, J.-H. Global stability for a virus dynamics model with nonlinear incidence of infection and removal. *SIAM J. Appl. Math.* **2007**, *67*, 337–353. [CrossRef]

49. Tang, Y.; Huang, D.; Ruan, S.; Zhang, W. Coexistence of limit cycles and homoclinic loops in a SIRS model with a nonlinear incidence rate. *SIAM J. Appl. Math.* **2008**, *69*, 621–639. [CrossRef]

50. Riaño, G. Epidemic models with random infectious period. *medRxiv 2020*. [CrossRef]

51. Plemmons, R.J. M-matrix characterizations. I—Nonsingular M-matrices. *Linear Algebra Its Appl.* **1977**, *18*, 175–188. [CrossRef]

52. Hurtado, P.J.; Krosingh, A.S. Generalizations of the 'linear chain trick': Incorporating more flexible dwell time distributions into mean field ODE models. *J. Math. Biol.* **2019**, *79*, 1831–1883. [CrossRef]

53. Hyman, J.M.; Li, J.; Stanley, E.A. The differential infectivity and staged progression models for the transmission of hiv. *Math. Biosci.* **1999**, *155*, 77–109. [CrossRef]

54. Fall, A.; Iggidr, A.; Sallet, G.; Tewa, J.-J. Epidemiological models and Lyapunov functions. *Math. Model. Nat. Phenom.* **2007**, *2*, 62–83. [CrossRef]

55. Beretta, E.; Takeuchi, Y. Global stability of an sir epidemic model with time delays. *J. Math. Biol.* **1995**, *33*, 250–260. [CrossRef]

56. Andò, A.; Breda, D.; Gava, G. How fast is the linear chain trick? A rigorous analysis in the context of behavioral epidemiology. *Math. Biosci.* **Eng. 2020**, *17*, 5059–5085. [CrossRef]

57. Champredon, D.; Dushoff, J. Intrinsic and realized generation intervals in infectious-disease transmission. *Proc. R. Soc. B Biol. Sci.* **2015**, *282*, 20152026. [CrossRef] [PubMed]

58. Horn, A.R.; Johnson, C.R. *Matrix Analysis*; Cambridge University Press: New York, NY, USA, 2012.

59. Shuai, Z.; van den Driessche, P. Global stability of infectious disease models using Lyapunov functions. *SIAM J. Appl. Math.* **2013**, *73*, 1513–1532. [CrossRef]

60. De la Sen, M.; Ibeas, A. On an SE(Is)(Ih)AR epidemic model with combined vaccination and antiviral controls for COVID-19 pandemic. *Adv. Diff. Eq.* **2021**, *2021*, 1–30. [CrossRef]

61. Shah, K.; Din, R.U.; Deebani, W.; Kumam, P.; Shah, Z. On nonlinear classical and fractional order dynamical system addressing COVID-19. *Results Phys.* **2021**, *24*, 104069. [CrossRef] [PubMed]

62. Tang, T.-Q.; Jan, R.; Bonyah, E.; Shah, Z.; Alzahrani, E. Qualitative analysis of the transmission dynamics of dengue with the effect of memory, reinfection, and vaccination. *Comput. Math. Methods Med.* **2022**. [CrossRef]

63. Lotka, A.J. *Analyse Démographique Avec Application Particulière à L'espèce Humaine*; Hermann: Paris, France, 1939.

64. Bacaër, N. *Mathématiques et Epidémies*; HAL: Lyon, France, 2021.

65. Claude Kamgang, J.; Sallet, G. Computation of threshold conditions for epidemiological models and global stability of the disease-free equilibrium (DFE). *Math. Biosci.* **2008**, *213*, 1–12. [CrossRef]
66. Dietz, K. The estimation of the basic reproduction number for infectious diseases. *Stat. Methods Med. Res.* **1993**, *2*, 23–41.
[CrossRef]

67. Perasso, A. An introduction to the basic reproduction number in mathematical epidemiology. *ESAIM Proc. Surv.* **2018**, *62*, 123–138. [CrossRef]

68. Roberts, M.G. The pluses and minuses of O. *J. R. Soc. Interface* **2007**, *4*, 949–961. [CrossRef]

69. Li, J.; Blakeley, D. The failure of R0. *Comput. Math. Methods Med.* **2011**, *2011*, 527610. [CrossRef]

70. Guo, H.; Li, M.Y. Global dynamics of a staged progression model for infectious diseases. *Math. Biosci. Eng.* **2006**, *3*, 513. [CrossRef] [PubMed]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.