Constructing a stochastic SEIR model

The stochastic SEIR model adds a latent progression channel and shows how event-vector programming scales cleanly to more compartments.

Three biological events

Event	Change in \((S,E,I,R)\)	Rate
Infection	\((-1,+1,0,0)\)	\(a_1=\beta SI/N\)
Progression	\((0,-1,+1,0)\)	\(a_2=\sigma E\)
Recovery	\((0,0,-1,+1)\)	\(a_3=\gamma I\)

Four coupled SDEs

\[ dS=-a_1dt-\sqrt{a_1}dW_1, \] \[ dE=(a_1-a_2)dt+\sqrt{a_1}dW_1-\sqrt{a_2}dW_2, \] \[ dI=(a_2-a_3)dt+\sqrt{a_2}dW_2-\sqrt{a_3}dW_3, \] \[ dR=a_3dt+\sqrt{a_3}dW_3. \]

The three Brownian motions are independent event-channel noises. Within a channel, the same increment is used with opposite signs in the source and destination compartments.

What the exposed stage changes

Infection adds to \(E\), not directly to \(I\). Progression transfers exposed people to infectiousness after a random latent period. Thus infection noise first affects \(S,E\), while progression noise links \(E,I\).

Interactive Python laboratory

Interactive PythonStochastic SEIR construction

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

N = 1000
initial = np.array([985.0, 10.0, 5.0, 0.0])
beta = 0.30
sigma = 0.20
gamma = 0.10
T = 140
dt = 0.01
seed = 63

steps = int(round(T / dt))
times = np.linspace(0, T, steps + 1)
rng = np.random.default_rng(seed)

path = np.zeros((steps + 1, 4))
path[0] = initial
corrections = 0

# Columns are infection, progression and recovery events
V = np.array([
    [-1,  0,  0],
    [ 1, -1,  0],
    [ 0,  1, -1],
    [ 0,  0,  1]
], dtype=float)

for n in range(steps):
    S, E, I, R = path[n]

rates = np.array([
        beta * S * I / N,
        sigma * E,
        gamma * I
    ])
    rates = np.maximum(rates, 0)

dW = np.sqrt(dt) * rng.normal(0, 1, size=3)

drift_change = V @ rates * dt
    diffusion_change = V @ (np.sqrt(rates) * dW)
    proposed = path[n] + drift_change + diffusion_change

if (proposed < 0).any():
        corrections += 1
        proposed = np.maximum(proposed, 0)
        proposed *= N / proposed.sum()

path[n + 1] = proposed

S, E, I, R = path.T
totals = path.sum(axis=1)

summary = pd.DataFrame({
    "quantity": [
        "Maximum population error",
        "Minimum compartment value",
        "Peak exposed",
        "Time of exposed peak",
        "Peak infectious",
        "Time of infectious peak",
        "Boundary corrections"
    ],
    "value": [
        np.max(np.abs(totals - N)),
        path.min(),
        E.max(),
        times[np.argmax(E)],
        I.max(),
        times[np.argmax(I)],
        corrections
    ]
})

print(summary.round(5).to_string(index=False))
print()
print("State-change matrix V")
print(pd.DataFrame(
    V,
    index=["S", "E", "I", "R"],
    columns=["infection", "progression", "recovery"]
).to_string())

plt.figure(figsize=(9.5, 5.2))
plt.plot(times, S, label="Susceptible", color="#2563eb")
plt.plot(times, E, label="Exposed", color="#d97706")
plt.plot(times, I, label="Infectious", color="#b42318")
plt.plot(times, R, label="Removed", color="#39815e")
plt.xlabel("Time (days)")
plt.ylabel("Population")
plt.title("One stochastic SEIR diffusion path")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Understand the matrix code

Expression	Meaning
`V @ rates`	Adds the mean effect of all event channels to each compartment.
`np.sqrt(rates) * dW`	Creates one diffusion contribution per event.
`V @ (...)`	Distributes each event’s random contribution with the correct signs.

Population conservation

Every column of \(V\) sums to zero. Therefore both \(Va(X)dt\) and \(V\sqrt{\operatorname{diag}(a(X))}dW\) sum to zero. This is the structural reason the uncorrected SDE conserves \(S+E+I+R=N\).

Interpretation and conditions

The mean latent period is \(1/\sigma\), but the diffusion does not impose a fixed delay.
Extinction requires both \(E=0\) and \(I=0\).
Continuous states approximate counts and are most credible away from small-count boundaries.
The declared projection is a numerical safeguard, not part of the derived SDE.

What this lesson adds

You can now construct an SEIR diffusion with three event channels, use a state-change matrix, preserve population structurally, and explain how latent progression separates infection from infectiousness.