OK, the tweet below was using a very simplistic heuristic; now let me write a thread about what the (hardly less simplistic!) SIR epidemiological model says about the attack rate of an epidemic. ⤵️ •1/23
The SIR model is a first-order nonlinear ordinary differential equation in two variables written s,i,r. (Yes, that makes three variables, but s+i+r=1 so one is redundant.) They represent proportions of the population that is susceptible (s), infected (i) and recovered (r) •2/23
(here, “susceptible” means “not yet infected”). If I write x′ for the time derivative of x, the SIR equations are s′ = −β·i·s and i′ = β·i·s − γ·i (so r′ = γ·i), for two constants β and γ describing the epidemic (and having dimensions of inverse time): … •3/23
… here, β·i·s represents new infections, and γ·i represents recoveries (either neglecting deaths or counting them as “recoveries”). I.e., we model recoveries as a first order kinetic process (essentially, the infected population, left alone, recovers at a rate γ), … •4/23
… and infections as a kinetic process of first order in each of i and s (so, second order overall kinetic, where “order” here means degree, it's still a first-order ODE, terminology sucks 🙄). Initially, r=0, i is very small but nonzero, and s=1−i is very close to 1. •5/23
Of particular importance is the ratio β/γ, which is the “basic reproduction number” κ (or R₀) of the infection. Note that this is the only free parameter in the equation, because we can rescale time to multiply β and γ by a common constant. •6/23
This “basic reproduction number” is usually denoted R₀, which is absolutely awful notation 😒 when doing a model called SIR where R stands for Recovered. So I'm writing r for “recovered” and I'll write κ := β/γ for this basic reproduction number. •7/23
Initially, of course, s remains close to 1, so i grows exponentially like exp(βt). This is what we're currently seeing with #Covid19. We can measure β as the logarithmic derivative of number of cases and γ from individual patient recovery (1/γ = expected time), … •8/23
… so we can compute the basic reproduction number κ = β/γ which is essentially the number of people each infected person infects in turn. Current estimates for #Covid19 are around β ~ 0.2/d, 1/γ ~ 15d, so κ ~ 3. Very roughly. Anyway. •9/23
Typical behavior is shown in the following graphs taken from Wikimedia Commons (commons.wikimedia.org/wiki/File:Sirs…), where s is in blue, i in green and r in red, and s+i+r here is 500 because some people don't know how to divide to normalize to 1, but no matter. •10/23 
What interests me is what s and r converge to as t→∞. One might think from the above graphs (and I did) that s→0, i.e., the epidemic stops when there are no susceptible individuals left, but this isn't true! Actually, s approaches a nonzero limit which I'll call s[∞]. •11/23
In other words, not everyone gets infected, only a proportion r[∞] = 1−s[∞] do (the “final attack rate”). (Of course, i[∞]=0 as is obvious from the γ·i term in the equations.) Now, how do we compute s[∞]? Well, from s′ = −β·i·s and r′ = γ·i we get, by dividing: … •12/23
… the relation s′/r′ = −κ·s; but s′/r′ is the derivative of s with respect to r, so this can be solved as s = s₀·exp(−κ·r). In particular, s[∞] = s₀·exp(−κ·r[∞]), but since also r[∞] = 1−s[∞], we have the equation s[∞] = s₀·exp(−κ·(1−s[∞])). •13/23
Assuming we start the epidemic with an infinitesimal number of cases, s₀ is essentially 1, so we're trying to solve the equation s[∞] = exp(−κ·(1−s[∞])) in the indeterminate s[∞] (proportion never infected) with parameter κ (reproduction number). •14/23
For κ<1, the solution is 1 (i.e., the epidemic gets nowhere). So let's take κ>1. The equation s = exp(−κ·(1−s)) is a transcendental equation, it doesn't have a solution in closed form using standard functions but it has one using a special function W: … •15/23
… namely s = −W(−κ·exp(−κ))/κ where W is Lambert's transcendental W function (=product logarithm), see en.wikipedia.org/wiki/Lambert_W… for more about this function but basically w := W(z) satisfies w·exp(w) = z. •16/23
(This solution comes from rewriting s = exp(−κ·(1−s)) as −κ·s·exp(−κ·s) = −κ·exp(−κ), hence −κ·s = W(−κ·exp(−κ)). Of course you need to check that we have the right branches of W, but I don't want to get into details.) •17/23
To summarize, the SIR model predicts that the proportion of the population that remains uninfected to the end is s = −W(−κ·exp(−κ))/κ where W is Lambert's function and κ>1 is the basic reproduction number. How does this function behave? •18/23
And how does it compare with the estimate s = 1/κ given by the simplistic reasoning given in the thread quoted in the first tweet of this thread? Well, in short, it's worse. 😕 (Essentially because SIR tracks the epidemic passed its “turning point” κ·s~1.) •19/23
For κ = 1 + h a basic reproduction number ever so slightly above 1, we have s = −W(−κ·exp(−κ))/κ = 1 − 2·h + (8/3)·h² + O(h³): compared to 1/κ which is 1 − h + h² + O(h³), we're basically predicting twice as many cases. •20/23
And for κ large, using the expansion W(x) = x − x² + O(x³) we find s = −W(−κ·exp(−κ))/κ = exp(−κ) + κ·exp(−2κ) + O(κ²·exp(−3κ)), so the fraction s of uninfected drops exponentially with the basic reproduction number κ. •21/23
For κ~3, SIR predicts s~6% remain uninfected, i.e., r~94% attack rate. I'm a mathematician, not an epidemiologist, so lest anyone claim I wrote #Covid19 will infect 94% of the Earth's population, let me emphasize that this is a highly simplistic model showing its limits! •22/23
Here's a plot showing the prediction s = −W(−κ·exp(−κ))/κ of uninfected in the SIR model, in blue, compared to s = 1/κ predicted by the other reasoning I mentioned, in red, in function of the basic reproduction number κ in abscissa. •23/23
✱ Correction: in tweet 8/23, I should have written “i grows exponentially like exp((β−γ)t) (so long as s~1)” since i′=(β−γ)i, not like exp(βt). My estimates for β should be correspondingly updated, but they are just orders of magnitude anyway. •24/(23+1)
