Programming a CTMC SEIR model

The SEIR model separates becoming infected from becoming infectious. A newly infected person first enters a random latent period in the exposed compartment.

Scenario: infection with a latent period

A closed population contains 300 people. Initially, 294 are susceptible, 4 are exposed and 2 are infectious. Exposed individuals carry the infection but are not yet infectious in this basic model.

Susceptible
\(S\)

infection →

Exposed
\(E\)

progression →

Infectious
\(I\)

recovery →

Removed
\(R\)

Important distinction: an infection event increases \(E\), not \(I\). Infectious prevalence rises later, only when a progression event moves an exposed person into \(I\).

Three events and their rates

1. Infection

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

\[a_1=\beta\frac{SI}{N}.\]

2. Progression

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

\[a_2=\sigma E.\]

3. Recovery

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

\[a_3=\gamma I.\]

The total rate is \(a_0=a_1+a_2+a_3\). At each jump, Gillespie’s algorithm selects one and only one of these events.

Interpret the progression parameter

Each exposed individual progresses at rate \(\sigma\) per day. Under the Markov assumption, that person’s latent time is exponentially distributed:

\[T_E\sim\operatorname{Exp}(\sigma),\qquad \mathbb E[T_E]=\frac{1}{\sigma}.\]

For \(\sigma=0.20\) per day, the mean latent period is 5 days. This is a mean, not a fixed five-day delay: individual progression times vary randomly.

When has an SEIR epidemic ended?

Condition	Can infection continue?	Reason
\(I=0\) but \(E>0\)	Yes	Exposed people can progress and create new infectious people.
\(E=0\) but \(I>0\)	Yes	Infectious people can still infect or recover.
\(E=0\) and \(I=0\)	No	All three event rates are zero.

Therefore, the correct biological stopping condition is E == 0 and I == 0, not merely I == 0.

Interactive Python laboratory

Run a complete SEIR epidemic. The program prints event identities that verify the compartment bookkeeping, plots all four state paths, and compares the cumulative infection and progression processes.

Interactive PythonCTMC SEIR model

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

# Initial state
S, E, I, R = 294, 4, 2, 0
initial_S = S
initial_E = E
initial_I = I
N = S + E + I + R

# Per-day parameters
beta = 0.32
sigma = 0.20
gamma = 0.10

time = 0.0
seed = 21
rng = np.random.default_rng(seed)

times = [time]
S_values = [S]
E_values = [E]
I_values = [I]
R_values = [R]
event_rows = []

while E > 0 or I > 0:
    # Recalculate all three event rates
    infection_rate = beta * S * I / N
    progression_rate = sigma * E
    recovery_rate = gamma * I

rates = np.array([
        infection_rate,
        progression_rate,
        recovery_rate
    ])
    total_rate = rates.sum()

if total_rate <= 0:
        raise ValueError(
            "Active infection exists but total rate is not positive."
        )

# Generate the next jump time
    waiting_time = rng.exponential(1 / total_rate)
    time += waiting_time

# Select an event from cumulative rate proportions
    cumulative = np.cumsum(rates / total_rate)
    u = rng.random()

if u < cumulative[0]:
        S -= 1
        E += 1
        event = "infection"
    elif u < cumulative[1]:
        E -= 1
        I += 1
        event = "progression"
    else:
        I -= 1
        R += 1
        event = "recovery"

# Validate the new state
    if min(S, E, I, R) < 0:
        raise ValueError("A compartment became negative.")

if S + E + I + R != N:
        raise ValueError("Population conservation failed.")

times.append(time)
    S_values.append(S)
    E_values.append(E)
    I_values.append(I)
    R_values.append(R)

event_rows.append({
        "time": time,
        "waiting time": waiting_time,
        "event": event,
        "S": S,
        "E": E,
        "I": I,
        "R": R
    })

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

event_counts = event_log["event"].value_counts()
n_infections = int(event_counts.get("infection", 0))
n_progressions = int(event_counts.get("progression", 0))
n_recoveries = int(event_counts.get("recovery", 0))

