Epidemiological modeling

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)

Join the Conversation

12 Comments

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

    Like

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

    Like

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

        Like

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

    Like

Leave a comment

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

Create your website at WordPress.com
Get started
%d bloggers like this: