Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Parity-Dependent Quantum Phase Transition in the Quantum Ising Chain in a Transverse Field
Next Article in Special Issue
A Transcriptome- and Interactome-Based Analysis Identifies Repurposable Drugs for Human Breast Cancer Subtypes
Previous Article in Journal
Mapping Topology of Skyrmions and Fractional Quantum Hall Droplets to Nuclear EFT for Ultra-Dense Baryonic Matter
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

New Results and Open Questions for SIR-PH Epidemic Models with Linear Birth Rate, Loss of Immunity, Vaccination, and Disease and Vaccination Fatalities

1
Laboratoire de Mathématiques Appliquées, Université de Pau, 64000 Pau, France
2
Département des Mathématiques, Université Ibn-Tofail, Kenitra 14000, Morocco
3
Department of Mathematics and Informatics, Polytechnic University of Bucharest, 062203 Bucharest, Romania
*
Author to whom correspondence should be addressed.
Symmetry 2022, 14(5), 995; https://doi.org/10.3390/sym14050995
Submission received: 31 March 2022 / Revised: 5 May 2022 / Accepted: 7 May 2022 / Published: 12 May 2022

Abstract

:
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 “ 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 R 0 < 1 , where 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 R 0 > 1 , there exists a unique fixed endemic point, and 2. the endemic point is locally stable when 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 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.

1. Introduction

Motivation. One of the hardest challenges facing epidemic models is dealing with models within which death is possible, in which the total population N varies, and in which the infection rates depend on N (as is the case in reality, except for a short period of time at the start of an epidemic). As these features confront the researcher with challenging behaviors (the first being that the uniqueness of the fixed endemic point may stop holding), sometimes hard to explain epidemiologically, it seems natural to attempt to identify the simplest class of realistic models for which a theory may be developed. The natural choice is “standard incidence rates”—see (1), since models with nonlinear infection rates are quite complex—see for example [1,2,3,4,5] for the very complex dynamical behaviors which may arise otherwise.
The next issue is choosing the type of birth function b ( N ) to work with. The easiest case is when b ( N ) is constant, but this corresponds to immigration rather than birth, and so our favorite are linear birth rates b ( N ) = Λ N . A bonus for this choice, as well known, is that normalization by the total population leads to a model with constant parameters, which looks similar to classic constant population models, but involves some extra nonlinear terms—see (2). The well-studied classic models may then be recovered via a heuristic “first approximation” (FA) of ignoring the extra terms. This approximation, which deserves to be investigated rigorously via slow-fast/singular perturbation/homogenization techniques, has the merit of putting under one umbrella constant and varying population models with linear birth rates.
At this point, let us mention that we believe that epidemic models should be ideally parameterized by the two matrices F , V , which intervene in the next generation matrix approach, which has been called disease-carrying and state evolution matrices [6,7]. A foundational paper in this direction is [8], Arino shows that further simplifications arise for models having only one susceptible class and disease-carrying matrix of rank one. The first fundamental question, the uniqueness of the endemic point when R 0 > 1 , may be resolved explicitly for the “FA approximation” (and hence also for “small perturbations”, numerically). This has motivated us to propose [9,10] to develop the theory of this class of models, which we call “Arino” or “SIR–PH” models.
Contributions and contents. The goal of our paper is to draw attention to two interesting open problems, for the SIR-PH, 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 of [11,12].
The SAIRS model (2) is presented in Section 2. The history of the problem and some oversights and errors in the literature are recalled in Section 2.1. The basic reproduction number R 0 and the weak R 0 alternative for the DFE equilibrium are established via the next generation matrix approach in Section 2.2. The local stability of the endemic point when the basic reproduction number satisfies R 0 > 1 for the FA model is established in Section 2.3.
A review of the theory of SIR-PH models is provided in Section 3, and some new results in Section 4.
The scaled SAIRS model is revisited in the Section 5, where some previous results in the literature are corrected and completed.

2. The SAIRS Model with Linear Birth Rate

In this paper, we consider a ten parameters SAIR (also called SEIR in the classic literature [13]) epidemic model inspired by [12,14,15,16,17,18,19,20], which we call SAIR/V+S (or SAIR for short) as it groups together immunized people in an R/V compartment. The letter A (from asymptomatic) stands for the fact that the individuals in this compartment may infect the susceptibles. This important feature, already present in [13], was further studied in [11,12,21]. The model studied in [12] is the most complete in the sense that it misses only one of the parameters of interest, namely the important extra death rate due to the disease—see for example [22]. The explanation for this omission in [12] lies probably in the fact that this paper follows the tradition of the “short term constant population epidemics”, in which the total population N is assumed constant, and the endemic point is unique and easy to find explicitly. Our paper investigates, for the model of [12], the topic of the first six papers cited, i.e., we attempt to deal with the variation of N ( t ) during epidemics that may last for a long time and which may never be totally eradicated. We generalize the uniqueness of the endemic point and the local stability results of these six papers and, at the same time, draw attention to certain unnecessary assumptions and mistakes. The hardest issue of global stability is only illustrated via some numerical simulations.
SAIR Model. Letting S ( t ) , E ( t ) , I ( t ) , R ( t ) , D ( t ) and D e ( t ) represent Susceptible individuals, Exposed individuals, Infective individuals, Recovered individuals, naturally Dead individuals and Dead individuals due to the disease, respectively, the model we consider is:
S ( t ) = Λ N ( t ) S ( t ) β i N ( t ) I ( t ) + β e N ( t ) E ( t ) + γ s + μ + γ r R ( t ) , E ( t ) = S ( t ) β i N ( t ) I ( t ) + β e N ( t ) E ( t ) ( γ e + μ ) E ( t ) , γ e = e i + e r I ( t ) = e i E ( t ) I ( t ) γ + μ + δ , R ( t ) = γ s S ( t ) + e r E ( t ) + γ I ( t ) ( γ r + μ ) R ( t ) , D ( t ) = μ ( S ( t ) + E ( t ) + I ( t ) + R ( t ) ) : = μ N ( t ) , D e ( t ) = δ I ( t ) , N ( t ) = S ( t ) + E ( t ) + I ( t ) + R ( t ) N ( t ) = ( Λ μ ) N ( t ) δ I ( t ) .
Epidemiologic meaning of the parameters:
1.
Λ R 0 and μ R 0 denote the average birth and death rates in the population (in the absence of the disease), respectively.
2.
The parameters β i R 0 and β e R 0 denote the infection rates for infective and exposed individuals, respectively;
3.
γ s R 0 is the vaccination rate, γ e R 0 is the rate at which the exposed individuals become infected or recovered, γ 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.
4.
e i R 0 , e r R 0 , are rates of transfer from E to I and R, respectively.
5.
δ R 0 is the extra death rate in the infected compartment due to the disease.
Some particular cases.
1.
If β e = e r = 0 , we obtain a SEIRS type model. If furthermore γ r = γ s = 0 , we arrive to the model masterly studied in [15]; if only one of these parameters is zero, we arrive to the models studied in [16,17,18,19].
2.
If β i = γ r = γ s = 0 = δ , we obtain the SIQR model [23,24].
Remark 1.
Note the notation scheme employed above, which could be applied to any compartmental model. A linear rate of transfer from compartment m to compartment c is denoted by m c , and the total linear rate out of m is denoted by γ m , which implies c m c = γ m , for example γ e = e i + e r . Extra death rates due to the epidemics in a department c are denoted by δ c . An exception is made though for the infectious compartment i , where we simply use the classic notations γ , δ , instead of γ i , δ i . Our scheme would simplify a lot, if adopted, perusing the rather random notations used in the literature.
The scaled model. It is convenient to reformulate (1) in terms of the normalized fractions s = S N , e = E N , i = I N , r = R N . Using N ( t ) = ( Λ μ ) N ( t ) δ I ( t ) , this yields the following nine parameters SEIRS epidemic model (the common death rate μ simplifies, and an extra δ i c appears in the equation of each compartment c–see for example [25] for similar computations).
s ( t ) = Λ s ( t ) β i i ( t ) + β e e ( t ) + γ s + Λ + γ r r ( t ) + δ s ( t ) i ( t ) e ( t ) i ( t ) = [ s ( t ) β e β i 0 0 + ( γ e + Λ ) 0 e i Λ + γ + δ + δ i t I d ] e ( t ) i ( t ) r ( t ) = γ s s ( t ) + e r e ( t ) + γ i ( t ) ( γ r + Λ ) r ( t ) + δ i ( t ) r ( t ) .
Remark 2.
Note that we have written the “infectious” middle equations to emphasize first the factorized form, similar to that encountered for Lotka-Volterra networks—see for example [26]. Secondly, for the factor appearing in these equations, we have emphasized the form
F V + δ i I 2 : = s B V + δ i I 2 .
Note that the matrices F , V are featured in the famous the (Next Generation Matrix) approach [13], that some authors refer to them as “new infections” and “transmission matrices”, that [7] call them the disease-carrying and state evolution matrices, and that ([27], Ch 5) gives a way to define an associated stochastic birth and death model associated to these matrices. The matrix B is further useful in defining and studying more general SIR-PH models—see [9,10] and below, and see also [28], ([29], (2.1)), [6] for related works.
Since the computations will become soon very cumbersome, we will start using from now on the following notations:
v 0 = Λ + γ r + γ s , v 1 = Λ + γ e , v 2 = Λ + γ + δ .
Note that v z , w i t h z = { 1 , 2 } , are the diagonal elements of V, and that v 0 appears as denominator in the DFE (9).
Figure 1 compares the qualitative behavior and equilibrium points of the ( s , i ) coordinates of the three variants of a SIR-type example (discussed in detail in [25]).
We restrict our study to the biological feasible region
D 4 = ( s , e , i , r ) R + 4 , s + e + i + r = 1 .
Its invariance is ensured by the fact that
( s + e + i + r ) = ( δ i Λ γ r ) ( s + e + i + r 1 ) ,
hence s 0 + e 0 + i 0 + r 0 = 1 implies s + e + i + r = 1 , for all t > 0 .
We will work in three dimensions by eliminating r = 1 s e i . The first equation of (2) changes then, and the system becomes
s ( t ) = Λ + γ r γ r ( e ( t ) + i ( t ) ) s ( t ) = + + Λ + γ r + γ s + ( β i δ ) i ( t ) + β e e ( t ) , e ( t ) i ( t ) = [ β e s ( t ) v 1 β i s ( t ) e i v 2 + δ i ( t ) I 2 ] e ( t ) i ( t ) ( s , e , i ) D 3 : = { ( s , e , i ) R + 3 , s + e + i 1 } .
Concerning fixed points, we note first, following ([15], (4.2)), ([17], (13)), the following necessary condition at a fixed endemic point:
s + e + i = 0 e ( γ e e i ) + γ i + s γ s = γ r r + Λ r δ i r   ( 1 e + i + s ) ( Λ + γ r δ i ) = e γ e e i + γ i + s γ s i < i c : = min [ 1 , ( Λ + γ r ) / δ ] ,
where the last implications follows since the right hand side of the last equation in (7) is positive, due to e r = γ e e i 0 .
Thus, the search for fixed points may be reduced to the domain
D f : = D 3 { i : i < i c = min [ 1 , ( Λ + γ r ) / δ ] } .
As a quick preview of the next generation matrix “factorization” idea, we note now that the disease free equation s = Λ s ( Λ + γ s ) + γ r ( 1 s ) (defined by e = i = 0 , r = 1 s ) has fixed point
s d f e = Λ + γ r Λ + γ r + γ s r d f e = 1 s d f e = γ s Λ + γ r + γ s ,
whenever Λ + γ r + γ s 0 , which will be assumed from now on. Note this formula contains only three of the parameters.
The other 6 parameters γ , δ , β i , β e , e i , γ e intervene in the “infection equations” for e , i , and will determine the basic reproduction number (12)
R 0 = s d f e Λ + γ e = + + β e + β i e i Λ + γ + δ : = s d f e R ,
where R is precisely the R 0 which would be obtained in the absence of vaccination.
When δ = 0 , our model (6) reduces to the eight parameters constant population model studied in [12], which has a unique endemic point, with coordinates expressible in terms of R 0 ; for example s e e = s d f e R 0 : = 1 R . Its global stability is, however, hard to prove ([12], Conjecture 15), and the authors resolve only the case β e < γ .

