A review of matrix SIR Arino epidemic models

Many of the models used nowadays in mathematical epidemiology, in particular in COVID-19 research, belong to a certain sub-class of compartmental models whose classes may be divided into three"(x, y, z)"groups, which we will call respectively"susceptible/entrance, diseased, and output"(in the classic SIR case, there is only one class of each type). Roughly, the ODE dynamics of these models contain only linear terms, with the exception of products between x and y terms. It has long been noticed that the basic reproduction number R has a very simple formula (3.3) in terms of the matrices which define the model, and an explicit first integral formula (3.8) is also available. These results can be traced back at least to [ABvdD+07] and [Fen07], respectively, and may be viewed as the"basic laws of SIR-type epidemics"; however many papers continue to reprove them in particular instances (by the next-generation matrix method or by direct computations, which are unnecessary). This motivated us to redraw the attention to these basic laws and provide a self-contained reference of related formulas for (x, y, z) models. We propose to rebaptize the class to which they apply as matrix SIR epidemic models, abbreviated as SYR, to emphasize the similarity to the classic SIR case. For the case of one susceptible class, we propose to use the name SIR-PH, due to a simple probabilistic interpretation as SIR models where the exponential infection time has been replaced by a PH-type distribution. We note that to each SIR-PH model, one may associate a scalar quantity Y(t) which satisfies"classic SIR relations", see (3.8). In the case of several susceptible classes, this generalizes to (5.10); in a future paper, we will show that (3.8), (5.10) may be used to obtain approximate control policies which compare well with the optimal control of the original model.

Finally, Section 5 reviews briefly the case of several classes of susceptibles. This topic requires further development; we include it however due to the recognized importance of heterogeneity factors.

The classic Kermack-McKendrick SIR epidemic model
The SIR process (S(t), I(t), R(t), t ≥ 0) divides a population of size N undergoing an epidemic into three classes called "susceptibles, infectives and removed". One may also work with the corresponding fractions s(t) = S(t) N , i(t) = I(t) N , and r(t) = 1− s(t)− i(t). It is assumed that only susceptible individuals can get infected. After having been infectious for some time, an individual recovers and may not become susceptible again. "Viewed from far away", this yields the SIR model with demography [KM27b,BCCF19] S (t) = − β N S(t)I(t) + ξ (N − S(t)) , where 1. N is the total, constant population size.
2. R (t), the number of removed per unit time, is the only quantity which is clearly observable, at least in the easy case when the removed are dead, as was the case of the original study of the Bombay plague [KM27b].
3. ξ is the population death rate, assumed to equal the birth rate.
4. γ is the removal rate of the infectious, which equals 1/duration of the infection (under the stochastic model of exponential infection durations, this is the reciprocal of the expected duration).
5. β, the infection rate, models the probability that a contact takes place between an infected and a susceptible, and that it results in infection.
Note that 1. The sum S + I + R = N is conserved and each value is positive, so the values of S, I, R remain in the interval [0, N ].
2. This system has a unique solution, since (given the boundedness of S, I, and R), the RHS above is Lipschitz.
From now on, we will assume that ξ = 0, and work with the fractions s, i, r, which satisfy We will call this the classic SIR model. Note that 1. s(t) is monotonically decreasing and r(t) is monotonically increasing, to, say, s ∞ , r ∞ ; therefore convergence to some fixed stable point (s ∞ , i ∞ , r ∞ ) must hold.

