Constructing a stochastic SIR model

The stochastic SIR model has two independent event channels, but their noise contributions must be shared between the compartments changed by each event.

Start from two event channels

Event	State-change vector \((S,I,R)\)	Rate
Infection	\((-1,+1,0)\)	\(a_1=\beta SI/N\)
Recovery	\((0,-1,+1)\)	\(a_2=\gamma I\)

Each event produces both a mean contribution and a random contribution in its own state-change direction.

The stochastic SIR equations

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

\(W_1\) and \(W_2\) are independent Brownian motions representing infection and recovery fluctuations.

Why the same increment must be shared

An infection removes exactly what it adds: the term \(-\sqrt{a_1}dW_1\) in \(S\) is paired with \(+\sqrt{a_1}dW_1\) in \(I\). Likewise, recovery noise leaves \(I\) and enters \(R\). Therefore

\[dS+dI+dR=0.\]

Using separate random numbers for the two sides of one event would destroy population conservation.

Interactive Python laboratory

Interactive PythonStochastic SIR construction

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

N = 1000
S0, I0, R0 = 990.0, 10.0, 0.0
beta = 0.28
gamma = 0.10
T = 100
dt = 0.01
seed = 52

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

S = np.zeros(steps + 1)
I = np.zeros(steps + 1)
R = np.zeros(steps + 1)
S[0], I[0], R[0] = S0, I0, R0
corrections = 0

for n in range(steps):
    infection_rate = beta * S[n] * I[n] / N
    recovery_rate = gamma * I[n]

dW_infection = np.sqrt(dt) * rng.normal()
    dW_recovery = np.sqrt(dt) * rng.normal()

infection_change = (
        infection_rate * dt
        + np.sqrt(max(infection_rate, 0)) * dW_infection
    )
    recovery_change = (
        recovery_rate * dt
        + np.sqrt(max(recovery_rate, 0)) * dW_recovery
    )

proposed = np.array([
        S[n] - infection_change,
        I[n] + infection_change - recovery_change,
        R[n] + recovery_change
    ])

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

S[n + 1], I[n + 1], R[n + 1] = proposed

totals = S + I + R

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

print(summary.round(6).to_string(index=False))

plt.figure(figsize=(9.5, 5.2))
plt.plot(times, S, label="Susceptible", color="#2563eb")
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 SIR diffusion path")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Matrix viewpoint

With state-change matrix \(V\) and rate vector \(a(X)\), the event-derived diffusion can be written

\[dX=Va(X)dt+V\sqrt{\operatorname{diag}(a(X))}dW.\]

This compact form becomes especially useful when SEIR adds a third event channel.

Conditions and interpretation

Rates must be non-negative in the biological state space.
The two Brownian increments are independent, but each increment is shared across compartments affected by its event.
Continuous SDE states are approximations to integer counts.
The displayed non-negativity projection is declared and counted; it is not part of the original SDE.

What this lesson adds

You can now construct a multivariable SIR diffusion from event vectors, preserve population through shared noise, and distinguish independent event channels from independent compartment noise.