This notebook contains material from CBE40455-2020; content is available on Github.
2.2 Campus SEIR Modeling¶
2.2.1 Campus infection data¶
The following data consists of new infections reported since August 3, 2020, from diagnostic testing administered by the Wellness Center and University Health Services at the University of Notre Dame. The data is publically available on the Notre Dame Covid-19 Dashboard.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from scipy.integrate import solve_ivp
from scipy.optimize import minimize
from datetime import timedelta
data = [
["2020-08-03", 0],
["2020-08-04", 0],
["2020-08-05", 0],
["2020-08-06", 1],
["2020-08-07", 0],
["2020-08-08", 1],
["2020-08-09", 2],
["2020-08-10", 4],
["2020-08-11", 4],
["2020-08-12", 7],
["2020-08-13", 10],
["2020-08-14", 14],
["2020-08-15", 3],
["2020-08-16", 15],
["2020-08-17", 80],
]
df = pd.DataFrame(data, columns=["date", "new cases"])
df["date"] = pd.to_datetime(df["date"])
fig, ax = plt.subplots(figsize=(8,4))
ax.bar(df["date"], df["new cases"], width=0.6)
ax.xaxis.set_major_locator(mdates.WeekdayLocator(byweekday=mdates.MO))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %d"))
plt.title("Reported New Infections")
plt.grid()
2.2.2 Fitting an SEIR model to campus data¶
Because of the limited amount of data available at the time this notebook was prepared, the model fitting has been limited to an SEIR model for infectious disease in a homogeneous population. In an SEIR model, the progression of an epidemic can be modeled by the rate processes shown in the following diagram.
$$\text{Susceptible} \xrightarrow {\frac{\beta S I}{N}} \text{Exposed} \xrightarrow{\alpha E} \text{Infectious} \xrightarrow{\gamma I} \text{Recovered} $$which yeild the following model for the population of the four compartments
$$\begin{align*} \frac{dS}{dt} &= -\beta S \frac{I}{N} \\ \frac{dE}{dt} &= \beta S \frac{I}{N} - \alpha E \\ \frac{dI}{dt} &= \alpha E - \gamma I \\ \frac{dR}{dt} &= \gamma I \\ \end{align*}$$The recovery rate is given by $\gamma = 1/t_{recovery}$ where the average recovery time $t_{recovery}$ is estimated as 8 days.
| Parameter | Description | Estimated Value | Source |
|---|---|---|---|
| $N$ | campus population | 15,000 | estimate |
| $\alpha$ | 1/average latency period | 1/(3.0 d) | |
| $\gamma$ | 1/average recovery period | 1/(8.0 d) | literature |
| $\beta$ | infection rate constant | tbd | fitted to data |
| $I_0$ | initial infectives on Aug 3, 2020 | tbd | fitted to data |
| $R_0$ | reproduction number | ${\beta}/{\gamma}$ |
N = 15000 # estimated campus population
gamma = 1/8.0 # recovery rate = 1 / average recovery time in days
alpha = 1/3.0
def model(t, y, beta):
S, E, I, R = y
dSdt = -beta*S*I/N
dEdt = beta*S*I/N - alpha*E
dIdt = alpha*E - gamma*I
dRdt = gamma*I
return np.array([dSdt, dEdt, dIdt, dRdt])
def solve_model(t, params):
beta, I_initial = params
IC = [N - I_initial, I_initial, 0.0, 0.0]
soln = solve_ivp(lambda t, y: model(t, y, beta), np.array([t[0], t[-1]]),
IC, t_eval=t, atol=1e-6, rtol=1e-9)
S, E, I, R = soln.y
U = beta*S*I/N
return S, E, I, R, U
def residuals(df, params):
S, E, I, R, U = solve_model(df.index, params)
return np.linalg.norm(df["new cases"] - U)
def fit_model(df, params_est=[0.5, 0.5]):
return minimize(lambda params: residuals(df, params), params_est, method="Nelder-Mead").x
def plot_data(df):
plt.plot(df.index, np.array(df["new cases"]), "r.", ms=20, label="data")
plt.xlabel("days")
plt.title("new cases")
plt.legend()
def plot_model(t, params):
print("R0 =", round(beta/gamma, 1))
S, E, I, R, U = solve_model(t, params)
plt.plot(t, U, lw=3, label="model")
plt.xlabel("days")
plt.title("new cases")
plt.legend()
plot_data(df)
beta, I_initial = fit_model(df)
plot_model(df.index, [beta, I_initial])
2.2.3 Fitted parameter values¶
from tabulate import tabulate
parameter_table = [
["N", 15000],
["I0", I_initial],
["beta", beta],
["gamma", gamma],
["R0", beta/gamma]
]
print(tabulate(parameter_table, headers=["Parameter", "Value"]))
2.2.4 Short term predictions of newly confirmed cases¶
Using the fitted parameters, the following code presents a short term projection of newly diagnosed infections. Roughly speaking, the model projects a 50% increase per day in newly diagnosed cases as a result of testing sympotomatic individuals.
The number of infected but asympotomatic individuals is unknown at this time, but can be expected to be a 2x multiple of this projection.
# prediction horizon (days ahead)
H = 1
# retrospective lag
K = 6
fig, ax = plt.subplots(1, 1, figsize=(12, 4))
for k in range(0, K+1):
# use data up to k days ago
if k > 0:
beta, I_initial = fit_model(df[:-k])
P = max(df[:-k].index) + H
c = 'b'
a = 0.25
else:
beta, I_initial = fit_model(df)
P = max(df.index) + H
c = 'r'
a = 1.0
# simulation
t = np.linspace(0, P, P+1)
S, E, I, R, U = solve_model(t, [beta, I_initial])
# plotting
dates = [df["date"][0] + timedelta(days=t) for t in t]
ax.plot(dates, U, c, lw=3, alpha=a)
ax.plot(df["date"], df["new cases"], "r.", ms=25, label="new infections (data)")
ax.xaxis.set_major_locator(mdates.WeekdayLocator(byweekday=mdates.MO))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %d"))
ax.grid(True)
ax.set_title(f"{H} day-ahead predictions of confirmed new cases");
