Programming a DTMC SIR model

The SIR model describes an infection followed by removal with lasting immunity. Recovery is irreversible, so the final removed population records how many people experienced infection during the modelled outbreak.

The new biological idea: irreversible recovery

Susceptible\(S\)

→

Infectious\(I\)

→

Removed\(R\)

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

A removed person does not return to susceptibility. Depending on the application, “removed” may mean recovered and immune, isolated, or otherwise no longer able to transmit.

How SIR differs from SIS

Feature	SIS	SIR
After recovery	\(I\to S\)	\(I\to R\)
Reinfection	Allowed	Not allowed in the basic model
State reduction	One count: \(I\)	Two counts are needed, such as \((S,I)\)
Disease-free states	Only \(I=0,S=N\)	Many states with \(I=0\)
Final epidemic size	No permanent removed total	Recorded by the final removed count

The SIR state and state space

For a closed population:

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

We may store the full state \((S_n,I_n,R_n)\) because it is biologically easy to read. Mathematically, only two counts are independent because:

\[ R_n=N-S_n-I_n. \]

Using \((S,I)\), the state space is the finite triangular set:

\[ \mathcal S= \{(s,i):s\ge0,\ i\ge0,\ s+i\le N\}. \]

Not every pair from \(\{0,\ldots,N\}^2\) is possible. For example, \(S=150,I=100\) is impossible when \(N=200\) because their sum already exceeds the population.

Possible transitions from ((s,i,r))

Infection

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

One susceptible person becomes infectious.

Recovery

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

One infectious person is removed permanently from transmission.

No change

\[ (s,i,r)\longrightarrow(s,i,r). \]

No event occurs during the fixed step.

Impossible reverse movements

The basic model has no \(R\to S\), \(R\to I\) or \(I\to S\) transition.

Transition probabilities

At current state \((s,i,r)\), the short-step probabilities are:

\[ \begin{aligned} p_{\mathrm{inf}}(s,i) &=\beta\frac{si}{N}\Delta t,\\ p_{\mathrm{rec}}(i) &=\gamma i\Delta t,\\ p_{\mathrm{stay}}(s,i) &=1-p_{\mathrm{inf}}(s,i)-p_{\mathrm{rec}}(i). \end{aligned} \]

The removed count does not appear directly in these probabilities. It still matters because population conservation restricts how large \(s\) and \(i\) can be.

The absorbing boundary is a whole set of states

If \(i=0\), then infection and recovery probabilities are both zero. Therefore every state:

\[ (s,0,N-s),\qquad s=0,1,\ldots,N, \]

is absorbing.

These states represent different completed outbreak histories. A large final susceptible count describes a small outbreak; a small final susceptible count describes a large outbreak.

Quantities obtained from one completed SIR path

Peak size\(I_{\max}\): the largest infectious count in this realisation.

Peak timeThe first time at which this path reaches \(I_{\max}\).

Final epidemic sizeThe number ever infected by the end of this realisation.

If initially \(R_0=0\), everyone who has experienced infection is either still infectious or removed. At extinction \(I_\infty=0\), so:

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

If some people are already removed initially, the number infected during the simulated outbreak is:

\[ Z=R_\infty-R_0=S_0-S_\infty+I_0. \]

The attack proportion is \(Z/N\). It is a realised outcome for one trajectory, not yet an expected value.

Interactive Python laboratory

The code generates one SIR outbreak, checks all stored states, calculates the peak and final epidemic size, and produces both a time trajectory and a state-space path.

Interactive PythonDTMC SIR outbreak

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

def sir_probabilities(S, I, N, beta, gamma, dt):
    p_infection = beta * S * I / N * dt
    p_recovery = gamma * I * dt
    p_no_change = 1 - p_infection - p_recovery

return p_infection, p_recovery, p_no_change

def one_sir_transition(S, I, R, beta, gamma, dt, rng):
    N = S + I + R
    p_infection, p_recovery, p_no_change = sir_probabilities(
        S, I, N, beta, gamma, dt
    )

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

u = rng.random()

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

# Model parameters and initial state
beta = 0.30
gamma = 0.10
dt = 0.01
max_days = 200
seed = 9

S, I, R = 197, 3, 0
N = S + I + R
initial_S = S
initial_I = I
initial_R = R

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

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

# A sufficient whole-state-space check:
# infection rate <= beta*N/4 and recovery rate <= gamma*N
upper_probability_bound = (
    beta * N / 4 + gamma * N
) * dt

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

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

max_steps = int(max_days / dt)

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

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

if I == 0:
        break

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

# Validate the complete stored path
totals = trajectory["S"] + trajectory["I"] + trajectory["R"]

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