2.1. Some History of the SEIRS and SAIRS Varying Population Models

A first important paper on the varying population SEIRS model (recall that both notations SEIRS and SAIRS have been used already for the same model) is [14], where β e = e r = 0 and where a proportion q of vaccines is allocated to the newborn. We took for simplicity q = 0 , as this parameter does not modify in an essential way the mathematics involved. Note that δ , Λ , γ e , γ s and γ r are denoted in [14] by α , r, σ , p and ϵ , respectively.
Besides the typical weak R 0 alternative, ([14], Thm 2.3(i)) establishes also global asymptotic stability (GAS) of disease free equilibrium when R 0 < 1 and when either (I) β i γ e Λ + γ e δ , (II) γ s = 0 , or a certain non-explicit condition holds. A second result ([14], Thm 2.3(ii)) establishes the uniqueness of the endemic point when R 0 > 1 , and , γ r min [ γ e , γ ] .
Stimulated by the many open cases left above, in [15], the authors considered the particular case β e = e r = γ r = γ s = 0 .As the authors explain, the difficulties met, especially in the global stability problem, forced them to devise new ingenious methods. (12) becomes now R 0 = β i γ e γ e + Λ Λ + γ + δ (note that β , γ e , Λ , δ and R 0 are denoted in [15] by λ , ϵ , b, α and σ , respectively, and that s d f e = 1 ). Building on previous results in the limiting cases δ = 0 for SEIR models with constant population [30] and γ e for SIR models with varying populations [31], they establish the uniqueness and global stability of the endemic point ([15], Cor. 6.2, Thm6.5) when R 0 > 1 .
Note that the Corollary is under the assumption of non-existence of non-constant periodic solutions and the Theorem requires the additional assumption δ < γ e ; thus, the complementary case was left open.
Next, the reference [16] studied the case of waning immunity, leaving again open many cases. Note that the so-called “geometric approach to stability” method, initiated in [15,16,32], was used many times afterward —see for example [19].
Eleven years later, the authors in [17] re-attacked the [14] problem with β e = e r = γ r = 0 and vaccination γ s 0 (denoted by a σ ). We note, by adding the two infection rates in ([17], Figure 1), which yields λ ( 1 a σ ) , that the model depends only on a σ , and so introducing two parameters for γ s unnecessarily complicates the mathematics. Furthermore, the infection rate λ ( 1 a σ ) might as well be denoted by one parameter, and we chose the classical β . Finally, one can relate their results to ours by substituting β with γ and setting σ = 1 (i.e., giving up the second source of infections), and λ = β / ( 1 γ s ) . Then, ([17], (8))
R 0 = s d f e β i γ e Λ + γ e Λ + γ + δ : = R s d f e ,
where
R = β i γ e Λ + γ e Λ + γ + δ and s d f e = Λ Λ + γ s .
See also (12).
Remark 3.
Let us note two open problems left in [17].
1.
([17], Thm 2.1) proves the global stability of the DFE when
R 1
(by using the Lyapunov function γ e e + ( γ e + Λ ) i , which is different from the Lyapunov function used in [33], e ( γ + Λ + δ ) + γ e i . This leaves open the case R 0 < 1 < R , when the DFE is locally, but perhaps not globally stable, and also the question of choosing the “most performant” Lyapunov function.
2.
([17], Thm 2.1) establishes global stability of a unique endemic point in certain cases and suggests that outside those cases, “there may exist stable periodic solutions”.
Seven years later, the reference [18] revisited the [16] problem with γ s = 0 and waning immunity γ r 0 (denoted by ρ ). The authors remark that in this particular case, there are still open questions: “We have not proved, but we strongly believe that if R 0 > 1 , then the system (2) has one and only one endemic equilibrium and that this equilibrium is globally asymptotically stable”. The [16] problem was revisited then in [19], who claims to have also removed the restriction δ < γ e of [15]; however, their crucial equation ([19], (3.28)) is wrong, see Section 5. The fact that the papers [17,19] contain unnecessary conditions, mistakes, and wrong conjectures suggests that for complicated epidemic models with infection rates depending on N, ingenuity stops being enough, and it must be accompanied by verifications via symbolic software, rendered public, as is the case with our paper (in line with the so-called “reproducible research” [34]).
Finally, the authors in [12] apply the geometric approach to global stability, taking into account vaccination and loss of immunity, and also the possibility that the exposed, or rather the asymptomatic, are infectious, but for a simplified constant population model with δ = 0 , Λ 0 . The basic reproduction number (12) becomes now ([12], (5))
R 0 = [ β e + β i e i Λ + γ ] s d f e γ e + Λ
(with the correspondence Λ μ , β e β A , e i α , e r δ A , γ δ I , γ r γ , γ s δ ). The global stability is proven only in the case β e < γ .
These multiple but only partially successful attempts at solving the [14] problem motivated us to revisit it. (Yet another motivation was the literature on bifurcation analysis of SIR and SEIR-type models [1,3,5,35,36,37], and [2,20] which contains yet another plethora of open problems.)

2.2. Warm-Up: The Weak R 0 Alternative for the DFE Equilibrium

Note that the DFE (9) is a stable point for the disease-free equation
s ( t ) = Λ + γ r s ( t ) Λ + γ r + γ s .
Then, we may check that the quadratic—linear decomposition of
e ( t ) i ( t ) = β e s ( t ) + δ i ( t ) ( γ e + Λ ) β i s ( t ) e i δ i ( t ) v 2 e ( t ) i ( t ) = F ( e , i ) V ( e , i ) ,
where
F ( e , i ) = β e s e + β i s i + δ e i δ i 2 , V ( e , i ) = ( Λ + γ e ) e e i e + i v 2
satisfies the splitting assumptions of the NGM (Next Generation Matrix) method [13].
The partial derivatives at the DFE of F ( e , i ) , V ( e , i ) are
F = s d f e β e β i 0 0 , V = v 1 0 e i v 2 ,
and
F V 1 = β e s d f e β i s d f e 0 0 1 v 1 0 e i v 1 v 2 1 v 2 = s d f e [ β e v 1 + β i e i v 1 v 2 ] s d f e β i v 2 0 0 .
We may conclude then by the well-known next generation matrix result of [13] that:
Proposition 1.
The SAIR/V+S basic reproductive number is
R 0 : = s d f e β t v 1 v 2 ,
where v i were defined in (4) and where
β t : = β i e i + β e v 2 .
Furthermore, the weak R 0 alternative (The strong R 0 alternative [33] also covers he case R 0 = 1 .) holds, i.e.,
1.
If R 0 < 1 ( Λ + γ r ) β t < v 0 v 1 v 2 then the disease-free equilibrium is locally asymptotically stable.
2.
If R 0 > 1 ( Λ + γ r ) β t > v 0 v 1 v 2 , the disease-free equilibrium is unstable.
Remark 4.
The formula for R 0 follows also from the [8] formula for SIR-PH models with new infection matrix of rank one F = s α β —see [9,10], and below, for the definition of this class. This formula becomes, after including demography parameters Λ , δ ,
R 0 = β V 1 α , V = A + D i a g ( δ + Λ 1 ) ,
In the particular SAIRS case, this formula may be applied with the parameters defining the model, which are
β = β e , β i , α = 1 0 , A = γ e 0 e i γ , δ = 0 , δ .
Finally, we will write
a : = 1 ( A ) = γ e e i , γ in this example ,
for a quantity which will appear often below.—See for example (24) where we used the standard notation in the theory of phase-type distributions.
A direct proof of Proposition (1) is also easy here. Indeed, after reducing (2) to a three order system by r = 1 s e i , the Jacobian matrix is
J = β i + δ i Λ γ r γ s s β e γ r s ( δ β i ) γ r β i s β e + i δ v 1 e δ + β i s 0 e i γ + ( 2 i 1 ) δ Λ ,
(where β i = β e e + β i i ), and
J D F E = v 0 β e s d f e γ r ( δ β i ) s d f e γ r 0 β e s d f e v 1 β i s d f e 0 e i v 2 .
Note the evident block structure { s } , { e , i } (which is the driving idea behind the next generation matrix method), with one negative eigenvalue v 0 .
Now the “infectious” determinant v 1 v 2 s d f e β t is positive iff
v 1 v 2 ( 1 s d f e β t v 1 v 2 ) = v 1 v 2 ( 1 R 0 ) > 0 R 0 < 1 .
This is also the stability condition, since it may be shown that s d f e β t < v 1 v 2 implies also the trace condition β e s d f e v 1 v 2 < 0 , and so the R 0 alternative holds.
Remark 5.
As usual, it is useful to introduce critical vaccination and critical “total contact” (recall β t = β i e i + β e v 2 ) parameters, as the unique solutions of R 0 = 1 with respect to γ s and β t . The critical values are
γ s * = ( Λ + γ r ) [ β t v 1 v 2 1 ] : = ( Λ + γ r ) [ R 1 ] , β t * = v 0 v 1 v 2 Λ + γ r = v 1 v 2 s d f e .
These are particular cases of the SIR-PH Formulas (29).

2.3. The Endemic Point for the FA Approximation, and the Determinant Formula

We consider now the following “first approximation” (FA) of the SIR-PH-FA varying population dynamics
i ( t ) = [ s ( t ) B V ] i ( t ) s ( t ) = Λ s ( t ) β i ( t ) Λ + γ s s ( t ) + γ r r ( t ) r ( t ) = a i ( t ) + s ( t ) γ s r ( t ) γ r + Λ ( s ( t ) , i ( t ) , r ( t ) ) R + 4 ,
where i ( t ) : = e ( t ) i ( t ) , β a n d a are defined in Remarks 4, and V in (11). Here the extra deaths δ due to the disease in state i are kept, but the quadratic interaction terms involving δ i were neglected. Under this approximation, the endemic point and the determinant of the Jacobian have elegant formulas in terms of R 0 , some already discovered and others hidden in the particular cases of [11,12].
Remark 6.
For this approximation, the sum of the variables is constant only if δ = 0 ; therefore, r may not be eliminated, and we must work in four dimensions.
Lemma 1.
(a) For the SAIR model (15) with extra deaths δ, put i e e = e e e i e e and
v 0 : = Λ + γ r + γ s , v 1 : = Λ + γ e , v 2 : = Λ + γ + δ , v 3 : = ( v 1 + γ r ) v 2 + γ r e i .
Then, the following formulas hold at the endemic point:
s e e = s d f e R 0 = 1 R , e e e i e e = Λ + γ r Λ v 3 + γ r δ e i ( 1 1 R 0 ) v 2 e i , β i e e = β t Λ + γ r Λ v 3 + γ r δ e i ( 1 1 R 0 ) = Λ ( R 0 1 ) v 0 v 1 v 2 Λ v 3 + γ r δ e i = : E I , a i e e = Λ + γ r Λ v 3 + γ r δ e i ( 1 1 R 0 ) = + + γ e γ + ( Λ + δ ) e r , r e e = γ s Λ + v 2 + e i ( Λ + δ ) + Λ ( R 1 ) v 2 e r + γ e i × 1 R Λ v 3 + γ r δ e i .
(b) All the coordinates are positive iff R 0 > 1 .
(c) The vector i e e checks the general SIR-PH normalization Formula (36).
Proof. 
(a) One can do a direct computation or apply ([10], Prop. 2), a particular case of which is included for completeness as Proposition 3, Section 3. That result is expressed in terms of the Perron Frobenius eigenvector of the matrix M = 1 R B V (for the 0 eigenvalue); in our case, this is v 2 e i , see second equation in (16).
(b) Is obvious.
(c) May be easily checked, combining the last two rows in (16). □
Remark 7.
Several particular cases of this problem have been studied in the literature. The case δ = 0 is studied in [12], and the reference [38] considered the case β e = γ r = e r = 0 , δ > 0 , where
α = 1 0 , A = γ e 0 γ e γ , a t : = ( A ) t 1 = 0 γ , β = 0 β i , δ = 0 δ .
The normalization Formula (36) reduces in this case to
β i = Λ [ R 1 s d f e ] = Λ β i γ e ( Λ + γ + δ ) γ e + Λ 1 γ s .
Lemma 2.
The following remarkable identity holds
D e t ( J e e ) = Λ v 0 v 1 v 2 ( 1 R 0 ) = D e t ( J d f e ) ,
where v i are defined in (4).
Proof. 
See the proof of Proposition 3. □
Remark 8.
We conjecture that the endemic point is always locally stable when 1 < R 0 . We have attempted to apply the classic Routh–Hurwitz–Lienard–Chipart–Schur–Cohn–Jury (RH) methods [39,40,41], which are formulated in terms of the coefficients of the characteristic polynomial D e t ( J z I n ) = ( z ) n + a 1 z n 1 + . . . + a n , and of certain Hurwitz determinants H i ([39], (15.22)). At order four, D e t ( J z I 3 ) = z 4 T r ( J ) z 3 + z 2 M 2 ( J ) z M 3 ( J ) + D e t ( J ) , where M 2 , M 3 are the sums of the second and third order principal leading minors of J, and one ends up with ([39], pg. 137)
T r ( J ) < 0 , M 2 > 0 , M 3 < 0 , D e t ( J ) > 0 , 0 < T r ( J ) M 2 M 3 T r ( J ) D e t ( J ) M 3 2 .
Now in our example the determinant is positive when R 0 > 1 by Lemma 2, and the trace, given by
γ r e i ( γ s + Λ ) ( Λ + δ ) + Λ v 2 v 1 ( R ( Λ + γ r ) + Λ ) + γ r ( Λ + γ s ) e i γ r ( Λ + δ ) + Λ v 2 ( γ r + v 1 ) v 2 ( Λ + γ r ) e i β i R v 2
is negative. However, the check of the sum of the second and third-order principal leading minors of the Jacobian at EE and the additional Hurwitz criterion seemed to exceed our machine power. (At order three, RH becomes T r ( J ) < 0 , T r ( J ) M ( J ) < D e t ( J ) < 0 , where M is the sum of the second-order principal leading minors of J, which is considerably simpler.)

3. A Review of Arino and Rank-One SIR-PH Models

3.1. SIR-PH Models with Demography, Loss of Immunity, Vaccination and One Susceptible and One Removed Classes

The fundamental concept of basic reproduction number R 0 can be only defined (as the spectral radius of the next generation matrix) for epidemic models to which the next generation matrix assumptions apply. It seems more practical, therefore, to restrict to “Arino models” where R 0 may be explicitly expressed in terms of the matrices that define the model [8,42,43,44].
The idea behind these models is to further divide the noninfected compartments into (Susceptible) (or input) classes, defined by producing “new non-linear infections”, and output R classes (like D , D e in our first 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 with the linear force of infection as it is known that non-linear forces of infection may lead to very complex dynamics [1,2,5,36,45,46], which are not always easy to interpret epidemiologically. This is in contrast with the Arino models, where typically one may establish the absence of periodic solutions (closed orbits, homoclinic loops, and oriented phase polygons) [47,48].
It is convenient to restrict even further to the case of one removed class (w.l.o.g. ) and only one susceptible class (a significant simplification).
Definition 1.
A “SIR-PH epidemic” of type n, with demography parameters ( Λ , μ ) (scalars), loss of immunity and vaccination parameters γ r , γ s , is characterized by two matrices A , B of dimensions n × n and a column vector of extra death rates δ . This model contains one susceptible class S, one removed state R (healthy, vaccinated, etc), and a n-dimensional vector of “disease” states I (which may contain latent/exposed, infective, asymptomatic, etc). The dynamics are:
S ( t ) = Λ N S ( t ) N β I ( t ) ( γ s + μ ) S ( t ) + R ( t ) γ r , β = 1 B , I ( t ) = [ S ( t ) N B + A D i a g ( δ + μ 1 ) ] I ( t ) , R ( t ) = a I ( t ) + γ s S ( t ) R ( t ) ( γ r + μ ) , a = 1 ( A ) N ( t ) = S ( t ) + 1 I ( t ) + R ( t ) = ( Λ μ ) N I ( t ) δ , D ( t ) = μ ( S ( t ) + 1 I ( t ) + R ( t ) ) , D e ( t ) = I ( t ) δ .
Here,
1.
I ( t ) R n is a row vector whose components model a set of disease states (or classes).
2.
R ( t ) accounts for individuals who recovered from the infection.
3.
B is a n × 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 β the vector containing the sum of the entries in each row of B, namely, β = 1 B .
4.
A is a n × 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 j , satisfies A i , j 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 0 and describes the rate at which individuals in the disease class i leave towards non-infectious compartments. Alternatively, A is a non-singular M-matrix [8,49]. (Recall that an M-matrix is a real matrix V with v i j 0 , i j , and having eigenvalues whose real parts are nonnegative [50].)
5.
δ R n is a row vector describing the death rates in the disease compartments, which are caused by the epidemic.
6.
γ r is the rate at which individuals lose immunity (i.e., transition from recovered states to the susceptible state).
7.
γ s is the rate at which individuals are vaccinated (immunized).
Remark 9.
1.
Note that a t : = ( A ) t 1 is a vector with a well-known probabilistic interpretation in the theory of phase-type distributions: it is the column vector that completes a matrix with negative row sums to a matrix with zero row sums.
2.
A particular but revealing case is that when the matrix B has rank 1 and is necessarily hence of the form B = α β , where α R n is aprobability column vectorwhose components α j represent the fractions of susceptibles entering into the diseased compartment j when infection occurs. We call this case “rank one SIR-PH”, following Riano [49], who emphasized its probabilistic interpretation—see also [51], and see [52] for an early appearance of such models.
It is convenient to reformulate (18) in terms of the fractions normalized by the total population
s = S N , i = 1 N I , r = 1 N R , N = s + 1 i + r .
The reader may check that the following equations hold for the scaled variables:
s ( t ) = Λ Λ + γ s s ( t ) + r ( t ) γ r s ( t ) β δ i ( t ) i ( t ) = [ s ( t ) B + A D i a g δ + Λ 1 + δ i ( t ) I n ] i ( t ) r ( t ) = s ( t ) γ s + a i ( t ) r ( t ) γ r + Λ + r ( t ) δ i ( t ) s ( t ) + 1 i ( t ) + r ( t ) = 1 ,
and the Jacobian, using δ i ( t ) i ( t ) = i ( t ) δ + δ i ( t ) I n , is
J = β i ( Λ + γ s ) s β δ γ r B i s B V + i δ 0 γ s a + r δ ( γ r + Λ ) + δ i I n + 2 .
By letting n : = s + 1 i + r , we have
n ( t ) = ( Λ δ i ( t ) ) ( 1 n ( t ) ) = 0 ;
The above equation guarantees that if s ( t 0 ) + i ( t 0 ) + r ( t 0 ) = 1 for some t 0 R 0 , then s ( t ) + i ( t ) + r ( t ) = 1 for all t t 0 . Accordingly, in what follows we will always assume that n ( t 0 ) = 1 , which guarantees that n ( t ) = 1 , t .
The following definition puts in a common framework the dynamics for the scaled process and two interesting approximations.
Definition 2.
Let Φ s , Φ i , Φ r { 0 , 1 } and let
s ( t ) = Λ Λ + γ s s ( t ) + r ( t ) γ r s ( t ) i ( t ) β + Φ s s ( t ) δ i ( t ) , i ( t ) = [ s ( t ) B + A D i a g δ + Λ 1 ) + Φ i δ i ( t ) ] i ( t ) , r ( t ) = s ( t ) γ s + a i ( t ) r ( t ) D i a g γ r + Λ 1 + Φ r r ( t ) δ i ( t ) , s ( t ) + 1 i ( t ) + r ( t ) = 1 .
1.
The model (22) with Φ s = Φ i = Φ r = 1 will be called scaled model (SM).
2.
The model (22) with Φ s = Φ i = Φ r = 0 will be called first approximation (FA).
3.
The model (22) with Φ s = Φ r = 1 , Φ i = 0 will be called intermediate approximation (IA).
Example 1.
The classic SEIRS model
s ( t ) = Λ s ( t ) β i i ( t ) + γ s + Λ + γ r r ( t ) e ( t ) i ( t ) = ( γ e + Λ ) β i s ( t ) γ e γ + Λ + δ e ( t ) i ( t ) r ( t ) = γ s s ( t ) + γ i ( t ) r ( t ) ( γ r + Λ ) .
is a particular case of SIR-PH-FA model obtained when
α = 1 0 , A = γ e 0 γ e γ , a = 1 ( A ) = 0 γ , β = 0 β i , s o B = 0 β i 0 0 , δ = 0 δ .
The SAIR is obtained by modifying the parameters to
α = 1 0 , A = γ e 0 e i γ , a = 1 ( A ) = e r γ , β = β e β i , s o B = β e β i 0 0 , δ = 0 δ .

3.2. The Eigenstructure of the Jacobian for the SIR-PH Scaled Model

For the scaled model, we can eliminate r. Then, the system becomes then n + 1 dimensional:
s ( t ) = Λ + γ r Λ + γ r + γ s s ( t ) γ r 1 i ( t ) s ( t ) β δ i ( t ) i ( t ) = [ s ( t ) B + A D i a g δ + Λ 1 + δ i ( t ) I n ] i ( t ) s ( t ) + 1 i ( t ) 1 .
The Jacobian matrix of the scaled model is given by
J = ( Λ + γ r + γ s ) β δ i γ r 1 s β δ B i s B + A D i a g ( δ + Λ 1 ) + δ i + δ i I n : = ( Λ + γ r + γ s ) β δ i γ r 1 s β δ B i s B V + δ i + δ i I n ,
and
J d f e = v 0 γ r 1 s β δ 0 s d f e B V ,
where V : = D i a g ( δ + Λ 1 ) A .
Remark 10.
Note the block structure (which suggested probably the next generation matrix approach), and that
D e t ( J d f e ) = v 0 D e t ( s d f e B V ) .
We highlight next a simple but important consequence of the fact that V = D i a g ( δ + Λ 1 ) A is an invertible matrix, especially when B is assumed to have rank 1.
Lemma 3.
(a) When B = α β is of rank 1, the matrix B V 1 has precisely one non-zero eigenvalue.
(b) The remaining eigenvalue equals the trace
T r ( B V 1 ) = β V 1 α : = R
Hence, the Perron-Frobenius eigenvalue of B V 1 is
λ P F B V 1 = R .
Proof. 
(a) Since B = α β has rank 1, the same holds for B V 1 , and the ”rank-nullity theorem”  r a n k ( B V 1 ) + n u l l i t y ( B V 1 ) = n [53] implies that ( n 1 ) of the eigenvalues of B V 1 are zero.
(b) Using the invariance of the trace under cyclic permutations, we conclude that the trace of α β V 1 equals β V 1 α . Since V 1 has only nonnegative entries, this value must be positive and hence the Perron-Frobenius eigenvalue. □

3.3. The Basic Reproduction Number for SIR-PH, via the Next Generation Matrix Method

We follow up here on a remark preceding ([8], Thm 2.1), and via the next generation matrix method [54,55], and show in the following proposition that their simplified formula for the basic reproduction number still holds when the loss of immunity and vaccination are allowed, provided that B = α β has rank one.
Proposition 2.
Consider a SIR-PH model (20), with parameters
( A , B , Λ , δ , γ s , γ r ) .
1.
The unique disease-free equilibrium is ( s d f e , 0 , r d f e ) = Λ + γ r Λ + γ r + γ s , 0 , γ s Λ + γ r + γ s .
2.
The DFE is locally asymptotically stable if R 0 < 1 and is unstable if R 0 > 1 , where
R 0 = λ P F ( F V 1 ) ,
where F = s d f e B , V = D i a g ( δ + Λ 1 ) A (see 26), and λ P F denotes the (dominant) Perron Frobenius eigenvalue.
3.
For B : = α β of rank one, we further have
(a) 
R 0 = s d f e R , where R = α V 1 β .
(b)  
The critical vaccination defined by solving R 0 = 1 with respect to γ s is given by
γ s * : = ( Λ + γ r ) α V 1 b 1 = ( Λ + γ r ) R 1 .
Proof. 
1. The disease free system (with i = 0 , r = 1 s ) reduces to
s ( t ) = Λ ( γ s + Λ ) s ( t ) + ( γ r s ( t ) ) ( 1 s ( t ) ) .
2. It is enough to show that the conditions of ([13], Thm 2) hold.
The DFE and its local stability for the disease-free system have already been checked in the SAIR/V+S example.
We provide now a splitting for the infectious equations:
i ( t ) = [ s ( t ) B + δ i ( t ) I n ] i ( t ) [ D i a g ( δ + Λ 1 ) A ] i ( t ) : = F ( s , i ) V ( i )
(where r = 1 s i 1 ). The corresponding gradients at the DFE i = 0 are
F = F ( X ( D F E ) ) i = s B V = V ( X ( D F E ) ) i = D i a g ( δ + Λ 1 ) A .
We note that F has non-negative elements and that V is an M-matrix, and, therefore, V 1 exists and has non-negative elements, Λ , δ . We may check that the next generation matrix conditions [13] are satisfied.
For example, the last non-negativity condition
i ( t ) [ D i a g ( δ + Λ 1 ) A ] 1 0 , i D ,
is a consequence of A being a M-matrix, which implies A 1 0 , componentwise.
3. (a) Using Lemma 3 and the obvious linearity in s , we may conclude that
R 0 = λ P F ( F V 1 ) = s d f e λ P F ( B V 1 ) = s d f e R .
3. (b) May be easily verified □

3.4. The Endemic Point of the SIR-PH-FA Model

In this section, we give more explicit results for the endemic equilibrium of the following approximate model, referred to as SIR-PH-FA
i ( t ) = [ s ( t ) B V ] i ( t ) s ( t ) = Λ s ( t ) β i ( t ) Λ + γ s s ( t ) + γ r r r ( t ) = a i ( t ) + s ( t ) γ s r ( t ) γ r + Λ ( s ( t ) , i ( t ) , r ( t ) ) R + n + 2 .
Remark 11.
For this approximation, the sum of the variables is constant only if δ = 0 ; therefore, r may not be eliminated.
If R 0 > 1 , then (33) may have a second fixed point within its forward-invariant set. This endemic fixed point must be such that the quasi-positive matrix s e e B V is singular, and that i is a Perron-Frobenius positive eigenvector.
Let i 0 denote an arbitrary positive solution of
( s e e B V ) i 0 : = M i 0 = 0 ,
and let
i = Λ [ 1 s e e 1 s d f e ] β 1 s e e γ r γ r + Λ a i 0 i 0
denote the unique vector of disease components i e e which satisfies also the normalization:
β 1 s e e γ r γ r + Λ a i = Λ [ 1 s e e 1 s d f e ] .
Proposition 3.
Consider a SIR-PH-FA model (33) with parameters ( Λ , A , B , δ , γ s , γ r ) , where A is assumed irreducible, with R 0 > 1 . Then:
1.  
There exists a unique endemic fixed point within its forward-invariant set R + n + 2 iff
β R a γ r γ r + Λ i 0 > 0 .
This fixed point satisfies
(a)  
s e e = 1 R ,
(b)  
that i given by (35) is a Perron-Frobenius positive eigenvector of the quasi-positive singular matrix s e e B V ,
(c)  
the normalization (36).
2.  
The determinant identity D e t ( J d f e ) = D e t ( J e e ) holds.
Proof. 
Recall the fixed point system
0 = ( s e e B V ) i 0 = Λ Λ + γ s s e e + r e e γ r s e e β i , 0 = a i + s e e γ s r e e ( γ r + Λ ) r e e = 1 γ r + Λ a i + γ s s e e
and note that r e e is positive if i , s e e are.
1. Let us examine the two cases which arise from factoring the disease equations. More precisely, we will search separately in the disease free set  { i = 0 } and in its complement. Then: (A) either i { i = 0 } and solving
Λ = Λ + γ s s r γ r 0 = s γ s r ( γ r + Λ ) = s r Λ + γ s γ s γ r ( γ r + Λ )
for s , r yields the unique DFE, or
(B) the determinant of the resulting homogeneous linear system for i 0 must be 0, which implies that s = s e e satisfies
d e t = + + s e e B V = 0 .
Using that V = D i a g ( δ + Λ 1 ) A is an invertible matrix, and d e t ( U U ) = d e t ( U ) d e t ( U ) , (37) may be written as
D e t s e e B V V 1 = 0 D e t B V 1 1 s e e I = 0 .
Thus:
(a) 1 s e e must equal the Perron Frobenius eigenvalue R (recall Lemma 3). Note that s e e < 1 follows from R R 0 > 1 .
(b) i is a Perron-Frobenius eigenvector of the quasi-positive matrix s e e B V , and hence may be chosen as positive.
(c) To determine the proportionality constant, it remains to solve the second equation:
β s e e a γ r γ r + Λ i = Λ Λ + γ s s e e + γ r γ r + Λ γ s s e e β a R γ r γ r + Λ i = Λ R Λ + γ s γ r γ r + Λ 1 = Λ R Λ γ s Λ γ r + Λ = Λ R γ s + γ r + Λ γ r + Λ ,
yielding (36).
2. At the DFE, the infectious equations decouple, and the triangular block structure implies
D e t ( J d f e ) = Λ v 0 D e t ( s d f e B V ) .
For the EE, we will compute the determinant of the Jacobian matrix:
J = s B V B i 0 s β β i ( Λ + γ s ) γ r a γ s ( γ r + Λ ) ,
after applying simplifying row and column operations which preserve the determinant (“Neville eliminations” [56]), to be denoted by ∝.
However, first, we will take a detour through the more explicit SAIRS-FA model (15), where i = e i ,
J = β e s v 1 β i s β e e + β i i 0 e i v 2 0 0 β e s β i s ( Λ + γ s ) β i i β e e γ r e r γ γ s ( Λ + γ r ) s β e v 1 β i s β e e + β i i 0 e i v 2 0 0 v 1 0 ( Λ + γ s ) γ r e r γ γ s ( Λ + γ r )
(here we added row one to row three), and D e t ( J d f e ) simplifies to
D e t ( J d f e ) = Λ v 0 ( v 1 v 2 s d f e β t ) = Λ v 0 v 1 v 2 ( 1 R 0 ) .
The Jacobian at the endemic point is
J e e = β e R v 1 β i R β e e e e + β i i e e 0 e i v 2 0 0 v 1 0 ( Λ + γ s ) γ r e r γ γ s ( Λ + γ r )
requires more work, and we will start by asking Mathematica for the LUDecomposition J = P e r L U D e t ( J ) = D e t ( P e r ) D e t ( U ) . The first factor P e r is a permutation matrix with determinant 1 in our case, the second, L, is lower triangular with one on the diagonal, and the second is upper triangular
U = e i v 2 0 0 0 v 2 γ e e i ( Λ + δ ) e i γ s γ s v 0 0 0 Λ v 0 v 1 v 2 ( 1 R 0 ) γ r e i ( Λ + δ ) + v 2 Λ ( γ r v 1 ) 0 0 0 0 γ r e i ( Λ + δ ) + v 2 Λ ( γ r v 1 ) e i ( Λ + δ ) γ e v 2 ,
with determinant
D e t ( U ) = Λ v 0 v 1 v 2 ( R 0 1 ) = D e t ( J e e ) .
We provide now a second derivation based on determinant preserving transformations
J e e β e R v 1 β i R β e e e e + β i i e e 0 e i v 2 0 0 v 1 0 v 0 γ r e r γ v 0 ( Λ + γ r ) β e R v 1 β i R β e e e e + β i i e e 0 e i v 2 0 0 v 1 0 v 0 γ r e r v 1 γ 0 Λ ,
where we substracted column four from three, and then added row three to row four and where e e e and i e e are defined in (16). We develop now by third column:
D e t ( J e e ) = E I e i v 2 0 v 1 0 γ r e r v 1 γ Λ + v 0 β e R v 1 β i R 0 e i v 2 0 e r v 1 γ Λ = E I e i v 2 0 v 1 0 γ r e r v 1 γ Λ v 0 Λ β e R v 1 β i R e i v 2 = E I e i v 2 0 v 1 γ r 0 γ r 0 ( Λ + δ ) Λ v 0 Λ β e R v 1 β i R e i v 2 = E I e i γ r ( Λ + δ ) + v 2 Λ ( v 1 + γ r ) + Λ v 0 v 2 ( β e R v 1 ) e i β i R = E I Λ v 3 + δ e i γ r
where E I is defined in (16), and the third equality follows by substracting column three from column one and then adding row one to row three. Then using β e = R v 1 v 2 e i β i v 2 , the last product in the equality four cancels, and recall E I = β e e e e + β i i e e = Λ ( R 0 1 ) v 0 v 1 v 2 Λ v 3 + γ r δ e i , this yields
D e t ( J e e ) = Λ v 0 v 1 v 2 ( R 0 1 ) = D e t ( J d f e ) .
Or, by developing the first row, the determinant reads
D e t ( J e e ) = ( β e R v 1 ) v 2 0 0 0 ( Λ + γ s ) γ r γ γ s ( Λ + γ r ) β i R e i 0 0 v 1 ( Λ + γ s ) γ r e r γ s ( Λ + γ r ) + E I e i v 2 0 v 1 0 γ r e r γ ( Λ + γ r ) = ( β e R v 1 ) Λ v 0 v 2 β i R Λ e i v 0 + E I γ γ r e i + v 2 ( v 1 ( Λ + γ r ) e r γ r )     = E I γ γ r e i + v 2 ( v 1 ( Λ + γ r ) e r γ r ) ( b y u s i n g β e = R v 1 v 2 e i β i v 2 ) = E I γ γ r e i v 2 γ r ( v 1 Λ e i ) + v 1 v 2 ( Λ + γ r ) = E I v 1 v 2 Λ + e i γ r ( v 2 γ ) + Λ γ r v 2 = E I Λ ( v 1 v 2 + e i γ r + γ r v 2 ) + e i γ r δ = E I Λ v 3 + e i γ r δ ,
thus, D e t ( J e e ) = Λ v 0 v 1 v 2 ( R 0 1 ) = D e t ( J d f e ) .

4. A Glimpse of the Intermediate Approximation for the SIR-PH Model

The intermediate approximation associated to (20) is
s ( t ) = Λ Λ + γ s s ( t ) + r ( t ) γ r s ( t ) β δ i ( t ) i ( t ) = s ( t ) B V i ( t ) r ( t ) = s ( t ) γ s + a i ( t ) r ( t ) Λ + γ r δ i ( t ) ) ;
Proposition 4.
1.  
The DFE points of the scaled, the intermediate approximation, and the FA are equal, given by ( Λ + γ r Λ + γ r + γ s , 0 , γ s Λ + γ r + γ s ) .
2.  
An endemic point must satisfy that i e e is a positive eigenvector of the matrix s e e B V for the eigenvalue 0 (same as for theFA), that
s e e = 1 R , r e e = γ s R + a i e e Λ + γ r δ i e e > 0 ,
and that
Λ ( R 1 ) γ s + γ s γ r Λ + γ r δ i e e = [ ( β δ ) R γ r Λ + γ r δ i e e a ] i e e .
Since this equation is quadratic (see (41)), we may have a priori two, one, or zero endemic points.
Proof. 
1.
The equations determining the DFE for the three models coincide.
2.
s e e = 1 R has been established in the proof of Proposition 3, and r e e = γ s R + a i e e Λ + γ r δ i e e follows immediately from the last equation in (39). By susbtitution of r e e and s e e into the first equation in (39), we obtain
Λ ( R 1 ) γ s + γ s γ r Λ + γ r δ i e e = [ ( β δ ) R γ r Λ + γ r δ i e e a ] i e e ,
which yields the result.
Remark 12.
1.
Under the substitution δ i e e 0 , (40) reduces to Formula (36).
2.
The existence of positive solutions to (40) requires studying a quadratic equation X λ 2 + Y λ + Z . Indeed, putting i e e = λ i 0 , (40) may also be written as
( Λ + γ r δ λ i 0 ) ( Λ ( R 1 ) γ s ) + γ s γ r = ( Λ + γ r δ λ i 0 ) ( β δ ) λ i 0 R γ r a λ i 0 γ s γ r = ( Λ + γ r δ λ i 0 ) ( β δ ) λ i 0 + γ s Λ ( R 1 ) R γ r a λ i 0 ,
and we find
X = δ i 0 ( β δ ) i 0 , Y = ( γ s Λ ( R 1 ) ) δ + R γ r a + ( Λ + γ r ) ( δ β ) i 0 , Z = γ s γ r + ( Λ + γ r ) ( Λ ( R 1 ) γ s ) .
Example 2. The intermediate approximation of the SAIRS modelis
e ( t ) i ( t ) = ( γ e + Λ ) + β e s ( t ) β i s ( t ) e i γ + Λ + δ e ( t ) i t s ( t ) = Λ s ( t ) β i i t + β e e ( t ) + γ s + Λ δ i t + γ r r ( t ) r ( t ) = γ s s ( t ) + e r e ( t ) + γ i t r ( t ) ( γ r + Λ δ i t ) .
is a particular case of SIR-PH-IA model ( α , A , a , β and δ were defined in (24)).
For the SAIRS model (38) with extra deaths δ, putting i e e ( 1 , 2 ) = e e e ( 1 , 2 ) i e e ( 1 , 2 ) , the following formulas hold at the endemic points:
s e e = s d f e R 0 = 1 R e e e ( 1 , 2 ) = v 2 e i v 0 ( δ R ( Λ + δ ) ) + Λ δ + Λ R Λ + γ s + δ R γ s + v 2 2 Λ R v 0 v 1 + Λ + γ s 2 δ e i δ e i v 1 v 2 R ± v 2 2 v 0 e i ( δ R ( Λ + δ ) ) + δ e i Λ + R γ s + Λ R e i Λ + γ s + v 2 Λ R v 0 v 1 + Λ + γ s 2 + 4 Λ δ e i R γ s v 0 ( R 1 ) v 1 v 2 R δ e i 2 δ e i δ e i v 1 v 2 R i e e ( 1 , 2 ) = v 2 e i v 0 ( δ R ( Λ + δ ) ) + Λ δ + Λ R Λ + γ s + δ R γ s + v 2 2 Λ R v 0 v 1 + Λ + γ s 2 v 2 δ δ e i v 1 v 2 R ± v 2 2 v 0 e i ( δ R ( Λ + δ ) ) + δ e i Λ + R γ s + Λ R e i Λ + γ s + v 2 Λ R v 0 v 1 + Λ + γ s 2 + 4 Λ δ e i R γ s v 0 ( R 1 ) v 1 v 2 R δ e i 2 v 2 δ δ e i v 1 v 2 R
Remark 13.
A yet another open problem is whether the endemic point must always exist for theFAandIAmodels..

