We revisit here a landmark five-parameter SIR-type model, which is maybe the simplest example where a complete picture of all cases, including non-trivial bistability behavior, may be obtained using simple tools. We also generalize it by adding essential vaccination and vaccination-induced death parameters, with the aim of revealing the role of vaccination and its possible failure. The main result is Theorem 1, which describes the stability behavior of our model in all possible cases.

Many of the models used nowadays in mathematical epidemiology, in particular in COVID-19 research, belong to a certain subclass of compartmental models whose classes may be divided into three “$(x,y,z)$” groups, which we will call respectively “susceptible/entrance, diseased, and output” (in the classic SIR case, there is only one class of each type). Roughly, the ODE dynamics of these models contains only linear terms, with the exception of products between x and y terms. It has long been noticed that the reproduction number $\mathcal{R}$ has a very simple Formula in terms of the matrices which define the model, and an explicit first integral Formula is also available. These results can be traced back at least to Arino, Brauer, van den Driessche, Watmough, and Wu (2007) and to Feng (2007), respectively, and may be viewed as the “basic laws of SIR-type epidemics”. However, many papers continue to reprove them in particular instances. This motivated us to redraw attention to these basic laws and provide a self-contained reference of related formulas for $(x,y,z)$ models. For the case of one susceptible class, we propose to use the name SIR-PH, due to a simple probabilistic interpretation as SIR models where the exponential infection time has been replaced by a PH-type distribution. Note that to each SIR-PH model, one may associate a scalar quantity $Y(t)$ which satisfies “classic SIR relations”,which may be useful to obtain approximate control policies.

The basic reproduction number $R_0$ is a concept which originated in population dynamics, mathematical epidemiology, and ecology and is closely related to the mean number of children in branching processes (reflecting the fact that the phenomena of interest are well approximated via branching processes, at their inception). Despite the very extensive literature around $R_0$ for deterministic epidemic models, we believe there are still aspects which are not fully understood. Foremost is the fact that $R_0$ is not a function of the original ODE model, unless we also include in it a certain $(F,V)$ gradient decomposition, which is not unique. This is related to the specification of the “infected compartments”, which is also not unique. A second interesting question is whether the extinction probabilities of the natural continuous time Markovian chain approximation of an ODE model around boundary points (disease-free equilibrium and invasion points) are also related to the $(F,V)$ gradient decomposition. We offer below several new contributions to the literature: (1) A universal algorithmic definition of a $(F,V)$ gradient decomposition (and hence of the resulting $R_0$). (2) A fixed point equation for the extinction probabilities of a stochastic model associated to a deterministic ODE model, which may be expressed in terms of the $(F,V)$ decomposition. Last but not least, we offer Mathematica scripts and implement them for a large variety of examples, which illustrate that our recipe offers always reasonable results, but that sometimes other reasonable $(F,V)$ decompositions are available as well.

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 $(\overrightarrow{\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 $\overrightarrow{\alpha }$ is to restrict it to the case where there is only one susceptible class $\mathsf{s}$, and when matrix B (given by the part of the new infection matrix, with respect to $\mathsf{s}$) is of rank one, with $B=b\overrightarrow{\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.).

Our paper reviews some key concepts in chemical reaction network theory and mathematical epidemiology, and examines their intersection, with three goals. The first is to make the case that mathematical epidemiology (ME), and also related sciences like population dynamics, virology, ecology, etc., could benefit by adopting the universal language of essentially non-negative kinetic systems as developed by chemical reaction network (CRN) researchers. In this direction, our investigation of the relations between CRN and ME lead us to propose for the first time a definition of ME models, stated in Open Problem 1. Our second goal is to inform researchers outside ME of the convenient next generation matrix (NGM) approach for studying the stability of boundary points, which do not seem sufficiently well known. Last but not least, we want to help students and researchers who know nothing about either ME or CRN to learn them quickly, by offering them a Mathematica package “bootcamp”, including illustrating notebooks (and certain sections below will contain associated suggested notebooks; however, readers with experience may safely skip the bootcamp). We hope that the files indicated in the titles of various sections will be helpful, though of course improvement is always possible, and we ask the help of the readers for that.

Our paper presents three new classes of models: SIR-PH, SIR-PH-FA, and SIR-PH-IA, and states two problems we would like to solve about them. Recall that deterministic mathematical epidemiology has one basic general law, the “${\mathcal{R}}_0$ alternative” of Van den Driessche and Watmough, which states that the local stability condition of the disease-free equilibrium may be expressed as ${\mathcal{R}}_0<1$, where ${\mathcal{R}}_0$ is the famous basic reproduction number, which also plays a major role in the theory of branching processes. The literature suggests that it is impossible to find general laws concerning the endemic points. However, it is quite common that 1. When ${\mathcal{R}}_0>1$, there exists a unique fixed endemic point, and 2. the endemic point is locally stable when ${\mathcal{R}}_0>1$. One would like to establish these properties for a large class of realistic epidemic models (and we do not include here epidemics without casualties). We have introduced recently a “simple” but broad class of “SIR-PH models” with varying populations, with the express purpose of establishing for these processes the two properties above. Since that seemed still hard, we have introduced a further class of “SIR-PH-FA” models, which may be interpreted as approximations for the SIR-PH models, and which include simpler models typically studied in the literature (with constant population, without loss of immunity, etc.). For this class, the first “endemic law” above is “almost established”, as explicit formulas for a unique endemic point are available, independently of the number of infectious compartments, and it only remains to check its belonging to the invariant domain. This may yet turn out to be always verified, but we have not been able to establish that. However, the second property, the sufficiency of ${\mathcal{R}}_0>1$ for the local stability of an endemic point, remains open even for SIR-PH-FA models, despite the numerous particular cases in which it was checked to hold (via Routh–Hurwitz time-onerous computations, or Lyapunov functions). The goal of our paper is to draw attention to the two open problems above, for the SIR-PH and SIR-PH-FA, and also for a second, more refined “intermediate approximation” SIR-PH-IA. We illustrate the current status-quo by presenting new results on a generalization of the SAIRS epidemic model.