peak_E = int(trajectory["E"].max())
peak_E_time = trajectory.loc[
    trajectory["E"].idxmax(), "time"
]
peak_I = int(trajectory["I"].max())
peak_I_time = trajectory.loc[
    trajectory["I"].idxmax(), "time"
]
final_size = N - S

print("CTMC SEIR summary")
print(f"Mean latent period: {1 / sigma:.2f} days")
print(f"Mean infectious period: {1 / gamma:.2f} days")
print(f"Epidemic duration: {time:.3f} days")
print(f"Peak E: {peak_E} at {peak_E_time:.3f} days")
print(f"Peak I: {peak_I} at {peak_I_time:.3f} days")
print(f"Final state [S, E, I, R]: [{S}, {E}, {I}, {R}]")
print(f"Final epidemic size: {final_size}")

print()
print("Event counts")
print(event_counts.to_string())

# Event-count identities provide useful checks
print()
print("Bookkeeping checks")
print(
    "Progressions = initial E + infections:",
    n_progressions == initial_E + n_infections
)
print(
    "Recoveries = initial I + progressions:",
    n_recoveries == initial_I + n_progressions
)
print(
    "Final R = recoveries:",
    R == n_recoveries
)

print()
print("First 10 events")
print(event_log.head(10).round(4).to_string(index=False))

# Figure 1: compartment paths
plt.figure(figsize=(9.5, 5.2))
plt.step(times, S_values, where="post",
         label="Susceptible", color="#2563eb")
plt.step(times, E_values, where="post",
         label="Exposed", color="#f59e0b")
plt.step(times, I_values, where="post",
         label="Infectious", color="#dc2626")
plt.step(times, R_values, where="post",
         label="Removed", color="#39815e")
plt.xlabel("Continuous time (days)")
plt.ylabel("Number of people")
plt.title("One complete CTMC SEIR epidemic")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

# Figure 2: cumulative infection and progression events
infection_indicator = (
    event_log["event"] == "infection"
).astype(int)
progression_indicator = (
    event_log["event"] == "progression"
).astype(int)

cumulative_infections = infection_indicator.cumsum()
cumulative_progressions = progression_indicator.cumsum()

plt.figure(figsize=(9.5, 4.8))
plt.step(event_log["time"], cumulative_infections,
         where="post", label="New infections",
         color="#2563eb")
plt.step(event_log["time"], cumulative_progressions,
         where="post", label="Progressions to infectious",
         color="#b45309")
plt.xlabel("Continuous time (days)")
plt.ylabel("Cumulative event count")
plt.title("Infection occurs before infectious progression")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Understand the three-event selection

Random interval	Selected event
\(0\le U<a_1/a_0\)	Infection
\(a_1/a_0\le U<(a_1+a_2)/a_0\)	Progression
\((a_1+a_2)/a_0\le U<1\)	Recovery

np.cumsum(rates / total_rate) creates the upper boundaries of these intervals. The final boundary is 1 apart from tiny floating-point rounding.

Why the event-count checks work

Every person entering \(E\) must eventually progress before the simulation ends. Therefore:

\[\text{progressions}=E(0)+\text{infection events}.\]

Likewise, every person entering \(I\) must eventually recover:

\[\text{recoveries}=I(0)+\text{progression events}.\]

These are not extra model assumptions. They are accounting identities for a completed closed SEIR epidemic and useful checks on the program.

Interpret the two plots

The compartment plot shows that the exposed peak can occur before the infectious peak. The exact timing and separation vary between trajectories.

In the cumulative-event plot, infection events usually lead progression events because people must first enter \(E\). The gap represents people currently waiting in the exposed compartment; initial exposed people explain the initial offset in the bookkeeping identity.

Try these experiments

Reduce sigma to 0.10. The mean latent period becomes 10 days; compare the peak times.
Increase sigma to 1.0. Observe how the model approaches faster SIR-like progression.
Use I = 0 but keep E = 4. Confirm that the program does not stop prematurely.
Try several seeds and compare whether the exposed or infectious peak changes more.

What is unique to this lesson?

You have separated infection from infectiousness, interpreted the exponential latent period, selected among three competing events, used the correct SEIR extinction condition, and verified the simulation with event-count identities. The next lesson moves from trajectory simulation to the CTMC generator matrix.