Thread by @gro_tsen

I've written a blog post about generalizing the SIR model to the case of variable (=heterogeneous) susceptibility. It's in French, so let me translate the key points in a Twitter thread because I think this is important and deserves more attention. •1/25 madore.org/~david/weblog/…

Just in case, I mention that Google Translate seems to do a fairly decent job on my blog post (though it does mess up the formula formatting): translate.google.com/translate?sl=f… — also, here's a Twitter thread in French which I'll mostly be translating below: •2/25

https://twitter.com/gro_tsen/status/1373628002786025472

SIR is the most basic model in epidemiology, and is still much used by some (though with variations not relevant here). It's the one which predicts things like a herd immunity threshold at 1 − 1/R₀ (n.b.: in what follows, I'll write “κ” for R₀). •3/25

This SIR model makes a lot of simplifying hypotheses, among these the fact that the population is homogeneous, so everyone has the same susceptibility to the infection (the same probability of being infected at a given time). •4/25

What I wish to describe is what happens if susceptibility is variable in the population (some people are more susceptible than others): we assume known its (initial) distribution s₀ and try to model the effect on epidemic evolution. •5/25

But let me emphasize first that here I'm talking about variations in susceptibility, NOT infectiousness (→some people infect others more easily, e.g., “superspreaders”). The latter does not have, per se, any influence on the epidemic as a whole, it just averages out. •6/25

The intuitive idea, which I've explained many times before, is that if some people are more susceptible, they will be infected — so made immune — earlier on, so as immunity accumulates, not only are there fewer susceptible people, they are also less susceptible on average. •7/25

So accumulation of immunity is made more efficient by this heterogeneity in susceptibility. And our goal is to model this in a mathematically precise way, and quantify the phenomenon. What's surprising is how well this works! •8/25

(All of this is essentially due to work by @mgmgomes1 and her coauthors, and probably all contained in the paper arxiv.org/abs/2008.00098 — I reworked it and reformulated it myself, but the ideas, as far as I know, are hers.) •9/25

Classical SIR in one tweet is this system of diff eqns:
‣ ds/dt = −β·i·s
‣ di/dt = β·i·s − γ·i
‣ dr/dt = γ·i
‣ (s+i+r=1)
with t being time, s,i,r proportions of susceptible, infected and recovered, β contagiousness, γ recovery rate and reproduction number κ := β/γ.
•10/25

It predicts in particular: an initial exponential growth with rate (=log slope) β, a herd immunity threshold at 1 − 1/κ, and a final attack rate of 1 + W(−κ·exp(−κ))/κ ≈ 1 − exp(−κ), where W is the Lambert transcendental function. •11/25

Moving to a heterogeneous susceptibility is remarkably easy: everything can be computed from the (initial) distribution s₀ of susceptibility through its Laplace transform φ (meaning φ(u) := 𝔼(exp(−u·X)) where 𝔼:=“expected value”, … •12/25

… and X is a random variable distributed following the initial susceptibility profile s₀ normalized by 𝔼(X)=1, viꝫ φ′(0)=−1). In fact, in SIR we just replace β·i·s (new infections term) by −β·i·φ′(φ⁻¹(s)) to get the heterogeneous model! •13/25

In other words, the equations for heterogeneous SIR are as follows:
‣ ds/dt = β·i·φ′(φ⁻¹(s))
‣ di/dt = − β·i·φ′(φ⁻¹(s)) − γ·i
‣ dr/dt = γ·i
‣ (s+i+r=1)
where φ is the (known) Laplace transform, φ⁻¹ its inverse function and φ′ its derivative.
•14/25

The classical SIR case, where the susceptibility profile s₀ is a δ distribution at 1, corresponds to φ(u) = exp(−u). Another important case is exponential distribution s₀(x) = exp(−x), giving φ(u) = 1/(1+u). I'll discuss the Γ distribution, generalizing this, below. •15/25

In fact, the epidemic evolution can be read from φ: it partially solves the differential equations (namely, we constantly have s = φ(κ·r), where κ := β/γ) and makes it possible to compute the herd immunity threshold and final attack rate, as I now explain. •16/25

The herd immunity threshold is the point where graph of φ has slope −1/κ (recall that κ := β/γ is the reproduction number), the final attack rate is the intersection of this graph with the line (also of slope −1/κ) through (0,1) and (κ,0) having abscissa >0. •17/25

(In the previous tweet, by “is the point where”, I mean that the point has coordinates (κ·r, s), and 1−s is the sought-after quantity. The original system can be described as the evolution of a point on the graph of φ with these coordinates.) •18/25

One important special case of susceptibility distribution is the Γ distribution, having relative variance v, say (inverse of the “shape” parameter). The exponential distribution is a special case of this (namely v=1), and the homogeneous case is the v→0 limit. •19/25

For this Γ distribution, we have φ(u) = (1+u·v)^(−1/v) and the equations are obtained by replacing β·i·s in classical SIR by β·i·s^(1+v), namely:
‣ ds/dt = −β·i·s^(1+v)
‣ di/dt = β·i·s^(1+v) − γ·i
‣ dr/dt = γ·i
‣ (s+i+r=1)
•20/25

What consequences does this have?
✱ Initially the epidemic still follows an exponential curve with rate β, as in the homogeneous case.
✱ But as a (small) amount of immunity accumulates, it is now 1+v times more effective at reducing the effective reproduction number, …
•21/25

… because when a proportion s<1 of population is susceptible, average susceptibility is now reduced to s^v (rel. to what it was initially).
✱ Herd immunity threshold, final attack rate and peak infection fraction are reduced w.r.t classical SIR (see graphs below).
•22/25

☞ Note: failing experimental data, I see no particular reason to believe that v=0 should be a better epidemic model than another value. Dimensional analysis suggests we assume v~1 (exponential) a priori. If, say, children are almost immune, clearly v=0 cannot hold! •23/25

Here are curves illustrating the dynamic for basic reproduction number κ=3 with v=0 left (=classical SIR =homogeneous) and v=1 right (exponential distribution). S is green, I is red, R is blue and effective reproduction number is black (bottom-right graph). •24/25

The herd immunity threshold is 1 − κ^(−1/(1+v)) (red below) and final attack rate is the >0 solution to (κ·v·r+1)^(−1/v) + r = 1 (blue). Here are their graphs in function of v for reproduction numbers κ=2 (solid), 3 (dashed), 4 (dotted). •25/25