Programming a DTMC SEIR model

SEIR separates infection from infectiousness. A newly infected person first enters an exposed compartment and becomes infectious only after a latent period.

The new biological idea: exposure before infectiousness

Susceptible\(S\)

→

Exposed\(E\)

→

Infectious\(I\)

→

Removed\(R\)

\[ S\xrightarrow{\text{infection}}E \xrightarrow{\text{progression}}I \xrightarrow{\text{recovery}}R. \]

An exposed person is infected but not yet infectious in this basic model. Exposure and progression are therefore different biological events.

Why not move directly from \(S\) to \(I\)?

In SIR, a new infection immediately increases \(I\). In SEIR, the same infection event changes:

\[ (S,E,I,R)\longrightarrow(S-1,E+1,I,R). \]

The infectious count changes later, when a progression event occurs:

\[ (S,E,I,R)\longrightarrow(S,E-1,I+1,R). \]

This separation introduces a delay between acquiring infection and becoming able to transmit.

The SEIR state space

For a closed population of size \(N\):

\[ S_n+E_n+I_n+R_n=N. \]

The full integer-valued state is \(X_n=(S_n,E_n,I_n,R_n)\). Its state space is:

\[ \mathcal S= \left\{(s,e,i,r)\in\mathbb Z_{\ge0}^{4}: s+e+i+r=N\right\}. \]

Only three compartment counts are mathematically independent because the fourth can be obtained from conservation. We nevertheless store all four because this makes biological interpretation and validation clearer.

Three events and four outcomes

Under the one-event-per-step assumption, three biological events are possible, plus no change:

Infection

\[(s,e,i,r)\to(s-1,e+1,i,r).\]

Progression

\[(s,e,i,r)\to(s,e-1,i+1,r).\]

Recovery

\[(s,e,i,r)\to(s,e,i-1,r+1).\]

No change leaves \((s,e,i,r)\) unchanged. These four outcomes must be mutually exclusive and their probabilities must sum to 1.

Transition probabilities

Let \(\beta\) be the transmission parameter, \(\sigma\) the exposed-to-infectious progression rate and \(\gamma\) the recovery rate. For a sufficiently short \(\Delta t\):

\[ \begin{aligned} p_{\mathrm{inf}} &=\beta\frac{SI}{N}\Delta t,\\ p_{\mathrm{prog}} &=\sigma E\Delta t,\\ p_{\mathrm{rec}} &=\gamma I\Delta t,\\ p_{\mathrm{stay}} &=1-p_{\mathrm{inf}}-p_{\mathrm{prog}}-p_{\mathrm{rec}}. \end{aligned} \]

Infection still depends on susceptible and infectious people—not exposed people—because this basic SEIR model assumes exposed individuals do not transmit.

Interpret \(\sigma\) correctly

Mean latent period

\[1/\sigma.\]

With \(\sigma=0.2\) per day, the mean latent period is 5 days.

Mean infectious period

\[1/\gamma.\]

With \(\gamma=0.1\) per day, the mean infectious period is 10 days.

These are mean waiting periods implied by constant progression and recovery rates. They do not mean that every individual progresses after exactly 5 days or recovers after exactly 10 days.

The crucial extinction condition

\(I=0\) is not sufficient

If \(I=0\) but \(E>0\), exposed people can progress into the infectious compartment. The epidemic process can continue.

\(E=0\) and \(I=0\)

When both infected compartments are empty, infection, progression and recovery probabilities are all zero. The state is absorbing.

\[ \text{Extinction condition:}\qquad E=0\quad\text{and}\quad I=0. \]

This is why the simulation must use if E == 0 and I == 0, not merely if I == 0.

Probability validity now includes progression

The time-step condition becomes:

\[ p_{\mathrm{inf}}+ p_{\mathrm{prog}}+ p_{\mathrm{rec}} \le1. \]

Adding a compartment adds an event and therefore another contribution to the total event probability. A time step valid for an SIR model is not automatically valid after SEIR progression is added.

A simple sufficient bound over the whole state space is:

\[ \left( \frac{\beta N}{4} +\sigma N +\gamma N \right)\Delta t\le1. \]

The bound is conservative but easy to calculate and guarantees non-negative no-change probability everywhere.

Interactive Python laboratory

Click Run code to simulate one SEIR trajectory. The output reports both exposed and infectious peaks, completed outbreak size, event counts and two graphs.

Interactive PythonDTMC SEIR epidemic

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def seir_probabilities(S, E, I, N, beta, sigma, gamma, dt):
    p_infection = beta * S * I / N * dt
    p_progression = sigma * E * dt
    p_recovery = gamma * I * dt
    p_no_change = (
        1 - p_infection - p_progression - p_recovery
    )

return (
        p_infection,
        p_progression,
        p_recovery,
        p_no_change
    )

def one_seir_transition(
    S, E, I, R, beta, sigma, gamma, dt, rng
):
    N = S + E + I + R