5. The Scaled SAIRS Model: Existence, Uniqueness, and Local Stability of the Endemic Point

5.1. Reduction to One Dimension and the Problem of Existence of Endemic Equilibrium

We will follow here the idea of the particular cases ([15], (4.3)), ([17], (14)) and ([19], (3.28)), in which the authors eliminate e , s in the fixed point equations
0 = Λ + γ r ( 1 e i ) s = + + v 0 + ( β i δ ) i + β e e , 0 0 = = + + β e s v 1 β i s e i v 2 + δ i I 2 e i
and study a resulting polynomial equation in i. An alternative to the successive eliminations suggested in [15,17,19] is to notice that a strictly positive endemic point must satisfy that the determinant
Δ = β e s + δ i v 1 β i s e i δ i v 2 = 0 .
After eliminating s from the first equation, the denominator of the determinant is the denominator of s , and the numerator of the determinant is a fifth-order polynomial in i , with factor i .
e = i ( v 2 i δ ) e i , s = Λ + γ r 1 e i v 0 + i ( β i δ ) + e β e .
From the equation of s in (44), we note a typo in ([17], (15)), which should be s = Λ i ( β i δ ) + Λ + γ s .
The next result shows that the existence of endemic points may be reduced for SAIRS to a one-dimensional problem.
Lemma 4.
An endemic point E e e satisfying 0 < i e e < min [ 1 , ( Λ + γ r ) / δ ] and
γ r ( e e e + i e e 1 ) Λ
will satisfy E e e D e .
Proof. 
By (44), 0 < i e e < 1 e e e > 0 . Also, the numerator of s e e , given by Λ + ( 1 i e e e e e ) γ r , is positive by (45). Furthermore, 0 < i e e < i c (recall (7)) implies that the denominator of s e e is positive.
Finally, 0 < 1 s e e e e e i e e follows from (7) and s e e > 0 , e e e > 0 , i e e > 0 . □
The next result shows that the endemic point is unique, without the unnecessary extra conditions assumed in [17].
Proposition 5.
R 0 > 1 ensures the existence and uniqueness of the endemic point for the [17] SIR-PH-SM problem.
Proof. 
For the [17] problem (45) holds trivially (since γ r = 0 , γ s 0 ). Thus, it only remains to show that the third order polynomial p ( i ) which results from the elimination has precisely one root in ( 0 , i c ) , when R 0 > 1 . (The polynomial is of fourth order in general, with a complicated formula, but the leading coefficient is δ 3 β e , and when β e = 0 and δ = 0 , the fourth order polynomial becomes of third and first order, respectively. )It turns out that this polynomial satisfies
p ( 0 ) = ( R 0 1 ) γ e v 1 v 2 Λ + γ s < 0 p ( i c ) = γ e 2 β i γ Λ + δ ( γ + δ ) γ s δ > 0 .
This implies that when R 0 > 1 , p ( i ) must have either one, two or three roots in ( 0 , i c ) .
The last case may be ruled out using an interesting algebraic identity ([15], (4.3)), ([17], (14))
p ( i ) = R 0 f S H ( i ) , f S H ( i ) = f L i ( i ) : = 1 + β i δ γ s + Λ i 1 i δ v 2 1 i δ v 1 .
This shows that solving the equation p ( i ) = 0 (which provides the endemic equilibrium values of i ) is equivalent to equating to R 0 a function with known roots f S H ( i ) .
As a sanity check, note that when i = 0 , this equation reduces to R 0 = 1 , which is consistent with the fact that the endemic point only appears at this threshold (equivalently, when R 0 = 1 , the polynomial p ( i ) = R 0 f S H ( i ) admits only one root, namely i = 0 ).
Refer now to the Figure 2 which plots the function f S H ( i ) . It may be easily checked that the roots of f S H ( i ) , given by i 1 = v 1 δ , i 2 = v 2 δ , i 3 = Λ β i δ , do not belong to the range ( 0 , i c = min [ 1 , Λ δ ] ) (recall that γ r = 0 ). Indeed, i 2 > 1 , i 1 > Λ δ , and R 0 > 1 β t > β t * β i > δ i 3 < 0 . It follows that the largest root of f S H ( i ) = R 0 must be outside the interval ( 0 , Λ δ ) , ending the proof. In ([17], Thm 3.3), the authors used a slightly different approach which leads to several unnecessary cases, as seen above. □
From the proof above, when R 0 > 1 , the line f = R 0 has exactly one intersection ( i e e , f S H ( i e e ) ) with the graph of f S H ( i ) that satisfies i e e [ 0 , i c ] ,– see Figure 2.
Figure 3a displays the bifurcation diagrams of the equilibrium value(s) of i for the [17] problem with respect to β ; we note the usual forward bifurcation diagram when R 0 reaches its critical value 1. Figure 3b displays the bifurcation diagrams of the equilibrium value(s) of i with respect to δ ; we notice that nothing happens in the critical points identified in ([17], Thm 3.3).
Figure 4 compares the qualitative behavior and equilibrium points of the ( s,e+i)-coordinates of the three variants of the SEIRS scaled model discussed in [17], and note that the endemic coordinates s E E F A = s E E I n = 1 / R , which is illustrated by a vertical green line.

