key: cord-0195417-2sn04nu3
authors: Avram, Florin; Adenane, Rim; Ketcheson, David I.
title: A review of matrix SIR Arino epidemic models
date: 2021-06-01
journal: nan
DOI: nan
sha: 18e0b34105785da7ffdfb76d0bdb4c6c46e71a9e
doc_id: 195417
cord_uid: 2sn04nu3

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 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.

2. the equilibrium set of stable points is (s, 0, 1 − s), s ∈ [0, 1].

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.

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). 

has the property that

and that z(t) = µ( s(t), Y (t)) := Y (t) + s(t) − 1 R ln[ s(t)], Z(t) = e −Rz(t) = s(t)e −R( s(t)+Y (t)) (3.10) are constant along the paths of the dynamical system (3.1).

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).

4. The final size of the susceptibles satisfies[ABvdD + 07, Thm.5.1]:

(3.14)

by the conservation of Y (t) + s(t) − 1 R ln( s(t)) between the times 0 and ∞; explicitly,

5. The integrated infectives I (a,b) =´b a i (s)ds satisfies

, (3.16) and the total integrated infectives I (∞) =´∞ 0 i (s)ds satisfies [ABvdD + 07, (6)] Note that this has been used to model the total cost of an epidemic [GJ72] .

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

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 .

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.

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.

(5.4) § Such a dynamics was first considered in [Gar68] .

It is convenient to reparametrize the model taking I as parameter, or, equivalently, by taking r = r1 = γI, γ := w1.

(5.5)

we find that the system has a family of first integrals which includes 1 R 1 log(s 1 /s 1 (0)) = ... = 1 R k log(s k /s k (0)) = ... where γ is defined in (5.5).

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 

A simple planning problem for COVID-19 lockdown

A final size relation for epidemic models

The final size of an epidemic and its relation to the basic reproduction number

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

Hamiltonian structure of compartmental epidemiological models

On the formulation of epidemic models (an appraisal of kermack and mckendrick)

Matrix-exponential distributions in applied probability

Exact and approximate analytic solutions in the SIR epidemic model

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

Micea T Sofonea, and Samuel Alizon, Optimal COVID-19 epidemic control until vaccine deployment, medRxiv

Optimal timing of oneshot interventions for epidemic control

A data driven analysis and forecast of an SEIARD epidemic model for COVID-19 in Mexico

Heterogeneous social interactions and the COVID-19 lockdown outcome in a multi-group SEIR model

An extended epidemic model on interconnected networks for COVID-19 to explore the epidemic dynamics

A light introduction to modelling recurrent epidemics, Mathematical epidemiology

Final and peak epidemic sizes for SEIR models with quarantine and isolation

A feedback SIR (fsir) model highlights advantages and limitations of infection-based social distancing

Optimal control of the transmission rate in compartmental epidemics

The mathematical analysis of an epidemic with two kinds of susceptibles

The cost of a general stochastic epidemic

Modeling the evolution of sars-cov-2 under non-pharmaceutical interventions

Data-driven control of the COVID-19 outbreak via non-pharmaceutical interventions: A geometric programming approach

Generalizations of the 'linear chain trick': incorporating more flexible dwell time distributions into mean field ode models

Balancing quarantine and self-distancing measures in adaptive epidemic networks

Stability analysis of a compartmental seihrd model for the ebola virus disease

Arif Ahmed Seikh, and Debashree Guha, A fuzzy dynamic optimal model for COVID-19 epidemic in india based on granular differentiability

Optimal control of an SIR epidemic through finite-time non-pharmaceutical intervention

Beyond just "flattening the curve": Optimal control of epidemics with purely non-pharmaceutical interventions

A contribution to the mathematical theory of epidemics

A contribution to the mathematical theory of epidemics

Early dynamics of transmission and control of COVID-19: a mathematical modelling study

Analytical solution of the SIRmodel for the temporal evolution of epidemics. part a: Time-independent reproduction factor

Strong approximation theorems for density dependent Markov chains

A divide and conquer strategy against the COVID-19 pandemic?!, medRxiv

An introduction to mathematical epidemiology

Generality of the final size formula for an epidemic of a newly invading infectious disease

Uziel Shemesh, ϑ-seihrd mathematical model of covid19-stability analysis using fast-slow decomposition

Estimating the basic reproductive number of COVID-19 cases in ghana

Lambert's W meets Kermack-McKendrick epidemics

An introduction to the basic reproduction number in mathematical epidemiology

A control approach to the covid-19 disease using a seihrd dynamical model, medRxiv

M-matrix characterizations. i-nonsingular mmatrices

Population modeling of early COVID-19 epidemic dynamics in french regions and estimation of the lockdown impact on infection rate

Optimal control of COVID-19 infection rate with social costs

A simple but complex enough ϑ-sir type model to be used with covid-19 real data. application to the case of italy

Epidemic models with random infectious period, medRxiv

On COVID-19 modelling

What the reproductive number r0 can and cannot tell us about covid-19 dynamics

Yannis Michalakis, and Samuel Alizon, Epidemiological monitoring and control perspectives: application of a parsimonious modelling framework to the COVID-19 dynamics in france

Evolving epidemiology and impact of non-pharmaceutical interventions on the outbreak of coronavirus disease