probabilities = seir_probabilities(
        S, E, I, N, beta, sigma, gamma, dt
    )
    p_infection, p_progression, p_recovery, p_no_change = (
        probabilities
    )

if p_no_change < -1e-12:
        raise ValueError("dt gives invalid transition probabilities.")

u = rng.random()

if u < p_infection:
        return S - 1, E + 1, I, R, "infection"
    elif u < p_infection + p_progression:
        return S, E - 1, I + 1, R, "progression"
    elif u < p_infection + p_progression + p_recovery:
        return S, E, I - 1, R + 1, "recovery"
    else:
        return S, E, I, R, "no change"

# Parameters and initial state
beta = 0.30
sigma = 0.20
gamma = 0.10
dt = 0.01
max_days = 250
seed = 11

S, E, I, R = 195, 3, 2, 0
N = S + E + I + R

initial_S = S
initial_E = E
initial_I = I
initial_R = R

if min(S, E, I, R, beta, sigma, gamma) < 0:
    raise ValueError("Counts and rates cannot be negative.")

if N <= 0 or dt <= 0:
    raise ValueError("N and dt must be positive.")

# Conservative probability check over the whole state space
upper_probability_bound = (
    beta * N / 4
    + sigma * N
    + gamma * N
) * dt

if upper_probability_bound > 1:
    raise ValueError("dt is not guaranteed valid over the state space.")

# Store the initial state
rng = np.random.default_rng(seed)
times = [0.0]
S_values = [S]
E_values = [E]
I_values = [I]
R_values = [R]
events = ["initial state"]

max_steps = int(max_days / dt)

for step in range(max_steps):
    S, E, I, R, event = one_seir_transition(
        S, E, I, R, beta, sigma, gamma, dt, rng
    )

times.append((step + 1) * dt)
    S_values.append(S)
    E_values.append(E)
    I_values.append(I)
    R_values.append(R)
    events.append(event)

# I = 0 alone is not an extinction condition in SEIR
    if E == 0 and I == 0:
        break

trajectory = pd.DataFrame({
    "time": times,
    "S": S_values,
    "E": E_values,
    "I": I_values,
    "R": R_values,
    "event": events
})

# Validate every stored state
totals = (
    trajectory["S"] + trajectory["E"]
    + trajectory["I"] + trajectory["R"]
)

if not (totals == N).all():
    raise ValueError("Population conservation failed.")

if (trajectory[["S", "E", "I", "R"]] < 0).any().any():
    raise ValueError("A compartment became negative.")

if not (trajectory["S"].diff().dropna() <= 0).all():
    raise ValueError("Susceptible count increased.")

if not (trajectory["R"].diff().dropna() >= 0).all():
    raise ValueError("Removed count decreased.")

# Calculate realised outcomes
final_time = trajectory["time"].iloc[-1]
final_S = trajectory["S"].iloc[-1]
final_E = trajectory["E"].iloc[-1]
final_I = trajectory["I"].iloc[-1]
final_R = trajectory["R"].iloc[-1]

peak_E = trajectory["E"].max()
peak_E_time = trajectory.loc[
    trajectory["E"].idxmax(), "time"
]
peak_I = trajectory["I"].max()
peak_I_time = trajectory.loc[
    trajectory["I"].idxmax(), "time"
]

new_transmission_infections = initial_S - final_S
people_infected_so_far = (
    initial_E + initial_I + new_transmission_infections
)

if final_E == 0 and final_I == 0:
    stopping_reason = "infection became extinct"
    final_epidemic_size = people_infected_so_far
    attack_proportion = final_epidemic_size / N
else:
    stopping_reason = "the time horizon was reached"
    final_epidemic_size = "not complete"
    attack_proportion = "not complete"

event_counts = (
    trajectory["event"]
    .value_counts()
    .reindex(
        ["infection", "progression", "recovery", "no change"],
        fill_value=0
    )
)

summary = pd.DataFrame({
    "outcome": [
        "Peak exposed count",
        "First exposed peak time",
        "Peak infectious count",
        "First infectious peak time",
        "Final susceptible",
        "Final exposed",
        "Final infectious",
        "Final removed",
        "New transmission infections",
        "Completed final epidemic size",
        "Attack proportion"
    ],
    "value": [
        peak_E,
        round(peak_E_time, 2),
        peak_I,
        round(peak_I_time, 2),
        final_S,
        final_E,
        final_I,
        final_R,
        new_transmission_infections,
        final_epidemic_size,
        (
            round(attack_proportion, 4)
            if final_E == 0 and final_I == 0
            else "not complete"
        )
    ]
})

print("Whole-state-space sufficient bound")
print(f"Upper bound on total event probability: "
      f"{upper_probability_bound:.4f}")