if (trajectory[["S", "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.")

# Realised outcomes
final_time = trajectory["time"].iloc[-1]
final_S = trajectory["S"].iloc[-1]
final_I = trajectory["I"].iloc[-1]
final_R = trajectory["R"].iloc[-1]

peak_I = trajectory["I"].max()
peak_time = trajectory.loc[
    trajectory["I"].idxmax(), "time"
]

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

if final_I == 0:
    final_epidemic_size = people_infected_so_far
    attack_proportion = final_epidemic_size / N
else:
    final_epidemic_size = "not complete"
    attack_proportion = "not complete"

if final_I == 0:
    stopping_reason = "infection became extinct"
else:
    stopping_reason = "the time horizon was reached"

summary = pd.DataFrame({
    "outcome": [
        "Peak infectious count",
        "First peak time (days)",
        "Final susceptible",
        "Final infectious",
        "Final removed",
        "New transmission infections",
        "People infected so far",
        "Completed final epidemic size",
        "Attack proportion"
    ],
    "value": [
        peak_I,
        round(peak_time, 2),
        final_S,
        final_I,
        final_R,
        new_transmission_infections,
        people_infected_so_far,
        final_epidemic_size,
        (
            round(attack_proportion, 4)
            if final_I == 0
            else "not complete"
        )
    ]
})

print("Whole-state-space sufficient bound")
print(f"Upper bound on 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("Last 8 stored states")
print(trajectory.tail(8).to_string(index=False))

# Figure 1: compartment trajectories
plt.figure(figsize=(9, 5))
plt.step(trajectory["time"], trajectory["S"],
         where="post", label="Susceptible", color="#2563eb")
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.axvline(peak_time, color="#7c3aed", linestyle="--",
            alpha=0.7, label="First peak time")
plt.xlabel("Time (days)")
plt.ylabel("Number of people")
plt.title("One stochastic SIR trajectory")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

# Figure 2: path through the two-dimensional state space
plt.figure(figsize=(7, 5.5))
plt.step(trajectory["S"], trajectory["I"],
         where="post", color="#7c3aed")
plt.scatter([initial_S], [initial_I],
            color="#16a34a", s=70, label="Initial state", zorder=3)
plt.scatter([final_S], [final_I],
            color="#dc2626", s=70, label="Final state", zorder=3)
plt.xlabel("Susceptible S")
plt.ylabel("Infectious I")
plt.title("The realised path through SIR state space")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Why use a sufficient probability bound?

Across the SIR state space:

\[ SI\le\frac{N^2}{4}, \qquad I\le N. \]

Therefore:

\[ p_{\mathrm{inf}}+p_{\mathrm{rec}} \le \left(\frac{\beta N}{4}+\gamma N\right)\Delta t. \]

This bound is deliberately conservative because both separate maxima need not occur at the same state. If the bound is at most 1, every state is guaranteed to have a non-negative no-change probability.

Why check monotonicity?

In the basic closed SIR model:

\[ S_{n+1}\le S_n, \qquad R_{n+1}\ge R_n. \]

The susceptible count can fall or stay unchanged but cannot increase. The removed count can rise or stay unchanged but cannot decrease. The code uses .diff() to calculate consecutive changes and checks these model-specific directions.

Read the state-space graph

The first graph uses time on its horizontal axis. The second graph removes time from the axes and plots each realised pair \((S_n,I_n)\).

An infection moves the state one unit left and one unit up.
A recovery moves it one unit down while \(S\) stays fixed.
No change leaves the point in the same position.
The path ends on the horizontal absorbing boundary \(I=0\) if extinction occurs.

Final epidemic size when the horizon is too short

If the simulation stops with \(I>0\), final_R is not yet the completed final epidemic size. Some currently infectious people may recover, and additional infections may occur later.

For this reason, the code reports the stopping reason. A completed final size should be claimed only when extinction has occurred or when an explicitly justified approximation is being used.

Understand the model-specific code

Code	Meaning
`initial_R = R`	Retains the starting removed count so infections occurring during the simulation can be separated from earlier removals.
`trajectory["S"].diff()`	Calculates each susceptible count minus the preceding susceptible count.
`.dropna()`	Removes the undefined first difference, because the first row has no preceding row.
`initial_S - final_S`	Counts new transmission infections, because only infection can reduce \(S\).
`initial_I + new_transmission_infections`	Counts everyone infected in this outbreak so far, including those infectious initially.
`final_epidemic_size / N`	Calculates the realised attack proportion only after extinction completes the outbreak.
`plt.scatter(...)`	Marks the initial and final points on the state-space path.

Try these experiments

Change the seed. Compare peak size, peak time and final epidemic size.
Set initial_I = 1 by using S, I, R = 199, 1, 0. Look for early extinction.
Use beta = 0.08 and gamma = 0.10. Compare the outbreak outcome when \(\beta/\gamma<1\).
Set max_days = 20. Explain why the displayed removed count may not be the completed final size.
Use S, I, R = 147, 3, 50. Confirm that new infections are measured relative to the initial removed count.

Biological interpretation

In this realisation, infection consumes susceptible people and recovery accumulates removed people. Early random events influence whether the outbreak disappears quickly or reaches a substantial peak.

Permanent removal gives the SIR process a direction that SIS does not have. Once susceptibility has been lost, it is not restored. This makes the final susceptible and removed counts informative summaries of the completed outbreak.

Peak size, peak time and final epidemic size from one run remain random outcomes. Their probability distributions require many independent trajectories, introduced in Lesson 12.

What this lesson has added

SIR recovery is irreversible and prevents reinfection.
The independent state can be represented by \((S,I)\) in a triangular state space.
Every state with \(I=0\) is absorbing.
\(S\) is non-increasing and \(R\) is non-decreasing.
A completed path gives a realised final epidemic size and attack proportion.
A state-space path displays transitions differently from a time trajectory.

Lesson 9 adds an exposed compartment and a new progression event to construct a DTMC SEIR model.