A basic SIR model for the epidemic.
Total population 327 million, with a single initiating infection. 12 day course of disease, initial R0 of 3.5 (as per the stats from China that Steve linked.) Assume 5% of cases are critical, 2% of those critical cases die with ICU care, 5% with general hospital care, and 50% with no care. 90,000 available ICU beds and 900,000 available hospital beds.
Run this model through to completion and it sweeps through the population in about a year, infecting most everyone (less about 9 million who escape), killing 5.4 million.
Now, suppose we impose infection controls on day 80, right about when there’s 1000 deaths from this thing. And then we vary how strong those controls are: from 0.35 (what the Chinese managed) up to nothing at all.
Here we see how the # of deaths varies with the strength of our controls. If we impose Chinese-style controls, we get away with only 5K deaths, or 1000-fold fewer than without the controls. But the inflection point isn’t just at 1.0. In particular if we can get the R0 below about 1.5 that gets us down under 500K, or a 10-fold reduction. At 1.0 we get down to about 50K.
If we take the slope of this plot, we can see the marginal impact of better controls at each point, in terms of # of deaths saved:
Why the discontinuity at R0 of just below 1.5? Well, here’s the number of critical cases over time for several different scenarios. The vertical line is when we impose controls, the horizontal lines are ICU bed capacity and total hospital bed capacity. Right below an R0 of 1.5, we dip the # of critical cases low enough that we never exceed hospital capacity. The trick here though is that until you get R0 below 1.0 you’re still sweeping through most of the population, just more slowly, slowly enough that the hospitals can handle it. And that means that, for example, if you get R0 to exactly 1.0 you have to keep those infection controls in place for ~5 years or else you’ll be back up close to 5million dead again.
Appendix: code to generate this model.
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from collections import namedtuple
PopState = namedtuple("PopState", [
"day",
"susceptible",
"infected",
"dead",
"recovered",
])
EpiParams = namedtuple("EpiParams", [
"days_infected",
"infection_rate",
"frac_critical",
"frac_dead_icu",
"frac_dead_hospital",
"frac_dead_home",
"start_day",
"controlled_infection_rate",
])
ICU_BEDS = 90000
HOSPITAL_BEDS = 900000
def update_state(state, params):
day = state.day
S = state.susceptible
I = state.infected
D = state.dead
R = state.recovered
N = S + I + D + R
days_betw_contacts = params.days_infected / params.infection_rate
if day >= params.start_day:
days_betw_contacts = params.days_infected / params.controlled_infection_rate
new_infections = (I * S) / (N * days_betw_contacts)
recovered_infections = (I / params.days_infected) * (1 - params.frac_critical)
# For critical cases, we split into ICU, hospital, and home care, in descending order
C = I * params.frac_critical;
C_icu = min(C, ICU_BEDS);
C_hospital = 0 if C < ICU_BEDS else min(C - ICU_BEDS, HOSPITAL_BEDS)
C_home = 0 if C < (ICU_BEDS + HOSPITAL_BEDS) else C - (ICU_BEDS + HOSPITAL_BEDS)
recovered_critical = (
((C_icu / params.days_infected) * (1 - params.frac_dead_icu)) +
((C_hospital / params.days_infected) * (1 - params.frac_dead_hospital)) +
((C_home / params.days_infected) * (1 - params.frac_dead_home))
)
dead_critical = (
((C_icu / params.days_infected) * (params.frac_dead_icu)) +
((C_hospital / params.days_infected) * (params.frac_dead_hospital)) +
((C_home / params.days_infected) * (params.frac_dead_home))
)
newS = S - new_infections
newI = I + new_infections - recovered_infections - recovered_critical - dead_critical
newR = R + recovered_infections + recovered_critical
newD = D + dead_critical
return PopState(day + 1, newS, newI, newD, newR)
def run_model(state, params):
infections = [state.infected]
deaths = [state.dead]
while (state.infected) >= 1:
state = update_state(state, params)
infections.append(state.infected)
deaths.append(state.dead)
return infections, deaths, state
total_population = 327 * 1e6
init_state = PopState(
day=1,
susceptible=total_population - 1,
infected=1,
dead=0,
recovered=0)
params = EpiParams(
days_infected=12,
infection_rate=3.5,
frac_critical=0.05,
frac_dead_icu=0.02,
frac_dead_hospital=0.05,
frac_dead_home=0.5,
start_day=79,
# Modify this parameter to vary the strength of infection controls
controlled_infection_rate=3.5,
)
infections_over_time, deaths_over_time, final_state = run_model(init_state, params)
Arguably Wrong 
in the last line of code, the argument ‘test_param’ is not defined.
LikeLike
Thanks. Fixed the error.
LikeLike
Thanks for sharing.
LikeLike
Does the number of ICU beds = 90000 exclude those that are normally used within the US?
That is, are they the total ICU beds available for corona-virus cases after all the other needs for ICU beds have been handled, or will they simply displace other needs for ICU beds?
LikeLike
90000 assumes about a 33% improvement in current capacity and displacement of existing ICU needs. Same for the ~900000 hospital beds.
LikeLike
So it sounds like we would need an additional 120,000 ICU beds.
LikeLike
Hi!
Based on what currently happens in Italy, I believe the frac_dead parameters are different than in this model. For frac_dead_icu we are probably close to 0.6 as opposed to 0.02. frac_dead_hospital is more difficult: It is not automatically higher than icu (probably they are in less serious conditions than icu patients). The same goes for frac_dead_home.
The Penn model suggests an R0 at 2.71 but again, based on what is happening in the last weeks, your 3.5 seems more correct.
Thank you for sharing. BTW, I am “translating” your python source to php.
LikeLike
Ah, and another input:
Yesterday it was reported that, on average, infected people shows symptoms for 4 days before entering ICU where they stay for another 4 days before decease. That probably also modifies the 12 day course parameter.
LikeLike
In fact, I think we must introduce more parameters.
days_infected should be split into two parts:
1. days_infected_spreader. This parameter deals with the period the infected person is spreading the virus. It affects new infections.
2. days_infected_quarantained. This parameter tells how many days the person is infected without the possibility to affect anybody else. It affects the death and recovery numbers. The assumption of zero spreading is theoretical since there is poof of hospital personnel and people visiting the home quarantained person being infected.
LikeLike
You might take a look at the updated model, it’s likely more tweakable to what you’re looking for: https://arguablywrong.home.blog/2020/03/27/updated-epidemiological-modeling/
LikeLike
The one thing this model doesn’t really take into account is that in the perfect petri dish that we have so far – the Diamond Princess cruise ship – the infection rate was less than 20%.
It overstates the case dramatically.
LikeLike