print()
print(f"Seed: {seed}")
print(f"Final time: {final_time:.2f} days")
print(f"Stopping reason: {stopping_reason}")
print()
print("Realised-outcome table")
print(summary.to_string(index=False))
print()
print("Event counts")
print(event_counts.to_string())
print()
print("Last 8 stored states")
print(trajectory.tail(8).to_string(index=False))

# Figure 1: all compartment trajectories
plt.figure(figsize=(9, 5.2))
plt.step(trajectory["time"], trajectory["S"],
         where="post", label="Susceptible", color="#2563eb")
plt.step(trajectory["time"], trajectory["E"],
         where="post", label="Exposed", color="#f59e0b")
plt.step(trajectory["time"], trajectory["I"],
         where="post", label="Infectious", color="#dc2626")
plt.step(trajectory["time"], trajectory["R"],
         where="post", label="Removed", color="#16a34a")
plt.xlabel("Time (days)")
plt.ylabel("Number of people")
plt.title("One stochastic SEIR trajectory")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

# Figure 2: exposed and infectious delay
plt.figure(figsize=(9, 4.8))
plt.step(trajectory["time"], trajectory["E"],
         where="post", label="Exposed", color="#f59e0b")
plt.step(trajectory["time"], trajectory["I"],
         where="post", label="Infectious", color="#dc2626")
plt.axvline(peak_E_time, color="#f59e0b", linestyle="--",
            alpha=0.7, label="First exposed peak")
plt.axvline(peak_I_time, color="#dc2626", linestyle=":",
            alpha=0.8, label="First infectious peak")
plt.xlabel("Time (days)")
plt.ylabel("Number of people")
plt.title("Exposure precedes infectiousness")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

How the random intervals change

With three events, the cumulative selection rule has three boundaries:

if u < p_infection:
    event = "infection"
elif u < p_infection + p_progression:
    event = "progression"
elif u < p_infection + p_progression + p_recovery:
    event = "recovery"
else:
    event = "no change"

The order could be changed if the cumulative boundaries were changed consistently. What matters is that interval lengths equal the corresponding probabilities and that the intervals do not overlap.

Interpret the two peak times

The exposed and infectious peaks are different random outcomes. In a typical growing outbreak, the exposed peak may occur before the infectious peak because individuals must progress through \(E\) before joining \(I\).

The code must calculate both times from the trajectory. We should not assume that their difference equals exactly \(1/\sigma\). The mean latent period describes individual waiting, whereas population peaks depend on the entire changing epidemic.

Final epidemic size in SEIR

At a completed extinction state, \(E=I=0\). If \(R_0=0\), the final epidemic size is:

\[ Z=R_\infty=N-S_\infty. \]

With initially exposed or infectious people, the size of the modelled outbreak includes them:

\[ Z=E_0+I_0+(S_0-S_\infty). \]

The expression \(S_0-S_\infty\) counts new transmission infections. Adding \(E_0+I_0\) includes the infections already present at time zero.

Understand the model-specific code

Code	Meaning
`p_progression = sigma * E * dt`	Calculates the probability of one exposed person becoming infectious.
`return S, E - 1, I + 1, R`	Applies progression without changing the total population.
`if E == 0 and I == 0`	Stops only when neither exposed nor infectious people remain.
`initial_E + initial_I`	Counts infections already present at the beginning.
`trajectory["E"].idxmax()`	Locates the first stored row attaining the exposed peak.
`.reindex([...], fill_value=0)`	Displays every event category even if one event never occurred in this realisation.

Try these experiments

Change the seed and compare the exposed peak, infectious peak and extinction time.
Set S, E, I, R = 197, 3, 0, 0. Confirm that progression can create infectious people even though \(I_0=0\).
Increase sigma to 0.5. Interpret the shorter mean latent period and changed peak timing.
Decrease sigma to 0.1. Compare the separation between exposed and infectious curves.
Increase dt until the conservative whole-state-space check rejects it.
Set max_days = 20. Explain why final epidemic size is reported as incomplete when \(E>0\) or \(I>0\).

Biological interpretation

The exposed compartment delays the appearance of infectious individuals after transmission. This delay can shift the infectious peak and allows hidden infected individuals to remain even when no infectious person is currently observed.

The SEIR model is appropriate only when a distinct latent stage is biologically meaningful. If exposed individuals can also transmit, the infection probability must be modified rather than silently using the basic SEIR assumption.

What this lesson has added

Infection moves people from \(S\) to \(E\), not directly to \(I\).
Progression is a third event with probability \(\sigma E\Delta t\).
The latent period separates infection from infectiousness.
\(I=0\) is not absorbing when exposed people remain.
Extinction requires both \(E=0\) and \(I=0\).
Time-step validity must include infection, progression and recovery probabilities.
Exposed and infectious peaks are distinct stochastic outcomes.

Lesson 10 changes viewpoint completely: instead of selecting one random path, it represents and updates the probability distribution over all states using a transition matrix.