solutions starting in the domain
4. The second equation of (2.1) implies the so-called threshold phenomenon: if then i(t) decreases always, without any intervention.
To avoid trivialities, we will assume R > 1 from now on.
5. When R > 1, the epidemic grows iff s > 1/R, i.e. until the susceptibles s(t) reach the immunity threshold after which infections decline. R is called basic reproduction number, and it models the number of susceptibles infected by one infectious (expected number, under more sophisticated stochastic, branching models).
An advantage of the classic SIR model is that it is essentially solvable explicitly: 1. We can eliminate r from the system using the invariant s + i + r = 1 (this is also possible for various generalizations like SIR with demography, as long as r does not appear explicitly in the rest of the equations).
2. It can easily be verified that is invariant, so that i is explicitly given by (2.7) 3. The maximal value of the infected i, achieved when s = 1/R, is (2.8) 4. By differentiating the right-hand side of (2.7), one finds that the maximal value of the "newly infected" − s = β i s is achieved when where the Lambert function L −1 is a real inverse of L(z) = ze z -see for example [Pak15,KS20,BS20]. Bounding s i is one interesting possibility for accomodating ICU constraints [Man20, (2.20)].
5. The infectious class converges to 0 and the susceptible and recovered converge monotonically to limits which may be expressed in terms of the "Lambert-W(right)" function L 0 .
Let us note that accurate numerical solutions of the evolution of the SIR or other compartmental epidemic may be obtained very quickly. 3 SIR-PH epidemics with one susceptible class (SIR epidemics with phase-type "disease time")) It has been known for a long while that R and the final size for many compartmental model epidemics may be explicitly expressed in terms of the matrices which define the model, and [ME06, ABvdD + 07, Fen07, And11] offer a quite general framework of "xyz" models which ensures this. We believe that these formulas have not received the attention they deserve (they keep being reproved), and decided therefore to review them below; we will call them matrix-SIR models.
A particular but revealing case is that when there is only one susceptible class, which we will call SIR-PH, following Riano [Ria20], who emphasized its probabilistic interpretation -see also [HK19].
Definition 3.1 A "SIR-PH ( α, V, β, W ) epidemic" contains a homogeneous susceptible class, but vector "diseased" state i (which may contain latent/exposed, infective , asymptomatic, etc) and vector removed states (healthy, dead, vaccinated, etc). It is defined by an ODE system: where 1. i (t) ∈ R n is a row vector whose components are fractions of diseased individuals of various types, which must satisfy i (∞) = 0.
2. β ∈ R n is a column vector whose components represent the relative transmission ability of the various disease classes.
3. α ∈ R n is a probability row vector with the components representing the fractions of susceptibles entering into the corresponding disease compartments, when infection occurs.

4.
A is a n × n Markovian sub-generator matrix describing rates of transition between the diseased classes i (i.e., a Markovian generator matrix for which the sum of at least one row is strictly negative). Alternatively, V := −A is a non-singular M-matrix. § 5. r(t) ∈ R p is a row vector which must satisfy r(∞) > 0, whose components represent (fractions of ) various classes which survive at the end of an infection.
6. W ∈ R n × R p is a matrix whose components represent the rates at which classes of diseased individuals become recovered. We assume that the matrix V = A W ∈ R m+n × R n has row sums 0, which implies mass conservation.
We turn now to a deceivingly simple particular example of the SIR-PH model, which explains its name.
Remark 3.2 Probabilistic interpretation of SIR-PH epidemics. For simplicity, let us group all the output classes of (3.1) into one r = r1, yielding: where we put a := W 1 = (−A)1.
(3.2) emphasizes the fact that SIR-PH models are in one to one correspondence with laws of phase-type ( α, A) [Ria20, (21)].
Let us recall now, as known essentially since [Kur78] -see also [?, Thm. 2.2.7]-that under proper scaling, the expected fractions s(t), i(t), r(t) of stochastic SIR § and more general compartmental models obey a "law of large numbers/fluid limit" which recovers the deterministic epidemic.
As an example, the SIR-PH model (3.1) may be derived as limit of a stochastic SIR model in which the exponential infection time has been replaced by a phase-type ( α, A) "dwell period" [HK19].
To illustrate the power of the SIR-PH formalism, consider now the case with two diseased states, latent and infectious, with phase-type dwell times, parametrized by ( α e , A e ) and ( α i , A i ), respectively. Using the well-known convolution formula -see for example [BN17, Thm. 3.1.26] we find that formulas like (3.3) (see other examples of such formulas below) still apply, with ( α, A, β) given by (3.4) § One such model stipulates that each infective j = 1, ..., J infects a randomly chosen susceptible, at encounter times which belong to independent Poisson processes P j (t), j = 1, ..., J, of rate β, and that infection durations are i.i.d. r.v.'s which are exponentially distributed at the end of which the individual recovers (or dies).
§ This can be also derived as the expected number of susceptibles infected during a dwell period, for the stochastic model (the so-called "survival method")-see [Per18] for an excellent review The "epidemic dwell strucure" ( α, A, β) of examples with more complicated network structures for the diseased may be constructed using Kronecker sums of the matrices defining each component.
Let us give now an example which does not in general belong to the SIR-PH class.
Example 1 The SIRV model (SIR with vaccination -see for example [BBGS20]) is defined by: This is of the form (3.1) with i = ( i), r = ( r, v) iff γ s = 0.
In the opposite case γ s = 0, one may still compute an invariant and for fixed s, putting γ = γ γs , β = β γs , µ 0 = µ 0 γs , it holds that i is explicitly given by When γ s > 0, the final size is s ∞ = 0.
We provide now a list of several formulas, obtained by replacing i in SIR by a scalar linear combination (3.8) [Fen07]. They are all easily proved; however the formula for the maximal value of the newly infected involves also a second linear combination y (3.12).
The solution of Z(s) = Z(0) with respect to s may be expressed in terms of is the principal branch of the Lambert-W function.
2. The derivative with respect to time is by the conservation of Y (t) + s(t) − 1 R ln( s(t)) between the time 0 and the time t 1/R of reaching the immunity threshold).
6. The final size of the removed satisfies: 7. The value of the infected combination Y when s = 1/R is (3.20) 8. The maximum size of the newly infected is achieved when Remark 3.6 Let us note that for control problems involving optimization objectives which only depend on Y (t), we are effectively optimizing a SIR model; this SIR approximation may be used to offer practical solutions for optimizing more complicated compartmental models.
Example 2 For SEIR, putting i = ( e, i), we may write

