key: cord-0530768-my628q8h authors: Avram, Florin; Adenane, Rim; Bianchin, Gianluca; Halanay, Andrei title: Stability analysis of an eight parameter SIR-type model including loss of immunity, and disease and vaccination fatalities date: 2021-12-16 journal: nan DOI: nan sha: e08e1144782480adf25644e98a250568a651795f doc_id: 530768 cord_uid: my628q8h We revisit here a landmark five parameter SIR-type model of [DvdD93, Sec. 4], 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 5, which describes the stability behavior of our model in all possible cases. Motivation. This paper has the dual purpose of providing a short guide to students of deterministic mathematical epidemiology, among which we count ourselves, and at the same time to illustrate the technical work one is faced with in an elementary, but not simple "exercise". Of course, one can easily find at least five mustread excellent textbooks and thesis surveying this field (with different emphasis: epidemiological, stability, or control), see for example [AAM92, ST11, Mar15, Thi18, BCCF19, Mon20, Bac21, DM21] , and also at least a hundred major papers which are a must read. We hope however that our little guide may help future students decide as to the order in which these materials must be assimilated. A bit of history. Deterministic mathematical modelling of diseases started with works of Ross on malaria, and imposed itself after the work of [KM27] on the Bombay plague of 1905-06. This was subsequently followed by works on measles, smallpox,chickenpox, mumps, typhoid fever and diphtheria, see for example [Ear08] , and recently the COVID-19 pandemic (see for example [Sch20, Bac20, Ket20, CELT20, DDMC + 20, SRE + 20, AAL20, HKT20, DLKM20, Fra20, Bak20, CGF + 20, CGF + 21], to cite just a few representatives of a huge literature). Note that at its beginning, mathematical epidemiology was a collection of similar examples dealing with current epidemics, and no one bothered to define precisely what is an epidemic, mathematically. One answer to this question, reviewed in Appendix B, may be found in the recent paper [KS08] . The deterministic mathematical epidemiology literature may be divided into three streams. 1. "Constant total population" models are the easiest to study. However, since death is an essential factor of epidemics, the assumption of constant total population (clearly a short term or very large population approximation) deserves some comment. One possible rigorous justification for deterministic constant population epidemiological models comes from slow-fast analysis [KK02, JKKPS21, Gin21] . This is best understood for models with demography (birth, death), which happen typically on a slower scale than the infectious phenomena. Here there is a natural partition of the compartments into a vector i(t) ∈ R − + of disease/infectious compartments (asymptomatic, infectious hospitalized,etc). These interact (quickly) with the other input classes like the susceptibles and output classes like the recovered and dead. For constant total population" models the total population N clearly plays no role, and one may use the assumption that the rate of infection is independent of N , of the form βS(t)I(t) (which was called pseudo, or simple mass-action incidence [dJDH95] ). Note that this simplification of the rate of infection (adopted in the majority of the literature) is inappropriate for varying population models. 2. Models with constant birth rate (in the analog stochastic model, this would correspond to emigration). These models include the previous class, and preserve some of its nice features, like the uniqueness of the endemic fixed point. They typically satisfy the "R 0 alternative", established via the next generation matrix approach, and also the fact that the endemic point exists only when R 0 > 1 -see [DlSNAQI19] for a recent nice review of several stability results for this class of models. 3. Finally, we arrive to the class our paper is concerned with: models with linear birth-rate ΛN ( or constant birth-rate per capita in the analog stochastic model), varying total population N , and "proportionate mixing" rate of infection βS(t) I(t) N (t) . As far as we know, this stream of literature was initiated in [BvdD90, BvdD93] , and allows the possibility of bi-stability when R 0 < 1 (absent from the previous models), even in a SIR example with five parameters [DvdD93] . ¶ . This last stream of literature is quite important, since constant birth-rate per capita is a natural assumption. Despite further remarkable works on particular examples -see for example [Gre97, LGWK99, Raz01, LM02, SH10, YSA10, LL18] (which preferred all direct stability analyses to the next generation matrix approach), the literature on models with varying total population, unlike the two preceding streams, has not yet reached general results. To understand this failure, it seemed to us a good idea to revisit an important SIR model with diseaseinduced deaths and loss of immunity [DvdD93, Raz01] . These important works illustrate already some of the complexities which may arise for varying population models, in particular the possibility of bi-stability when R 0 < 1. This surprising fact lead us to introduce the concept of strong global stability of the DFE (disease free equilibrium), and to find conditions for this to hold in our example, (see Proposition 4). Since the method used is simply linear programming, we hope to extend this in future work. ¶ This reveals that for an initial high number of infectives, the trajectory may lie in the basin of attraction of a stable endemic point instead of being eradicated. The discrepancy with what is expected from the corresponding stochastic model suggests that in this range the deterministic model is inappropriate Our model. We further added to the extension of [KM27] introduced in [DvdD93, Raz01] vaccinations. Grouping together the recovered and vaccinated, our "SIR/V+S" model is described by: It involves six states: S : R ≥0 → N >0 describing the number of susceptible individuals in the population, I : R ≥0 → N >0 describing the number of infections, R : R ≥0 → N >0 describing the number of recovered or vaccinated, D : R ≥0 → N >0 describing the number of natural deaths in the population, D e : R ≥0 → N >0 describing the number of deaths originated by the disease, and N : R ≥0 → N >0 describing the total number of (alive) individuals in the population. The parameters Λ ∈ R ≥0 and µ ∈ R ≥0 denote the average birth and death rates in the population (in the absence of the disease), respectively, γ s ∈ R ≥0 is the vaccination rate, γ r ∈ R ≥0 denotes the rate at which immune individuals lose immunity (this is the reciprocal of the expected duration of immunity), γ ∈ R >0 is the rate at which infected individuals recover from the disease, ν i ∈ R ≥0 is the extra death rate due to the disease, and ν r ∈ R ≥0 is an average death rate in the recovered-vaccinated compartment (due to e.g. deaths caused by the vaccine). Note that in what follows we use the notation γ c to denote the total rate at which individuals leave a certain compartment C towards other non-deceased compartments, and we use ν c to denote the rate at which individuals leave compartment C towards a deceased compartment. In (1.1), susceptible individuals become infected at rate β N (t) I(t) (thus moving to the I compartment), they are vaccinated at rate γ s (thus moving to the R compartment), and deaths occur at rate µ (thus moving to the D compartment); infected individuals recover at rate γ (thus moving to the R compartment), die of nondisease related causes at rate µ (thus moving to the D compartment), and die of disease-related causes at rate ν i (thus moving to the D e compartment); recovered individuals lose immunity at rate γ r (thus moving to the S compartment), die of non-disease related causes at rate µ (thus moving to the D compartment), die of diseaserelated causes at rate ν r , and become re-infected at rate βr N (t) I(t) (thus moving to the I compartment) (see Remark 1 for a discussion on re-infections). Note that D, D e are completely determined once the other classes are found. These "output classes" will not be mentioned further (since they are only relevant in control problems, which are outside our scope). The dynamics of I, the disease class, allows computing the basic reproduction number via the next generation matrix method. Finally, the input classes S, R determine by themselves the disease free equilibrium. Remark 1 Models of the form (1.1) that account for non-constant population sizes can be especially useful in two scenarios: (i) to study diseases that remain infectious for long periods of time with small disease mortality rate, where the natural death/birth rate of the population plays a central role (such as HIV/AIDS, malaria and tuberculosis), as well as (ii) to study diseases with short infectious periods but with a substantial disease mortality rate, where the death rate due to the disease plays a central role (such as measles, influenza, SARS/COVID). The two-way transfers between the recovered and infected compartments (recall that R βr γ I). can be used to account for multiple variants of the disease, whereby immunity to one variant does not guarantee immunity to all other variants. For instance, diseases such as HIV can develop resistance to medications, and such resistance can be transmitted to a partner. In this cases, even when the second party has recovered, it may become re-infected with a different variant. Notice that, an even more general case is considered in [Raz01] , where vaccinated individuals may transition to the infected compartment. Remark 2 In practical situations, certain parameter relations, for instance βr > β, νr > ν i , might seem "unreasonable from a medical point of view". In what follows, we choose not to assume any relationship among the parameters in (1.1) in order to highlight the fact surprising mathematical behaviors like bistability-see Figure 6d -may arise when "things go wrong". We make now several remarks, which serve as an appetizer for the rest of the paper. Remark 3 The critical value β = βr defines a model where both the sickness and the vaccination do not affect the infectivity (neither confers any immunity). The recovered class might be better viewed then as a susceptible class So of older individuals, with extra mortality rate νr > 0 and no births: A moment of reflection reveals that this particular case has two surprising features: a) after the infection is over, transfers only occur towards the old class So, and b) transfers between the two age groups occur; it is hard to make sense of this, without further imposing γs = γr = 0. In what follows, we will investigate the stability properties of an equivalent system to (1.1) obtained by investigating the normalized quantities s = S N , i = I N , r = R N and given by (see Proposition 1, Section 2 for a detailed derivation). Remark 5 Note this reduces to the classic SIR model [KM27] when Λ = ν i = νr = γr = γs = βr = 0. Remark 6 Factoring the second equation reveals the threshold after which the infection starts decreasing β s(t) + ν i i(t) + (νr + βr) r(t) < Λ + γ + ν i (1.2) herd immunity ( s, i) line, or max-line (line since r = 1 − s − i). Notice that, in contrast with the case of models with constant population size and no loss of immunity where the herd immunity condition depends only on the susceptible state, the above condition depends on three states ( s, i, r). Also, when ν i = νr = βr = 0, this reduces to the well-studied herd immunity threshold. 1) The inequality obtained at the DFE, when i = 0, βs df e + (νr + βr)r df e ≤ Λ + γ + ν i ⇔ βs df e + (νr + βr)r df e Λ + γ + ν i := R 0 < 1 (1.3) turns out to ensure the local stability of the disease free equilibrium (disease free equilibrium) -see section 3. 2) R 0 introduced above coincides with the famous basic reproduction number defined via the next generation matrix approach. 3) The disease free equilibrium (obtained by plugging i = 0, r = 1 − s in the fixed point equation) is such that Remark 7 A) The equality R 0 = 1 is linear in γs and yields the so called "critical vaccination" provided that the denominator doesn't blow-up. This formula is positive if either βr ≤ Λ + γ + ν i ≤ β, or β ≤ Λ + γ + ν i ≤ βr. When βr = 0, we recover a classical critical vaccination formula (1.6) B) The equality R 0 = 1 is also linear in β and yields a "critical contact rate " β * = (γ + ν i + Λ) (Λ + γr + γs) − βrγs Λ + γr . (1.7) It may be checked that at this critical value the value iee of the infectious at the lower endemic point crosses the i = 0 axis, and may reduce the number of endemic points from 3 to 2 -see Figure 7 for details. Contents. Section 2 reviews the dimension reduction available for the proportions of models with linear birth rate, and emphasizes the fact that the well studied deterministic model is an approximation of the model with linear birth rate studied here. A second, finer "intermediate approximation" is introduced as well. Section 3 computes R 0 via the next generation matrix approach (NGM), establishing thereby the wellknown weak R 0 alternative [VdDW08] ; it also introduces the concept of "strong global stability" of the DFE in Proposition 4, which may be useful for more advanced models. The endemic equilibria are discussed in Section 4. Section 5 identifies more precisely, in the particular case ν r = 0, the case when global stability does not hold. The results involve heavily the vaccination parameter γ s and its critical value. The increased complexity of the model forces a geometric approach, already hinted at in [DvdD93, Sec. 4 ]. This ends up in the consideration of 10 cases, one of which, Theorem 52.(c), remains only partly resolved. Section 6 discusses the particular case γ s = 0, generalizing and providing missing details of the results in [DvdD93, Sec. 4 ]. Section 7 gives the simpler results for the first approximation FA (actually for a slightly more general "classic/pedagogical model"). Finally, Section 8 provides the proof of a technical result, and Section 9 reviews the pillar of deterministic epidemic models: the definition of the basic reproduction number via the next generation matrix method. It is convenient to reformulate (1.1) in terms of the normalized fractions This process, sometimes called "non-dimensionalizing" (see for example [Ras21] ), allows us to provide the following equivalent representation of (1.1). Proposition 1 Let s, i, r be as defined in (2.1). Then, the dynamics (1.1) can be equivalently rewritten as: Proof By using 3) we obtain for the susceptibles: Similarly, and Finally, s(t)+ i(t)+ r(t) = 1 follows from N (t) = (S(t)+I(t)+R(t)) by substituting (2.1), which proves the equivalence between (1.1) and (2.2). Remark 8 Note that the natural death rate µ does not intervene in (2.2) , which is to be expected. Indeed, since this rate is the same for all the compartments, it has no effect on the fractions. Remark 9 1. Note that the conservation equation n := s + i + r = 1, in general, does not follow from the first three equations in (2.2). Indeed, by summing the right hand-sides we have: which shows that n (t) = 0 in general cases. However, the above differential equation guarantees that if is forward-invariant, since the flow along its boundaries is directed towards the interior -see [Raz01] § . Note that the conservation equation may replace either the last or the first equation in the dynamics, and allows reducing the computation of fixed points to dimension 2. 2. The Jacobian matrix of (2.2) is given by Following [LGWK99] , let us investigate when this system is order preserving in the interior of the convex feasible region D. The Jacobian 2) is order preserving with respect to the partial ordering defined by the orthant (s, i) ∈ R 2 , s ≥ 0, i ≥ 0 . 2) is order preserving with respect to the partial ordering defined by the orthant (s, i) ∈ R 2 , s ≤ 0, i ≥ 0 . The study of the dynamics in Proposition 1 is quite challenging, and it may be useful sometimes to consider also the two approximations introduced in the following definition. 1. The model (2.5) with Φ s = Φ i = Φ r = 1 will be called scaled model (SM). 2. The model (2.5) with Φ s = Φ i = Φ r = 0 will be called first approximation (FA) in the specific case Λ = µ. The FA is thus: 3. The model (2.5) with Φ s = Φ r = 1 and Φ i = 0 will be called intermediate approximation (IA) . Remark 10 Each model, in particular SIR/V+S, has a SM, FA and IA version, which will be denoted by SIR/V+S-SM, SIR/V+S-FA, SIR/V+S-IA. The FA is not a constant population model when νc > 0, for some compartment c. if and only if µ = Λ and ν i = νr = 0. Thus, the popular assumption of constant population size [Het76] applies only to epidemics without extradeaths, which contradicts the essence of most epidemics [BCCF19, Ch. 10.2]. Clearly, constant population papers have in mind some large N or short term approximation, but this is rather vague. On the other hand, the FA approximation (2.6) with µ = Λ, as well as IA, may be heuristically justified as approximations obtained by ignoring certain quadratic terms in (2.2). This justifies studying FA without restrictive assumptions like ν i = 0. Remark 13 A considerable part of the epidemics literature studies (1.1) with N (t) = 1 (this produces an analog of (2.6) with µ = Λ). The justification for studying this model is of course an assumption that N is "approximately constant". The purpose of our paper is not to assume that; however, we chose to include results about them, under the name "classic/pedagogic models" (PM), to be in line with this part of the literature. As explained, we need here only results on the FA model (which approximate the object of interest to us (2.2)), and these may be easily recovered by replacing µ with Λ. We conclude this section by illustrating in Figure 3 a comparison between the trajectories of the first approximation, of the intermediate approximation, and of the scaled model. Note that the approximate dynamics are an accurate approximation of the SM at the beginning of the epidemic (i.e., when i(t) ≈ 0). This period starts with the lower part in Figure 3 , and continues until the processes start turning towards their distinct endemic points -see [JKKPS21] Note that the epidemic will spend at first a long time (since births and deaths have slow rates as compared to the disease) in the vicinity of the manifold i(t) = 0, where the three processes are indistinguishable, before turning towards the endemic equilibrium point(s). It is convenient to eliminate r from s + i + r = 1, and work with the following two-dimensional scaled dynamic defined on the positively invariant region [Raz01] The fixed points are the solutions of The disease free system ( with i = 0, r = 1 − s) reduces to and its fixed points are such that s satisfies the equation One root is always in [0, 1] and will be denoted by s df e . Remark 14 s df e is continuous in νr, since for νr small, Remark 15 The other root in the quadratic case νr > 0 in which case it yields a second DFE point with s = 0 = i. Assumption 1 We will exclude from now on the particular boundary case (2.12), which may be resolved by elementary explicit eigenvalues computations -see [Raz01] § . Once this case is removed, the DFE is unique, and we may apply the next generation matrix method, the first step of which consists in checking the local asymptotic stability of s df e for the disease-free equation (2.9). In our case, this amounts to proving that This is automatic both when ν r = 0, and when ν r > 0, since by (2.10) λ − = − ∆ df e < 0. 3 Local DFE stability via the next generation matrix approach [vdDW02] Recall (cf. Assumption 1) that we exclude the case Λ = γ r = 0, ν r ≥ γ s (2.12), so that the DFE defined in (2.10) is unique and locally stable in the disease free space (2.13), so that the next generation matrix approach may be attempted. This proceeds in the following steps: 1. Separate the "infected equation" for i (2.2) as a difference of two terms, F, V, as follows . The necessary decomposition conditions on F, V (9.2) are immediate, and the gradients F = ∂F ∂ i i=0,s=s df e = βs df e + (ν r + β r )r df e , V = ∂V ∂ i i=0,s=s df e = γ + ν i + Λ, § Note however that while not necessarily interesting from an epidemics point of view, this case is remarkable however mathematically.(More precisely, the extra DFE point (0, 0) may be either source or saddle-point, and there are two endemic points, which may be either a sink and a saddle, or two sinks [Raz01] ; finally, for the general SIR/V+S model, both cases may be achieved as small perturbations of the particular case (2.12)) satisfyF ≥ 0, V > 0, as required. One may conclude then, by the NGM result [VdDW08, Thm 1], that: The DFE is locally stable if the Perron-Frobenius eigenvalue of the next generation matrix satisfies and is unstable if R 0 > 1. The equality R 0 = 1 is linear in γs and yields provided that the denominator doesn't blow-up. The parameter γ * s is called "critical vaccination". The daunting formula (3.2) simplifies and turns out to provide crucial help for stability analysis in the following particular cases ‡ νr = 0 =⇒ (3. 3) The first formula is positive if either βr ≤ Λ + γ + ν i ≤ β, or β ≤ Λ + γ + ν i ≤ βr, and a geometric consequence of this is provided in Section 5. Remark 17 (Direct stability analysis for the DFE) The Jacobian of (2.7) (with r eliminated) is We recover the same result as via the next generation matrix method. The DFE is locally stable iff λ − < 0, and λ P < 0 ⇔ s df e β + (νr + βr)r df e Λ + γ + ν i < 1. Proposition 3 If R 0 < 1 and DFE is the unique equilibrium point, then it is globally asymptotically stable in the invariant set D. Proof: Each solution starting in D is obviously bounded so its ω-limit set is not empty. The Poincare-Bendixon Theorem implies that this is the unique equilibrium point DFE, since else it would be a closed orbit ( see [HS74] ) and one may check similarly as in [DvdD93, Thm 3.1] that there exist no periodic orbits. Since the explicit conditions for the uniqueness of the equilibrium point are pretty complicated, we prefer to give them only in a particular case, in section 5. For the general case, we may find a simple criteria if we weaken the concept of global stability of the disease free equilibrium as follows then the DFE is "strongly globally stable", in the sense that the function L( s, i) = i is Lyapunov over the invariant In lay terms, we may say that "the epidemics never picks up" in this case. Proof: The non-negative function L( s, i) = i, with L(1, 0) = 0, is Lyapunov (i.e. the epidemics may never increase) iff The maximimization of i reduces thus to a simple linear programming problem. Thus, there must exist an extremal point of D where the maximum is attained, and max (s,i)∈D The endemic equilibrium set may be obtained algebraically by solving one of the variables from i / i = 0 in (2.8), and plugging into the other equation. Eliminating s using yields a quadratic equation . Since locating analytically the endemic points requires considerable effort, we will restrict starting with the next section to the case ν r = 0. We note however already that: 1. The equation is quadratic iff β = β r , β = ν i and β r + ν r = ν i . We will exclude at first these particular cases, but note that they suggest that different regimes occur when the respective thresholds are crossed -this will be further explained below. 2. Locating the endemic points may be achieved geometrically by studying the intersection of the vertical and horizontal isoclines: (a) the i = 0 isocline is given by i = 0 and by the herd immunity line which is a hyperbola. This passes by definition through the DFE, and intersects therefore necessarily the domain when s df e < 1 ⇔ γ s > 0. 5 The case ν r = 0, γ s > 0 Getting sharper stability results beyond the weak R 0 alternative requires locating the endemic points, and this seems quite difficult in general. Therefore, we will restrict from now on to the case ν r = 0, which avoids the necessity of handling the complications arising from the square root formula of s df e . The case γ s = 0, essentially covered in [DvdD93, Sec. 4], requires special treatment -see Section 6. The quadratic equation Ai 2 + Bi + C = 0 has coefficients . (5.1) The endemic points are still complicated: , where ∆ = 4 (β − βr) (β − ν i ) (Λ (−ν i + βr + νr) + (γ + µ)γr) (5.2) By solving C = 0, we may compute however an important critical value for β, above which one (the higher) of the endemic points crosses above the i = 0 axis, entering therefore D (equivalently, R 0 become larger than 1). Locating the endemic points may be attempted algebraically -see Lemma 2 in the Appendix seemed quite difficult. We resorted therefore to a geometric study of the isoclines in Theorem 5, attempting to explain geometrically all the possible cases. For example, the case R 0 > 1 turned out equivalent to the unicity of the endemic point, and to the fact that the immunity line intersection with the domain lie on both sides of the s = 0 isocline. Remark 20 Since the number of crossing points must be odd on one hand, and less than two on the other, identifying such a crossing is equivalent to the uniqueness of the endemic point. We turn now to listing some elementary geometric facts, in particular the coordinates of various intersections points. 1. When ν i = β, the s = 0 isocline becomes the hyperbola with asymptotes s = γr νi−β , i = Λ+γr+γs νi−β . Note that the center is in the first quadrant when β < ν i and in the third quadrant otherwise, and that the intersection E with the line s = 0 has i E = 1 + Λ γr , outside D. Also important will be the convexity of the branch which intersects D. From i (s df e ) = 2 s df e + k (s + k) 3 we find that our branch is convex when β > ν i and concave otherwise. The equality case ν i = β is analyzed in the following remark. which never belongs to the feasible region. Indeed iee > 0 if and only if ν i > βr, and this implies further , 0), β r = β (when β r = β, B goes to ∞). This point satisfies: The six cases listed above are the basis of our geometric analysis provided in Theorem 5. Remark 22 s B coincides with s df e iff γs = γ * s , which fits with the fact that γ * s ∈ (0, 1) iff one of these two cases occurs -recall Remark 16. 3. Another important point is the point A where the immunity line (4.1) intersects i = 1 − s, with coordinates It is easy to show that A ∈ D iff β > Λ + ν i + γ, that it moves then to the fourth quadrant when Λ + ν i + γ ≥ β > ν i , and that it jumps from the fourth to the second quadrant when β decreases below the threshold ν i . 4. The point where the immunity line intersects s = 0 is C(0, 1 + γ+Λ νi−βr ). This is negative when ν i < β r < Λ + γ + ν i , in the domain when β r ≥ Λ + γ + ν i , and bigger than 1 when β r ≤ ν i . More precisely, 5. The unique point D where the hyperbola intersects s + i = 1 within the domain has coordinates . Remark 23 When β = βr, the slope β−βr βr−νi of the immunity line becomes 0 and The dynamical system admits a unique endemic point and a unique EE given by which may be checked to belong to the feasible region iff Indeed, iee > 0 requires either β < ν i , which leads subsequently to a contradiction, or β > Λ + γ + ν i , which may be shown to imply s > 0, s + i < 1. The stability analysis of this case is elementary and left to the reader. The formula (5.4) and the subsequent Remark 22 suggest splitting the analysis according to the order of the three quantities β, Λ + γ + ν i , β r and on the orders of γ s , γ * s = (Λ + γ r ) β−(Λ+γ+νi) Λ+γ+νi−βr and of (ν i , β). We end up with ten cases, nine of which are fully resolved in Theorem 5, and one of which is left partly open. Note that, as proved in the Appendix, these ten cases form a disjoint decomposition of the parameter space, if only strict inequalities are allowed. Before stating Theorem 5, we provide graphical illustrations of the 10 cases. r ] < β r . For β < ν i we are in the case of Fig. 4 (b) (Thm. 5.2.b)), with the immunity line to the left of the hyperbola, and no endemic points. The same situation occurs for ν i < β < β 2 , where β 2 is the largest root of ∆(β) = 0, except that the hyperbola changes to convex -see Fig. 4 (c) (Thm. 5.2.c)(i)). After β 2 two endemic points emerge -see Fig. 4 (d) (Thm. 5.2.c)(ii)). The lower endemic point exits through the boundary i = 0 at β * , defined in (5.3), after which the remaining endemic point remains the only stable point. Theorem 5 Suppose νr = 0, and that neither two of the three parameters β, βr, ν i coincide. Then, one of the following cases must arise: 1. R 0 > 1 ⇔ precisely one endemic point lies in D, which may occur in one of the following four ways: (a) AC ∈ D crosses the hyperbola with positive slope Λ + γ + ν i ≤ β r < β -see Figure 3 (a); (b) AC ∈ D crosses the hyperbola with negative slope Λ + γ + ν i ≤ β < β r -see Figure 3 (b); (c) β r ≤ Λ + γ + ν i ≤ β, γ s < γ * s ⇔ BA ∈ D crosses the hyperbola - Fig. 3(c) ; (d) β ≤ Λ + γ + ν i ≤ β r , γ s ≥ γ * s ⇔ BC ∈ D crosses the hyperbola - Fig. 3(d) . In all these cases, the endemic point is a sink and DFE is a saddle point. 2. R 0 ≤ 1 is the union of the following six cases (similar to the above, but taking also into account the convexity of the hyperbola branch, in some cases): (a) β r < Λ + γ + ν i < β, ν i < β, when both the points A, B lie to the right of a convex hyperbola branch, and the isoclines do not intersect in D-see Figure 4 (a). (b) β < ν i < Λ + γ + ν i < β r , when both the points C, B lie to the left of a concave hyperbola branch, and the isoclines do not intersect in D-see Figure 4 (b). (c) When both the points C, B lie to the left of a convex hyperbola branch, we have two subcases: (i) When ν i < β < Λ + γ + ν i < β r , γ s < γ * s , ∆ < 0, the isoclines intersect D, but do not intersect each other -see Figure 4 (c). (ii) The isoclines intersect in D, yielding two endemic points -see Figure 4 (d). Necessary conditions for this are ν i < β< Λ + γ + ν i < β r , γ s < γ * s , ∆ ≥ 0, and we conjecture that the necessary and sufficient conditions are obtained by adding β r ≥ β (+) r , where β (+) r is defined in (5.8). In this case, the DFE is one of two sink points, whose attraction domains are separated by the separatrices of the third fixed saddle point. (d) The i /i = 0 isocline does not intersect the interior of D in the following two cases: (i) β ≤ β r ≤ Λ + γ + ν i , with hyperbole concave -see figure 5(a)-or convex, according to whether β > ν i or not; (ii) β r < β ≤ Λ + γ + ν i , with hyperbole convex -see figure 5 (b)-or convave, according to whether β < ν i or not; 3. In all the cases when the DFE is the unique fixed point within the domain, it is a globally stable sink. 1. As noted already, β > Λ + γ + ν i is equivalent to fact that the point A = ( γ+Λ β−νi , β−νi−γ−Λ β−νi ) lies in D, which applies in the cases 1.(a), 1.(b), 1.(c) . Furthermore, we may check that in this case the point A is always to the right (above) the hyperbola, i.e. To conclude the existence of a unique endemic point, it is enough then to find in these three cases a point of the immunity line below the hyperbola. Referring to Figures 3, we see that the following cases may arise: (a) AC ∈ D crosses the hyperbola, and s B < 0. We must then be in the case Λ + γ + ν i < β r < β, which implies β r ≥ Λ + ν i + γ, and so C ∈ D is below the hyperbola -see Figure 3 (a). (b) AC ∈ D crosses the hyperbola, and s B ≥ 1. We must then be in the case Λ + γ + ν i < β < β r , which implies again β r ≥ Λ + ν i + γ, and so C ∈ D is again below the hyperbola -see Figure 3 (b). (c) BA ∈ D crosses the hyperbola when R 0 > 1, β > Λ + γ + ν i > β r , s B ≤ s df e ⇔ γ s ≤ γ * s , since B is below the hyperbola -see Figure 3 (c). (d) BC ∈ D crosses the hyperbola (and A / ∈ D) when R 0 > 1, β r > Λ + γ + ν i > β, s B > s df e ⇔ γ s ≥ γ * s , since B is below the hyperbola -see Figure 3 (d). 2. We turn now to the case R 0 ≤ 1, when the point A lies in the fourth quadrant if ν i < β, and in the second quadrant otherwise. (a) In this case s B > s df e , s A > s D . Since the immunity line is to the right of a convex hyperbola branch, it is clear geometrically that they cannot intersect within the domain -see Figure 4 (a). (b) Similarly, the end points in D of the immunity line are to the left of a concave hyperbola branch, and so they cannot intersect within the domain -see Figure 4 (b). (c) For two endemic points to exist, it is necessary that s B = Λ+νi+γ−βr β−βr ∈ (0, s df e ), which requires that and this implies that the other intercept C = (0, Λ+γ+νi−βr νi−βr ) is also in D. The i /i = 0 isocline intersects then D, and may also intersect the hyperbola, when the discriminant satisfies ∆ ≥ 0. The inequality ∆ > 0 is quadratic in β r and maybe rewritten as We conjecture based on numerical evidence that the two endemic points belong to D only when β r is larger than the largest root β (+) r of ∆ = 0, defined implicitly in (5.8) * * When γs = 0, the largest root is reduced to (d) In the last case we must show that the i /i = 0 isocline does not intersect the interior of D. Equivalently, we must show that in each of the two subcases β < β r ≤ Λ + γ + ν i β r ≤ β ≤ Λ + γ + ν i none of the points A, B, C belongs to D. This is a tedious computation, not reproduced here. For a quick check, we offer a Mathematica file on our website (we rely mostly on FindInstance with empty output to show that certain cases do not exist, and on the command Reduce to decompose other cases into subcases). 6 The boundary case γ s = ν r = 0 The case γ s = ν r = 0, generalizes on the particular case γ s = ν r = 0 = γ r , which was called SIRI model in [DvdD93, Sec. 4 ]. Note first that in this case s df e = 1. It follows that the 4, 5 cases (1(d) and 2(a)) may not arise,and the case 6 (2(b)) becomes degenerate since s df e = 1. It becomes now possible that the hyperbola does not intersect the domain. This is equivalent to its slope at the DFE i (s df e ) = − Λ+γr β−νi+γr being either positive or less than −1, and further equivalent to β ≤ Λ + ν i . This further implies R 0 < 1, and, together with the absence of endemic points and of periodic solutions, leads to the fact that the DFE is the global attractor. After having dealt with this case, which includes also the concave case β > ν i , one may restrict to the case when the hyperbola does intersect the domain, which is equivalent to its slope at the DFE being such that −1 < i (s df e ) < 0 ⇔ β > Λ + ν i (note this implies that its center is in the third quadrant). The proof becomes simpler than in the previous section. For example, the point D of intersection of the hyperbola with i = 1 − s satisfies and hence belongs to D. The roots of ∆ = 0 simplify now to: The pedagogical model is defined as follows: (7.1) The following properties hold: may be shown to be positively invariant with respect to (7.1); therefore, this region must include an attractor set [MLH92, VDL09] . 2. The DFE equilibrium point of the pedagogic system is obtained plugging i df e = 0 in (7.1). Solving with respect to ( s, r) the remaining first and third equation , 0, Λγ s µγ r + (µ + ν r ) (µ + γ s ) . (7.3) In particular, the DFE for the FA model is obtained by substituting Λ = µ in (7.3): , 0, Λγ s Λγ r + (Λ + ν r ) (Λ + γ s ) . (7.4) Note the relation ν r γ s = 0 =⇒ i df e + r df e + s df e = 1. 3. When β r > 0 there may be two endemic equilibrium points (we omit their complicated expressions), but when β r = 0 there is a unique FA EE, obtained by plugging in the first and third equation, and solving the linear system where z := R γr(Λ+νi)+(Λ+νr)(γ+Λ+νi) . Note the intriguing simplification of i (EE) , which shows that and it may be shown that the endemic point belongs to the domain iff 1 ≤ Rs df e . This gives a pre-warning on the role of the parameter Rs df e := R 0 in the stability of the DFE. Example 1 In the particular case νr = γr = γs = 0 =⇒ s df e = 1, we recover the SIR-FA example, for which the sharp threshold property holds [SvdD13, (4.1)], i.e. the disease free equilibrium is globally stable iff R 0 = R s df e ≤ 1. The endemic point simplifies to We may observe that the endemic point X ( Lemma 1 If only strict inequalities between Λ + γ + ν i , βr and β are allowed, then the following ten cases β < Λ + γ + ν i < βr, γs ≥ γ * s βr < Λ + γ + ν i < β, ν i < β β < ν i < Λ + γ + ν i < βr ν i < β< Λ + γ + ν i < βr, γs < γ * s , ∆ < 0 ν i < β< Λ + γ + ν i < βr, γs < γ * s , ∆ ≥ 0 β < βr < Λ + γ + ν i βr < β < Λ + γ + ν i , form a disjoint decomposition of the parameter space. Proof: We note first that in all cases where equality is allowed, the equality case may be arbitrarily assigned to any of the two cases it separates. It is enough therefore to consider only strict inequalities in this Lemma. We want to show that the union of the ten cases equals the union of the six cases representing the possible orders of β, β r , Λ + γ + ν i , which are . Now the cases 1, 2, with Λ + γ + ν i < min[β, β r ], and the cases 9, 10, with Λ + γ + ν i > min[β, β e ] appear only once in the ten cases of theorem 5, as case 1.(a), 1(b), and 2(d)(i-ii). Next, we may check that the union of the cases 3 and 5 (1.(c) and 2.(a) in the Theorem) form together the permutation 3. This requires checking that the other two of the four formal cases taking into account the possible orders between β, ν i and γ s , γ * s are void; the tedious verification is included in the Mathematica file available on our website. To conclude, it remains to check that the cases 4,6,7,8 (i.e. 1.(d), 2.(b), 2.(c)(i) and 2.(c)(ii) in the Theorem) form together a partition of the permutation 4. Note first that cases 7 and 8 may be combined in ν i < β < Λ + γ + ν i < β r , γ s < γ * s . Next, we show in the Mathematica file that the case 4 β < Λ + γ + ν i < β r , γ s ≥ γ * s is incompatible with β < ν i , and so we can modify the case 4 to ν i < β < Λ + γ + ν i < β r , γ s ≥ γ * s . So, 4, 6 and 7-8 become      ν i < β < Λ + γ + ν i < β r , γ s ≥ γ * s β < ν i < Λ + γ + ν i < β r ν i < β< Λ + γ + ν i < β r , γ s < γ * s , whose union is clearly the permutation 4. Lemma 2 A) A necessary and sufficient condition for having precisely one endemic point with s ∈ (0, 1) is C(A+B +C) < 0, where A, B, C are defined in (5.1) . B) Necessary and sufficient conditions for having precisely two endemic points with s ∈ (0, 1) are ∆ > 0 and Proof: The conditions for having two roots bigger than 0 are B A < 0 < C A , and the conditions for having two roots smaller than 1 are obtained by applying these, after substituting y = 1 − x, yielding the result. Since expressing these simple conditions in terms of the parameters of the model turned out quite difficult, we did not finalize this approach. 9 Apendix B: A short review of deterministic epidemic models 9.1 What is a deterministic epidemic model? To put in perspective our work, we would like to start by a definition of deterministic epidemic models, lifted from [KS08] . Definition 2 A deterministic epidemic model is a dynamical system with two types of variables x(t) := ( i(t), z(t)) ∈ R N + , where 1. i(t) model the number (or density) in different compartments of infected individuals (i.e. latent, infectious, hospitalized, etc) which should ideally disappear eventually if the epidemic ever ends; 2. z(t) model numbers (or densities) in compartments of individuals who are not infected (i.e. susceptibles, immunes, recovered individuals, etc). The system must admit an equilibrium called disease free equilibrium (DFE), and hence a "quasi-triangular" linearization the form where D is some forward-invariant subset, where "quasi-triangular" refers to the fact that the functions A i,i , A z,i , Az,z depend on all the variables x(t), and where N is the dimension of x(t). As shown in [KS08] , any epidemic model admitting an equilibrium point ( 0, z df e ) admits the representation (9.1), under suitable smoothness assumptions. In what follows, we will call the point x df e = ( 0, z df e ) a disease free equilibrium (DFE). Remark 24 Note that the essential feature of (9.1) is the "factorization of the disease equations". One of the central objectives in mathematical epidemiology is to study the stability of DFE, i.e. the conditions which make possible eradicating the sickness. For simple models, this amounts, independently of the initial state, to verifying that a famous threshold parameter called basic reproduction number R 0 is less than 1 (and the fundamental problem of interest is, when R 0 > 1, to offer control strategies which force the epidemic to reach the DFE). The basic reproduction number or "net reproduction rate" R 0 is a pillar concept in demography, branching processes and mathematical epidemiology -see the introduction of the book [Bac21]. 1. The notation was first introduced by the father of mathematical demography Lotka [Lot39, Die93] . In epidemiology, the basic reproduction number models the expected number of secondary cases which one infected case would produce in a homogeneous, completely susceptible stochastic population, in the next generation. As well known in the simplest setup of branching process, having this parameter smaller than 1 makes extinction sure. The relation to epidemiology is that an epidemics is well approximated by a branching process at its inception, a fact which goes back to Bartlett and Kendall. 2. With more infectious classes, one deals at inception with multi-class branching processes, and stability holds when the Perron-Frobenius eigenvalue of the "next generation matrix " (NGM) of means is smaller than 1. 3. For deterministic epidemic models, it seems at first that the basic reproduction number R 0 is lost, since the generations disappear in this setup -but see [Bac21, Ch. 3 ], who recalls a method to introduce generations which goes back to Lotka, and which is reminiscent of the iterative Lotka-Volterra approach of solving integro-differential. At the end of the tunnel, a unified method for defining R 0 emerged only much later, via the "next generation matrix" approach [DHM90, vdDW02, VdDW08, DHR10, Per18] . The final result is that local stability of the disease free equilibrium holds iff the spectral radius of a certain matrix called "next generation matrix", which depends only on a set of "infectious compartments" i (which we aim to reduce to 0), is less than one. This approach works provided that certain assumptions listed below hold § . 1. The foremost assumption is that the disease-free equilibrium ( 0, z df e ) is unique and locally asymptotically stable within the disease-free space i = 0, meaning that all solutions of z (t) = ( z(t) − z df e )A z,z ( 0, z(t)), z(0) = z 0 must approach the point z df e when t → ∞). 2. Other conditions are related to an "admissible splitting" as a difference of two parts F, V, called respectively "new infections", and "transitions" Definition 3 A splitting i (t) = F( i(t), z(t)) − V( i(t), z(t)) will be called admissible if F, V satisfy the following conditions [VdDW08, SvdD13] : F( 0, z(t)) = V( 0, z(t)) = 0 F( i(t), z(t)) ≥ 0, ∀( i(t), z(t)) V j ( i(t), z(t)) ≤ 0, when i j = 0, n j=1 V j ( i(t), z(t)) ≥ 0, ∀( i(t), z(t)). , where the subscript j refers the j'th component. The splitting of the infectious equations has its origins in epidemiology. Mathematically, it is related to the "splitting of Metzler matrices"-see for example [FIST07] . Note however that the splitting conditions may be satisfied for several or for no subset of compartments (see for example the SEIT model, discussed in [VdDW08] , [Mar15, Ch 5] ). Unfortunately, for deterministic epidemic models, there is no clear-cut definition of R 0 [Rob07, LB, Thi18]. ‡ 3. We turn now to the last conditions, which concern the linearization of the infectious equations around the DFE. Putting L = A i,i ( 0, z df e ), and letting f denote the perturbation from the linearization, we may write: The "transmission and transition" linearization matrices F, V must satisfy that F ≥ 0 componentwise and V is a non-singular M-matrix, which ensures that V −1 ≥ 0. ‖ Under the assumptions above, the next generation matrix method concludes that R 0 := λ P F (F V −1 ) is a threshhold parameter in the following sense [VdDW08, Thm 1]: 1.(a) When R 0 < 1, the DFE is locally asymptotically stable, while (b) when R 0 > 1, the DFE is unstable. We will call this the "weak R 0 alternative". § And so R0 is undefined when these assumptions are not satisfied. ‡ A possible explanation is that several stochastic epidemiological models may correspond in the limit to the same deterministic model. ‖ 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 [DlSNAQI19, Thm. 1-3]. Also, global asymptotic stability of the DFE holds when R 0 ≤ 1, provided the "perturbation from linearity" f = i(F − V ) − F + V is non-negative On matrix-sir arino models with linear birth rate, loss of immunity, disease and vaccination fatalities, and their approximations A simple planning problem for COVID-19 lockdown Infectious diseases of humans: dynamics and control Optimal control of a sir epidemic with icu constraints and target objectives Un modèle mathématique des débuts de l'épidémie de coronavirus en france Reactive social distancing in a SIR model of epidemics such as COVID-19 Mathematical models in epidemiology Analysis of a disease transmission model in a population with varying size A method for proving the non-existence of limit cycles COVID-19 pandemic control: balancing detection policy and lockdown intervention under ICU sustainability How long should the COVID-19 lockdown continue? The optimal lockdown intensity for COVID-19 Optimal COVID-19 epidemic control until vaccine deployment, medRxiv On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations The construction of next-generation matrices for compartmental epidemic models The estimation of the basic reproduction number for infectious diseases How does transmission of infection depend on population size, Epidemic models: their structure and relation to data Optimal timing of one-shot interventions for epidemic control 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 Problemi di controllo in epidemiologia matematica e comportamentale A disease transmission model in a nonconstant population A light introduction to modelling recurrent epidemics Epidemiological models and Lyapunov functions A feedback SIR (fSIR) model highlights advantages and limitations of infection-based social distancing Slow invariant manifolds of slow-fast dynamical systems Global dynamics of a staged progression model for infectious diseases Hopf bifurcation in epidemic models with a latent period and nonpermanent immunity Qualitative analyses of communicable disease models Balancing quarantine and selfdistancing measures in adaptive epidemic networks Differential equations,dynamical systems and linear algebra A geometric analysis of the SIR, SIRS and SIRWS epidemiological models Optimal control of an SIR epidemic through finite-time non-pharmaceutical intervention Asymptotic analysis of two reduction methods for systems of chemical reactions A contribution to the mathematical theory of epidemics Computation of threshold conditions for epidemiological models and global stability of the disease-free equilibrium (DFE) The failure of R0 Global dynamics of a SEIR model with varying total population size Global asymptotic stability for the seirs models with varying total population size Qualitative analyses of SIS epidemic model with vaccination and varying total population size Analyse démographique avec application particulièreà l'espèce humaine An introduction to mathematical epidemiology Dynamic models of infectious diseases as regulators of population sizes Trends in biomathematics: Modeling cells, flows, epidemics, and the environment Optimal epidemic suppression under an ICU constraint An introduction to the basic reproduction number in mathematical epidemiology A model for a vector-borne disease with control based on mosquito repellents: A viability analysis Multiple equilibria for an SIRS epidemiological system The pluses and minuses of 0 On COVID-19 modelling Global analysis of an SEIR model with varying population size and vaccination Ramsès Djidjou-Demasse, Christian Selinger, Yannis Michalakis, and Samuel Alizon, Epidemiological monitoring and control perspectives: application of a parsimonious modelling framework to the COVID-19 dynamics in France Dynamical systems and population persistence Global stability of infectious disease models using Lyapunov functions Mathematics in population biology Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission Further notes on the basic reproduction number Constructions of Lyapunov functions for classic SIS, SIR and SIRS epidemic models with variable population size, Foro-Red-Mat: Revista electrónica de contenido matemático Global analysis for a general epidemiological model with vaccination and varying population We thank N. Bacaer for providing the references [Bac21, Lot39] , and we thank the referees for their thorough reviews and suggestions.