SEIR

SEIR()

A Susceptible-Exposed-Infectious-Recovered (SEIR) model is a more realistic extension of the standard SIR model in which individuals do not become instantly infectious upon exposure, but undergo an incubation period, the timescale of which is governed by the parameter, \(\alpha\):

\[\begin{split}\begin{aligned} \frac{\mathrm{d}S}{\mathrm{d}t} &= - \frac{\beta SI}{N} \\ \frac{\mathrm{d}E}{\mathrm{d}t} &= \frac{\beta SI}{N} - \alpha E \\ \frac{\mathrm{d}I}{\mathrm{d}t} &= \alpha E - \gamma I \\ \frac{\mathrm{d}R}{\mathrm{d}t} &= \gamma I \end{aligned}\end{split}\]

We use the flu-like parameters of the SIR model demonstration with an incubation period of 2 days.

from pygom import common_models
import matplotlib.pyplot as plt
import numpy as np

# Parameters
n_pop = 1e4
gamma = 1/4
alpha = 1/2
R0 = 1.3
beta = R0*gamma

model = common_models.SEIR(
    {
        'beta':beta,
        'gamma':gamma,
        'alpha':alpha,
        'N':n_pop
    }
)

# Output times
tmax = 365
dt = 1
t = np.arange(0, tmax, dt)

# Initial conditions
i0 = 1
x0 = [n_pop-i0, 0, i0, 0]
model.initial_values = (x0, t[0])

# Deterministic evolution
solution = model.solve_deterministic(t)

We also run an SIR model with the same parameters to compare the outputs

model = common_models.SIR(
    {
        'beta':beta,
        'gamma':gamma,
        'N':n_pop
    }
)

x0 = [n_pop-i0, i0, 0]
model.initial_values = (x0, t[0])

solution_sir = model.solve_deterministic(t)

We see that the SEIR model changes the profile of the epidemic as compared with an SIR model, but the overall final sizes are the same.

state_names = ["S", "E", "I", "R"]
sir_idx = {"S": 0, "I": 1, "R": 2}

fig, axes = plt.subplots(1, 4, layout="constrained", figsize=(10, 3))

for i, state_name in enumerate(state_names):
    axes[i].plot(t, solution[0].result.y[:, i], label="SEIR")

    if state_name in sir_idx:
        axes[i].plot(
            t,
            solution_sir[0].result.y[:, sir_idx[state_name]],
            label="SIR",
        )

    axes[i].set_title(state_name)
    axes[i].set_xlabel("Time")
axes[0].legend()

plt.show()