Examples of SIR-PH models used in COVID-19 modelling
We derive now R and Y from (3.3), (3.8), for some popular compartmental models. Note that we will be reformulating the original results (which, unfortunately, have already appeared several times with different notations), using a unifying notation.
Example 3 The SEIHRD model [Ivo17, PZL20, PF20, NHS20, RFVP + 21] has disease states i = ( e, i, h). We use here the version in [PF20] (we would rather call this SI 2 HRD model), defined by where we denoted by γ e , γ i , γ h the sum of the constant rates out of e, i, h, and by i h the rate out of i and reaching h, etc. Then, R = βe γe + e i γe β i γ i [PF20, 2], and Y = e+ i γeβ i βeγ i +e i β i . When β e = 0 = e r =⇒ e i γe = 1, we recover R = β i γ i [PZL20] and Y = e + i. Example 4 The SEIHCRD model of [KK20] has disease states i = ( e, i, h, c) and is defined by Then, Example 6 The SI aps QR model (with asymptomatic, pre-symptomatic and symptomatic infectious) [SK21,(3.2)], [GRH21].
The disease states are i = ( i a , i p , i s , q), and the model is defined by Then, [SK21], and, where R a := βa γa , R p := βp γp , and R s := βs γs .

S (m) Y R models with m groups of susceptibles
The SIR-compartment model makes the unrealistic assumption that the population through which the disease is spreading is well-mixed. However, differences in susceptibility and rates of contact between individuals strongly affect their likelihood of catching COVID-19. A model which attempts to capture this aspect is: Lemma 5.1 A disease free equilibrium (s 1 , s 2 , . . . , s m , 0, r 0 ) of (5.1) is asymptotically stable iff sR < 1, where s = k s k and is the spectral radius of the next generation matrix.
While the final size may also be obtained under this model [And11, Thm. 2.1], for transient behavior it is convenient to turn to a simpler model.

A generalization of heterogeneous SEIR [DT20]
Assume now that β k = β k β, where β k ∈ R + and W=β w, where w is a row vector. Putting y = i β, the system (5.1) becomes: and r(t) =ˆt 0 y(u)du w := I(t) w.
We conclude with some preliminary results on this model.
Lemma 5.2 a) [DT20] Put s(t) = k s k (t), p k = s k (0)/ s(0).The solution of (5.3) satisfies the time dependent SYR system: where a(t) = k β k s k (t) s(t) (5.9) is a positive non-increasing function with a(0) = k p k β k :=β. b) Y (t) defined in (3.8) with R = k p k R k , R k = β k αV −1 β, satisfies: and is unimodal, with a maximum on the immunity/recovery line k s k R k = 1. Proof. a) The derivative of a satisfies