5.2. Conjecture for the SEIRS Model

The uniqueness of the endemic point and its local stability when R 0 > 1 was also claimed in [19], and we conjecture that these results are true, even for SEIRS, with β e = e r = 0 , γ s 0 , γ r 0 .
However, this must be viewed as an open problem even in the case β e = e r = 0 = γ s , since the crucial analog of (46), the equation ([19], (3.28)) used intensively in their proof is wrong, see Figure 5.
Indeed, recall that the i equation may be written as
f L u ( i e e ) = R 0 + g L u ( i e e ) ,
where f L u ( i ) = f S H ( i ) (with γ s replaced by 0). The equation ([19], (3.30)) states that g L u ( i ) is rational, which implies that the polynomial p ( i ) is of fourth order, while the correct order of p ( i ) is 3. In fact, ([19], (3.30)) should be replaced by the second order polynomial
g L u ( i ) = R 0 γ r i γ r Λ + 1 γ r ( γ i δ + Λ + δ ) γ e + i ( β i δ ) + Λ Λ ( Λ + γ + δ ) v 1
We conjecture that their result holds true, since in this case also the third order polynomial resulting from the elimination when R 0 > 1 satisfies
p ( 0 ) = ( R 0 1 ) γ e Λ ( Λ + γ + δ ) v 1 < 0 , p ( i c ) = β i γ γ e 2 Λ δ > 0 ,
and so p ( i ) must still have either one, two or three roots in ( 0 , i c ) .

