The previous model I built was a toy model: the simplest possible dynamics that would give you a feel for what’s going on, and parameter selection that was educated guesswork.
There was sufficient interest in it, though, that I decided it would be worth building a more sophisticated and hopefully accurate model.
State dynamics and parameterization
The previous model was a simple SIR model with only three states. Individuals were either susceptible, meaning vulnerable to infection but not yet infected, infected as well as infectious to susceptible individuals, or recovered (or dead). This is very simple to calculate, but has some obvious limitations. It doesn’t show the progression of the disease in individuals and the delay from action. Critical cases, for example, were simply calculated as a static fraction of all infected states at any point, when in fact they should be time delayed (it takes time for an infection to progress to a critical state). It’s hard to use to examine the effects of specific controls; all I did was vary the effective 
For the new model, I’ve spent some time reading case reports and aggregated studies. Key references are:
- A clinical analysis of 1099 cases throughout China describing clinical course (Feb 28)
- A case study of 9 clustered German cases with daily testing using both PCR and viral culturing (Mar 8)
- A survey of hospital associations estimating total hospital bed and ICU bed capacity for the U.S. (Mar 17)
- A survey/modeling study in two Chinese cities counting contact numbers, including between different age groups (Mar 20)
- A CDC report of hospitalization rates, ICU admission rates, and mortality rates of confirmed American cases (Mar 26)
- A CDC report of testing results from two cruise ships, including one where all crew and passengers were tested (Mar 26)
- A detailed report of Chinese cases in Wuhan with estimation of transmission dynamics (Mar 26)
Based on my reading, it looks like the progression of the disease is as follows:
- An incubation period of 4.1-7.0 days after infection, during which the person is not yet infectious to others.
- A short pre-symptomatic period of ~2 days when the person is infectious but before onset of symptoms.
- Some fraction of people are fully asymptomatic, and can spread the disease until it is cleared and they recover, about 12 days after the end of the incubation period.
- The rest will develop noticeable upper respiratory symptoms, for a duration of 12 days before recovery or decline.
- Some fraction of symptomatic cases will progress to non-severe or severe support needs from a hospital (for non-severe cases) or ICU (for severe cases). Hospitalization lasts a further 11-13 days, depending on severity.
- Some fraction of hospitalized cases will die; the rest will recover and retain immunity to the disease.
Here is a diagram of the disease course:
Of particular note is that I have modeled hospitalized patients as no longer being infectious; this is due to the German study suggesting low viral shedding 10 days post symptom onset.
The basic
The asymptomatic fraction is estimated at 50% based on data from the Diamond Princess cruise ship. All passengers and crew on the ship were tested and ~50% of the positive test results had shown no symptoms.
The results of symptomatic cases are age-stratified based on the CDC report of American cases, including rates of hospitalization, ICU admission, and death. New infections are assumed to be randomly distributed among age brackets proportional to their fraction of the susceptible population.
The most speculative part of the model is the effects of health care saturation. While I have reasonable estimates of the available hospital and ICU beds, it’s necessarily speculative to estimate what happens to people with cases requiring hospitalization if that hospitalization is not available. I’ve made what I feel are likely conservative guesses. I’ve assumed triage of scarce beds to younger and less-severe patients.
Model with no controls
The first test is simply to run the model without attempting to impose infection controls and see what the result is. Here is a log-scale plot of total hospitalizations, ICU admissions, and deaths over the 2-year course of the epidemic.
The deaths top out at ~10 million total. Note that this is higher than my previous model, in part because of better estimation of available beds and in part because of longer estimated hospital stays. The death rate is strongly biased towards the elderly, in part due to intrinsically worse outcomes and in part due to triage of scarce resources.
Nearly 20% of 85+ Americans would be expected to die, and even 1 in 19 of 55-64 year olds, 1 in 47 of 45-54 year olds, and 1 in 164 20-44 year olds.
Without the impact of scarce hospital resources, the results look much better:
You can see that the red line no longer bends upward when the hospitalized cases hit the ICU or hospital bed capacity, reducing the total deaths almost 4-fold to 2.8 million. As before, the deaths are heavily skewed towards the elderly:
Unfortunately, avoiding running out of hospital capacity at the peak would require millions of available ICU beds simultaneously, which is rather unrealistic. We’ll look into the effects of imposed infection controls in my next post.
Appendix: code
PopState = namedtuple("PopState", [
"susceptible",
"incubating",
"presymptomatic",
"asymptomatic",
"symptomatic",
"nonsevere",
"severe",
"dead",
"recovered",
])
GeneralParams = namedtuple("GeneralParams", [
"incubating_days",
"presympt_days",
"asympt_days",
"sympt_days",
"nonsev_days",
"severe_days",
"R0",
"frac_asymptomatic",
])
AgeBracketParams = namedtuple("AgeBracketParams", [
"frac_nonsevere",
"frac_severe",
"nonsevere_dead_home",
"severe_dead_icu",
"severe_dead_hosp",
"severe_dead_home",
])
shared_params = GeneralParams(
5.0, # Days incubating before infectious
2.0, # Days infectious before symptomatic
10.0, # Days infectious for asymptomatic cases
12.5, # Days infectious for symptomatic cases
11.0, # Days hospitalized for non-severe cases
13.0, # Days hospitalized for severe cases
2.2, # R0, base infectivity rate
0.5, # Estimated fraction asymptomatic cases
)
# These are high-end estimates from the CDC publication
bracket_params = {
"age_0_19": AgeBracketParams(
0.025, 0.0,
0.10,
0.0, 0.0, 0.0,
),
"age_20_44": AgeBracketParams(
0.166, 0.042,
0.15,
0.05, 0.10, 0.15,
),
"age_45_54": AgeBracketParams(
0.179, 0.104,
0.18,
0.08, 0.16, 0.24,
),
"age_55_64": AgeBracketParams(
0.189, 0.112,
0.20,
0.23, 0.45, 0.85,
),
"age_65_74": AgeBracketParams(
0.247, 0.188,
0.25,
0.26, 0.55, 0.95,
),
"age_75_84": AgeBracketParams(
0.277, 0.310,
0.31,
0.34, 0.67, 1.0,
),
"age_85": AgeBracketParams(
0.413, 0.290,
0.35,
0.94, 1.0, 1.0,
),
}
bracket_states = {
"age_0_19": PopState(81.98e6, 0, 0, 0, 0, 0, 0, 0, 0),
"age_20_44": PopState(108.84e6, 10, 0, 0, 0, 0, 0, 0, 0),
"age_45_54": PopState(41.63e6, 0, 0, 0, 0, 0, 0, 0, 0),
"age_55_64": PopState(42.27e6, 0, 0, 0, 0, 0, 0, 0, 0),
"age_65_74": PopState(30.48e6, 0, 0, 0, 0, 0, 0, 0, 0),
"age_75_84": PopState(15.39e6, 0, 0, 0, 0, 0, 0, 0, 0),
"age_85": PopState(6.55e6, 0, 0, 0, 0, 0, 0, 0, 0),
}
brackets = [
"age_0_19",
"age_20_44",
"age_45_54",
"age_55_64",
"age_65_74",
"age_75_84",
"age_85",
]
ICU_BEDS = 31482
HOSPITAL_BEDS = 501470
def update_state(states, general, params, day=1, start_day=70, control_rate=2.2, icu_beds=ICU_BEDS, hospital_beds=HOSPITAL_BEDS):
def total_state(state):
return sum(states[age].__getattribute__(state) for age in brackets)
total_S = total_state("susceptible")
total_N = sum(sum(states[age]) for age in brackets)
# Generate the new states as dicts so they're mutable
new_states = {age: states[age]._asdict() for age in brackets}
# Adjust the controlled rate
R0 = control_rate if day >= start_day else general.R0
# Adjust for constant time between contacts
mean_infectious_days = general.presympt_days
mean_infectious_days += general.asympt_days * general.frac_asymptomatic
mean_infectious_days += general.sympt_days * (1 - general.frac_asymptomatic)
days_betw_contacts = mean_infectious_days / R0
# susceptible --> incubating
total_preI = total_state("presymptomatic")
new_infections_from_presympt = (total_preI * total_S) / (total_N * days_betw_contacts)
total_asymI = total_state("asymptomatic")
new_infections_from_asympt = (total_asymI * total_S) / (total_N * days_betw_contacts)
total_symI = total_state("symptomatic")
new_infections_from_sympt = (total_symI * total_S) / (total_N * days_betw_contacts)
new_infections = new_infections_from_presympt + new_infections_from_asympt + new_infections_from_sympt
for age in brackets:
# each age gets new infections proportionate to its fraction of the susceptible population
# (i.e., panmixia across all age groups)
transfer = (states[age].susceptible / total_S) * new_infections
assert transfer <= states[age].susceptible
assert transfer >= 0
new_states[age]["susceptible"] -= transfer
new_states[age]["incubating"] += transfer
# incubating --> presymptomatic
for age in brackets:
transfer = (states[age].incubating / general.incubating_days)
assert transfer <= states[age].incubating
assert transfer >= 0
new_states[age]["incubating"] -= transfer
new_states[age]["presymptomatic"] += transfer
# presymptomatic --> asymptomatic / symptomatic
for age in brackets:
transfer = (states[age].presymptomatic / general.presympt_days)
assert transfer <= states[age].presymptomatic
assert transfer >= 0
new_states[age]["presymptomatic"] -= transfer
new_states[age]["asymptomatic"] += transfer * general.frac_asymptomatic
new_states[age]["symptomatic"] += transfer * (1 - general.frac_asymptomatic)
# asymptomatic --> recovered
for age in brackets:
transfer = (states[age].asymptomatic / general.asympt_days)
assert transfer <= states[age].asymptomatic
assert transfer >= 0
new_states[age]["asymptomatic"] -= transfer
new_states[age]["recovered"] += transfer
# symptomatic --> recovered / nonsevere / severe
for age in brackets:
transfer = (states[age].symptomatic / general.sympt_days)
assert transfer <= states[age].symptomatic
assert transfer >= 0
new_states[age]["symptomatic"] -= transfer
new_states[age]["severe"] += transfer * params[age].frac_severe
new_states[age]["nonsevere"] += transfer * params[age].frac_nonsevere
new_states[age]["recovered"] += transfer * (1 - (params[age].frac_severe + params[age].frac_nonsevere))
# nonsevere --> recovered / dead
# and
# severe --> recovered / dead
free_hospital_beds = hospital_beds
free_icu_beds = icu_beds
for age in brackets:
# nonsevere patients get first pick of hospital beds for triage
transfer = (states[age].nonsevere / general.nonsev_days)
assert transfer <= states[age].nonsevere
assert transfer >= 0
if transfer > 0:
new_states[age]["nonsevere"] -= transfer
# youngest patients get first pick of hospital beds for triage
nonsevere = states[age].nonsevere
in_bed = min(nonsevere, free_hospital_beds)
at_home = nonsevere - in_bed
free_hospital_beds -= in_bed
dying = (at_home / nonsevere) * transfer * params[age].nonsevere_dead_home
assert dying <= transfer
new_states[age]["dead"] += dying
new_states[age]["recovered"] += transfer - dying
# severe patients get what's left over
transfer = (states[age].severe / general.severe_days)
assert transfer <= states[age].severe
assert transfer >= 0
if transfer > 0:
new_states[age]["severe"] -= transfer
# bed shortage
severe = states[age].severe
in_icu = min(severe, free_icu_beds)
free_icu_beds -= in_icu
in_bed = min(severe - in_icu, free_hospital_beds)
free_hospital_beds -= in_bed
at_home = severe - (in_icu + in_bed)
dying = transfer * (
(at_home / severe) * params[age].severe_dead_home +
(in_bed / severe) * params[age].severe_dead_hosp +
(in_icu / severe) * params[age].severe_dead_icu
)
assert dying <= transfer
new_states[age]["dead"] += dying
new_states[age]["recovered"] += transfer - dying
# Now cast everything back to states
final_states = {
age: PopState(
new_states[age]["susceptible"],
new_states[age]["incubating"],
new_states[age]["presymptomatic"],
new_states[age]["asymptomatic"],
new_states[age]["symptomatic"],
new_states[age]["nonsevere"],
new_states[age]["severe"],
new_states[age]["dead"],
new_states[age]["recovered"],
)
for age in brackets
}
return final_states, day + 1
def run_model(states, general, params):
def total_state(state):
return sum(states[age].__getattribute__(state) for age in brackets)
def total_infections():
return sum([
total_state("incubating"),
total_state("presymptomatic"),
total_state("asymptomatic"),
total_state("symptomatic"),
total_state("nonsevere"),
total_state("severe"),
])
infections = [total_infections()]
nonsevere = [total_state("nonsevere")]
severe = [total_state("severe")]
deaths = [total_state("dead")]
day = 1
while total_infections() >= 0.25:
states, day = update_state(states, general, params, day=day, start_day=230, control_rate=2.2,
icu_beds=ICU_BEDS, hospital_beds=HOSPITAL_BEDS
)
infections.append(total_infections())
nonsevere.append(total_state("nonsevere"))
severe.append(total_state("severe"))
deaths.append(total_state("dead"))
return np.array(infections), np.array(nonsevere), np.array(severe), np.array(deaths), states
I, H, C, D, states = run_model(bracket_states, shared_params, bracket_params)
plt.figure(figsize=(12, 8))
plt.grid()
plt.plot(H, label="All Hospitalizations", color='k')
plt.plot(C, label="ICU admissions", color='orange')
plt.plot(D, label="Deaths", color='red')
plt.legend()
plt.yscale('log')
plt.ylim(ymin=1)
plt.xlim((0, 600))
plt.hlines(ICU_BEDS, 0, 600, linestyle='dashed', color='g')
plt.hlines(HOSPITAL_BEDS + ICU_BEDS, 0, 600, linestyle='dashed', color='b')
plt.ylabel("# of People")
plt.xlabel("Days After Initial Infection")
death_rates = [100 * states[age].dead / (states[age].dead + states[age].recovered + states[age].susceptible) for age in brackets]
plt.figure(figsize=(8, 6))
plt.bar(range(len(brackets)), death_rates)
plt.grid()
labels = [
"0-19",
"20-44",
"45-54",
"55-64",
"65-74",
"74-85",
"85+",
]
plt.xticks(range(len(brackets)), labels, rotation=45)
plt.xlabel("Age Bracket")
plt.ylabel("% of Bracket Dying of COVID-19")
Arguably Wrong 
You say: “The asymptomatic fraction is estimated at 50% based on data from the Diamond Princess cruise ship. All passengers and crew on the ship were tested and ~50% of the positive test results had shown no symptoms.” This is based on one test per person, which may have been conducted before symptoms appeared. 50% is roughly the ceiling.
Another paper estimates that on the Diamond Princess the true asymptomatic rate (ie, those who never developed symptoms) was 15-20%. If they accounted for false negatives that would drop a bit, but probably not much.
https://www.eurosurveillance.org/content/10.2807/1560-7917.ES.2020.25.10.2000180
The guy leading the South Korean response said 20-25% of their cases are asymptomatic, but I haven’t seen the relevant numbers presented in the literature.
Of course, all of these estimates apply to adults, not to children.
LikeLike
Yup, quite possible. All these parameters are pretty poorly pinned down right now. If the asymptomatic rate is lower than 50%, we’ll see higher incidence of hospitalization/ICU use/death. I’d bet large sums the asymptomatic rate isn’t ~99%, the way that Oxford model was supposing.
LikeLike
Have you considered formalizing this in a system of differential equations published in an online stock and flow system like insightmaker.com?
LikeLike
The referenced “German study” study only looks at patients with a mild disease state. They are a special cohort known as the ‘Munich cluster’ which was epidemiologically linked and closely observed based on contact tracing resulting from a specific event. So I don’t think it supports the case that all hospital patients are generally not shedding 10 days after symptom onset.
Here’s a quote from the lead author, Wendtner, “In terms of scientific significance, our study benefited from the fact that all of the cases were linked to an index case, meaning they were not simply studied based on the presence of certain symptoms. In addition to getting a good picture of how this virus behaves, this also enabled us to gain other important insights, including on viral transmission.”
This article goes into more detail about this special group of patients.
https://www.reuters.com/article/us-health-coronavirus-germany-defences-i/pass-the-salt-the-minute-details-that-helped-germany-build-virus-defences-idUSKCN21R1DB
LikeLike
It’s quite possible that the mild symptoms of these patients with low shedding are qualitatively different from those of sicker hospitalized patients, but I just don’t know of any other studies where they actually tried to culture the virus to see if they’re shedding active virus. If you can point me to any data on viral shedding from other groups, I’d be very grateful.
LikeLike
Here’s a paper where they don’t verify that the detected RNA was infectious, but it does indicate that certain patients shed for different lengths of time. If sicker patients are shedding for longer than asymptomatic patients then I don’t think we necessarily have reason to believe those particles are not infectious.
Risk Factors for Viral RNA Shedding in COVID-19 Patients
https://erj.ersjournals.com/content/early/2020/05/07/13993003.01190-2020
“The median time for viral RNA shedding in patients with or without CHD were 21 days and 19 days”
LikeLike