6. Conclusions and Further Work

Our paper highlighted several open problems for SIR-PH, SIR-PH-IA, and SIR-PH-FA models. The following general directions seem worthy of further work.
1.
Investigate whether the beautiful determinant Formula (17) hidden in the papers of [11,12] (valid when δ = 0 ) may be extended to SIR-PH models, exploiting the partition (26).
2.
Study the case of two or more compartments susceptible to become infected (for example the SEIT model [13,57]).
3.
Study the scaled model as a perturbation of the FA model.
4.
Study stability via the geometric approach of Li, Graef, Wang, Karsai, Muldoney, and Lu [15].
5.
We hope that the use of more sophisticated and fast software will allow researchers in the future to progress with the interesting questions raised by models with higher dimensions. Here, exploiting symmetries may turn out helpful.

Author Contributions

Supervision, Florin Avram; Writing-review and editing, Florin Avram, Rim Adenane and Andrei Halanay. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

We thank Mattia Sensi and Sara Sottile for useful suggestions and references and for providing some of the codes used for producing the figures.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. 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. [Google Scholar] [CrossRef] [PubMed]
  2. 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. [Google Scholar] [CrossRef] [PubMed]
  3. Vyska, M.; Gilligan, C. Complex dynamical behaviour in an epidemic model with control. Bull. Math. Biol. 2016, 78, 2212–2227. [Google Scholar] [CrossRef] [Green Version]
  4. Roostaei, A.; Barzegar, H.; Ghanbarnejad, F. Emergence of Hopf bifurcation in an extended SIR dynamic. arXiv 2021, arXiv:2107.10583. [Google Scholar]
  5. Gupta, R.; Kumar, A. Endemic bubble and multiple cusps generated by saturated treatment of an SIR model through Hopf and Bogdanov–Takens bifurcations. Math. Comput. Simul. 2022, 197, 1–21. [Google Scholar] [CrossRef]
  6. 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. [Google Scholar] [CrossRef] [Green Version]
  7. De la Sen, M.; Ibeas, A.; Alonso-Quesada, S.; Nistal, R. On the Carrying and Evolution Matrices in Epidemic Models. J. Phys. Conf. Ser. 2021, 1746, 012015. [Google Scholar] [CrossRef]
  8. 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. [Google Scholar]
  9. Avram, F.; Adenane, R.; Ketcheson, D.I. A review of matrix SIR Arino epidemic models. Mathematics 2021, 9, 1513. [Google Scholar] [CrossRef]
  10. Avram, F.; Adenane, R.; Basnarkov, L.; Bianchin, G.; Goreac, D.; Halanay, A. On matrix-SIR Arino models with linear birth rate, loss of immunity, disease and vaccination fatalities, and their approximations. arXiv 2021, arXiv:2112.03436. [Google Scholar]
  11. Robinson, M.; Stilianakis, N.I. A model for the emergence of drug resistance in the presence of asymptomatic infections. Math. Biosci. 2013, 243, 163–177. [Google Scholar] [CrossRef] [PubMed]
  12. Ottaviano, S.; Sensi, M.; Sottile, S. Global stability of SAIRS epidemic models. Nonlinear Anal. Real World Appl. 2022, 65, 103501. [Google Scholar] [CrossRef]
  13. Van den Driessche, P.; Watmough, J. Further notes on the basic reproduction number. In Mathematical Epidemiology; Springer: Berlin/Heidelberg, Germany, 2008; pp. 159–178. [Google Scholar]
  14. Greenhalgh, D. Hopf bifurcation in epidemic models with a latent period and nonpermanent immunity. Math. Comput. Model. 1997, 25, 85–107. [Google Scholar] [CrossRef]
  15. 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. [Google Scholar] [CrossRef] [Green Version]
  16. Li, M.Y.; Muldowney, J.S. Dynamics of differential equations on invariant manifolds. J. Differ. Equ. 2000, 168, 295–320. [Google Scholar] [CrossRef] [Green Version]
  17. Sun, C.; Hsieh, Y.H. Global analysis of an SEIR model with varying population size and vaccination. Appl. Math. Model. 2010, 34, 2685–2697. [Google Scholar] [CrossRef]
  18. Britton, T.; Ouédraogo, D. SEIRS epidemics with disease fatalities in growing populations. Math. Biosci. 2018, 296, 45–59. [Google Scholar] [CrossRef] [Green Version]
  19. Lu, G.; Lu, Z. Global asymptotic stability for the SEIRS models with varying total population size. Math. Biosci. 2018, 296, 17–25. [Google Scholar] [CrossRef]
  20. Douris, P.; Markakis, M. Global Connecting Orbits of a SEIRS Epidemic Model with Nonlinear Incidence Rate and Nonpermanent Immunity. Eng. Lett. 2019, 27, 1–10. [Google Scholar]
  21. Ansumali, S.; Kaushal, S.; Kumar, A.; Prakash, M.K.; Vidyasagar, M. Modelling a pandemic with asymptomatic patients, impact of lockdown and herd immunity, with applications to SARS-CoV-2. Annu. Rev. Control 2020, 50, 432–447. [Google Scholar] [CrossRef]
  22. Graef, J.R.; Li, M.Y.; Wang, L. A Study on the Effects of Disease Caused Death in a Simple Epidemic Model; American Institute of Mathematical Sciences: Springfield, MI, USA, 1998; Volume 1998, p. 288. [Google Scholar]
  23. Ferreira, R.A.; Silva, C.M. A nonautonomous epidemic model on time scales. J. Differ. Equ. Appl. 2018, 24, 1295–1317. [Google Scholar] [CrossRef]
  24. Russo, G.; Wirth, F. Matrix measures, stability and contraction theory for dynamical systems on time scales. arXiv 2020, arXiv:2007.08879. [Google Scholar] [CrossRef]
  25. 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. [Google Scholar] [CrossRef]
  26. Goh, B.S. Global stability in many-species systems. Am. Nat. 1977, 111, 135–143. [Google Scholar] [CrossRef]
  27. Bacaër, N. Mathématiques et Épidémies. 2021. Available online: https://hal.archives-ouvertes.fr/hal-03331469 (accessed on 20 March 2022).
  28. Ballyk, M.M.; McCluskey, C.C.; Wolkowicz, G.S. Global analysis of competition for perfectly substitutable resources with linear response. J. Math. Biol. 2005, 51, 458–490. [Google Scholar] [CrossRef]
  29. Fall, A.; Iggidr, A.; Sallet, G.; Tewa, J.J. Epidemiological models and Lyapunov functions. Math. Model. Nat. Phenom. 2007, 2, 62–83. [Google Scholar] [CrossRef]
  30. Li, M.Y.; Muldowney, J.S. Global stability for the SEIR model in epidemiology. Math. Biosci. 1995, 125, 155–164. [Google Scholar] [CrossRef]
  31. 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. [Google Scholar] [CrossRef] [Green Version]
  32. Li, M.Y.; Muldowney, J.S. A geometric approach to global-stability problems. SIAM J. Math. Anal. 1996, 27, 1070–1083. [Google Scholar] [CrossRef]
  33. Shuai, Z.; van den Driessche, P. Global stability of infectious disease models using Lyapunov functions. SIAM J. Appl. Math. 2013, 73, 1513–1532. [Google Scholar] [CrossRef] [Green Version]
  34. LeVeque, R.J.; Mitchell, I.M.; Stodden, V. Reproducible research for scientific computing: Tools and strategies for changing the culture. Comput. Sci. Eng. 2012, 14, 13–17. [Google Scholar] [CrossRef]
  35. Wang, W. Backward bifurcation of an epidemic model with treatment. Math. Biosci. 2006, 201, 58–71. [Google Scholar] [CrossRef] [PubMed]
  36. 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. [Google Scholar] [CrossRef] [Green Version]
  37. Zhou, L.; Fan, M. Dynamics of an SIR epidemic model with limited medical resources revisited. Nonlinear Anal. Real World Appl. 2012, 13, 312–324. [Google Scholar] [CrossRef]
  38. Sun, S. Global dynamics of a SEIR model with a varying total population size and vaccination. Int. J. Math. Anal. 2012, 6, 1985–1995. [Google Scholar]
  39. Wiggers, S.L.; Pedersen, P. Routh–hurwitz-liénard–chipart criteria. In Structural Stability and Vibration; Springer: Berlin/Heidelberg, Germany, 2018; pp. 133–140. [Google Scholar]
  40. Anderson, B.; Jury, E. A simplified Schur-Cohn test. IEEE Trans. Autom. Control 1973, 18, 157–163. [Google Scholar] [CrossRef]
  41. Daud, A.A.M. A Note on Lienard-Chipart Criteria and its Application to Epidemic Models. Math. Stat. 2021, 9, 41–45. [Google Scholar] [CrossRef]
  42. Ma, J.; Earn, D.J. Generality of the final size formula for an epidemic of a newly invading infectious disease. Bull. Math. Biol. 2006, 68, 679–702. [Google Scholar] [CrossRef]
  43. Feng, Z. Final and peak epidemic sizes for SEIR models with quarantine and isolation. Math. Biosci. Eng. 2007, 4, 675. [Google Scholar]
  44. Andreasen, V. The final size of an epidemic and its relation to the basic reproduction number. Bull. Math. Biol. 2011, 73, 2305–2321. [Google Scholar] [CrossRef]
  45. Georgescu, P.; Hsieh, Y.H. Global stability for a virus dynamics model with nonlinear incidence of infection and removal. SIAM J. Appl. Math. 2007, 67, 337–353. [Google Scholar] [CrossRef] [Green Version]
  46. Hu, Z.; Ma, W.; Ruan, S. Analysis of SIR epidemic models with nonlinear incidence rate and treatment. Math. Biosci. 2012, 238, 12–20. [Google Scholar] [CrossRef] [PubMed]
  47. Razvan, M. Multiple equilibria for an SIRS epidemiological system. arXiv 2001, arXiv:math/0101051. [Google Scholar]
  48. 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. [Google Scholar] [CrossRef] [Green Version]
  49. Riaño, G. Epidemic Models with Random Infectious Period. medRxiv 2020. [Google Scholar] [CrossRef]
  50. Plemmons, R.J. M-matrix characterizations. I—Nonsingular M-matrices. Linear Algebra Its Appl. 1977, 18, 175–188. [Google Scholar] [CrossRef] [Green Version]
  51. Hurtado, P.J.; Kirosingh, 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. [Google Scholar] [CrossRef] [Green Version]
  52. 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. [Google Scholar] [CrossRef]
  53. Horn, R.A.; Johnson, C.R. Matrix Analysis; Cambridge University Press: Cambridge, UK, 2012. [Google Scholar]
  54. 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. [Google Scholar] [CrossRef]
  55. Diekmann, O.; Heesterbeek, J.; Roberts, M.G. The construction of next-generation matrices for compartmental epidemic models. J. R. Soc. Interface 2010, 7, 873–885. [Google Scholar] [CrossRef] [Green Version]
  56. Gasca, M.; Pena, J.M. Total positivity and Neville elimination. Linear Algebra Its Appl. 1992, 165, 25–44. [Google Scholar] [CrossRef] [Green Version]
  57. Martcheva, M. An Introduction to Mathematical Epidemiology; Springer: Berlin/Heidelberg, Germany, 2015; Volume 61. [Google Scholar]
Figure 1. 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 R 0 = 3.21 , β = 5 , γ = 1 / 2 , Λ = μ = 1 / 10 , γ r = 1 / 6 , γ s = 0.01 , δ i = 0.9 , δ r = 0 . E E S c , E E I n , E E F O A are the stable endemic points of the scaled model, intermediate model, and the FA model, respectively. The green vertical line denotes the immunity threshold 1 / R = s E E F O A = s E E I n . 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).
Figure 1. 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 R 0 = 3.21 , β = 5 , γ = 1 / 2 , Λ = μ = 1 / 10 , γ r = 1 / 6 , γ s = 0.01 , δ i = 0.9 , δ r = 0 . E E S c , E E I n , E E F O A are the stable endemic points of the scaled model, intermediate model, and the FA model, respectively. The green vertical line denotes the immunity threshold 1 / R = s E E F O A = s E E I n . 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).
Symmetry 14 00995 g001
Figure 2. The existence and uniqueness of endemic point i e e in the interval [ 0 , i c ] for the [17] problem, ([15], Figure 1), with f S H ( i ) as defined in (46), f S H ( Λ δ ) = ( γ + δ ) γ e β i Λ + δ γ s Λ δ ( Λ + γ + δ ) v 1 Λ + γ s , β = 18 , β e = γ r = 0 , Λ = 4 5 and δ = γ = e i = γ s = γ e = 1 .
Figure 2. The existence and uniqueness of endemic point i e e in the interval [ 0 , i c ] for the [17] problem, ([15], Figure 1), with f S H ( i ) as defined in (46), f S H ( Λ δ ) = ( γ + δ ) γ e β i Λ + δ γ s Λ δ ( Λ + γ + δ ) v 1 Λ + γ s , β = 18 , β e = γ r = 0 , Λ = 4 5 and δ = γ = e i = γ s = γ e = 1 .
Symmetry 14 00995 g002
Figure 3. Bifurcation Diagrams with respect to β . (a) Forward bifurcation diagram with respect to β , DFE stable from 0 to β * S H where the coordinates of the endemic point are negative and the endemic point is stable after the critical bifurcation parameter β * S H , with β * S H : = ( Λ + γ + δ ) v 1 Λ + γ s Λ γ e , and δ = 1 . Bifurcation Diagrams with respect to δ . Here, γ r = 0 , Λ = 4 5 and γ = γ e = e i = γ s = 1 . (b) Bifurcation Diagram with respect to δ , the endemic point is stable from 0 to δ * : = Λ β i γ e ( Λ + γ e ) ( Λ + γ s ) γ and the DFE is stable when δ δ * , with β = 18 . The four points before the yellow one ensure the uniqueness of the endemic point when R 0 > 1 , see ([17], Theorem 3.3).
Figure 3. Bifurcation Diagrams with respect to β . (a) Forward bifurcation diagram with respect to β , DFE stable from 0 to β * S H where the coordinates of the endemic point are negative and the endemic point is stable after the critical bifurcation parameter β * S H , with β * S H : = ( Λ + γ + δ ) v 1 Λ + γ s Λ γ e , and δ = 1 . Bifurcation Diagrams with respect to δ . Here, γ r = 0 , Λ = 4 5 and γ = γ e = e i = γ s = 1 . (b) Bifurcation Diagram with respect to δ , the endemic point is stable from 0 to δ * : = Λ β i γ e ( Λ + γ e ) ( Λ + γ s ) γ and the DFE is stable when δ δ * , with β = 18 . The four points before the yellow one ensure the uniqueness of the endemic point when R 0 > 1 , see ([17], Theorem 3.3).
Symmetry 14 00995 g003aSymmetry 14 00995 g003b
Figure 4. Parametric Plot of (s,e+i) of the SEIRS scaled model [17] “with slightly modified parameters” and its FA and intermediate approximations starting from a starting point SP= ( 0.0828333 , 0.0015 ) , with R 0 = 3.57143 , β = 57.1429 , γ r = 0 , δ = 2 , γ = γ e = γ s = Λ = 1 . EESc, EEIn, EEFA denote the intermediate points of the scaled model, interdemiate approximation, and FA approximation, respectively. The green vertical line denotes the immunity threshold 1 / R = s E E F A = s E E I n .
Figure 4. Parametric Plot of (s,e+i) of the SEIRS scaled model [17] “with slightly modified parameters” and its FA and intermediate approximations starting from a starting point SP= ( 0.0828333 , 0.0015 ) , with R 0 = 3.57143 , β = 57.1429 , γ r = 0 , δ = 2 , γ = γ e = γ s = Λ = 1 . EESc, EEIn, EEFA denote the intermediate points of the scaled model, interdemiate approximation, and FA approximation, respectively. The green vertical line denotes the immunity threshold 1 / R = s E E F A = s E E I n .
Symmetry 14 00995 g004
Figure 5. The existence and uniqueness of the endemic point i e e for the [19] problem “with slightly modified parameters", with i c = 7 12 and β i = 9 , Λ = 1 6 , δ = 2 , γ = 1 2 , γ e = 2 , γ r = 1 , e i = 2 .
Figure 5. The existence and uniqueness of the endemic point i e e for the [19] problem “with slightly modified parameters", with i c = 7 12 and β i = 9 , Λ = 1 6 , δ = 2 , γ = 1 2 , γ e = 2 , γ r = 1 , e i = 2 .
Symmetry 14 00995 g005
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Avram, F.; Adenane, R.; Halanay, A. New Results and Open Questions for SIR-PH Epidemic Models with Linear Birth Rate, Loss of Immunity, Vaccination, and Disease and Vaccination Fatalities. Symmetry 2022, 14, 995. https://doi.org/10.3390/sym14050995

AMA Style

Avram F, Adenane R, Halanay A. New Results and Open Questions for SIR-PH Epidemic Models with Linear Birth Rate, Loss of Immunity, Vaccination, and Disease and Vaccination Fatalities. Symmetry. 2022; 14(5):995. https://doi.org/10.3390/sym14050995

Chicago/Turabian Style

Avram, Florin, Rim Adenane, and Andrei Halanay. 2022. "New Results and Open Questions for SIR-PH Epidemic Models with Linear Birth Rate, Loss of Immunity, Vaccination, and Disease and Vaccination Fatalities" Symmetry 14, no. 5: 995. https://doi.org/10.3390/sym14050995

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop