Programming a CTMC SIR model

The SIR model changes one biological assumption from SIS: recovery gives lasting removal from infection. This creates an accumulating removed compartment and a well-defined final epidemic size.

Scenario: immunity after recovery

A closed community contains 200 people. Three are initially infectious and the other 197 are susceptible. An infectious person either infects a susceptible person or recovers. Recovered people do not return to susceptibility during the modelled period.

Susceptible
\(S\)

infection →

Infectious
\(I\)

recovery →

Removed
\(R\)

Key consequence: every infection permanently moves one person out of \(S\), and every recovery permanently moves one person into \(R\). Hence \(R\) never decreases and the epidemic cannot persist indefinitely through reinfection.

State changes and event rates

Infection

\[(S,I,R)\to(S-1,I+1,R),\qquad a_1=\beta\frac{SI}{N}.\]

One susceptible person becomes infectious. Total population is unchanged.

Recovery

\[(S,I,R)\to(S,I-1,R+1),\qquad a_2=\gamma I.\]

One infectious person enters the removed class permanently.

The total event rate is \(a_0=a_1+a_2\). When \(I=0\), both rates are zero, so the epidemic has ended.

Why SIR needs more state information than SIS

In SIS, \(S=N-I\), so \(I\) alone determines the state. In SIR, the same infectious count can occur with different susceptible counts and therefore different infection rates.

State	\(I\)	\(S\)	Infection rate
\((190,5,5)\)	5	190	\(\beta(190)(5)/N\)
\((50,5,145)\)	5	50	\(\beta(50)(5)/N\)

Thus we store \(S\), \(I\) and \(R\), although conservation \(S+I+R=N\) means only two are mathematically independent.

Two important outcomes

Peak infectious count

The largest value of \(I\) along one random trajectory. Its time is also random.

Final epidemic size

The total number infected during the outbreak:

\[N-S_{\mathrm{final}}=R_{\mathrm{final}}-R_0^{\mathrm{initial}}.\]

Here the initial removed count is zero, so final epidemic size equals final \(R\). Do not confuse the initial compartment value \(R(0)\) with the reproduction number, commonly also written \(\mathcal R_0\).

Interactive Python laboratory

Run one complete SIR epidemic. The first figure shows the compartment path; the second shows how the cumulative number ever infected grows until extinction.

Interactive PythonCTMC SIR model

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

# Initial population
S, I, R = 197, 3, 0
N = S + I + R

# Per-day parameters
beta = 0.25
gamma = 0.10

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

# Store the initial state
times = [time]
S_values = [S]
I_values = [I]
R_values = [R]
event_rows = []

while I > 0:
    # Rates must use the current S and I
    infection_rate = beta * S * I / N
    recovery_rate = gamma * I
    total_rate = infection_rate + recovery_rate

if total_rate <= 0:
        raise ValueError("A positive-I state needs a positive rate.")

# Time until the next event
    waiting_time = rng.exponential(1 / total_rate)
    time += waiting_time

# Select and apply exactly one event
    if rng.random() < infection_rate / total_rate:
        S -= 1
        I += 1
        event = "infection"
    else:
        I -= 1
        R += 1
        event = "recovery"

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

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

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

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

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

peak_I = int(trajectory["I"].max())
peak_time = trajectory.loc[
    trajectory["I"].idxmax(), "time"
]
final_size = N - S
attack_proportion = final_size / N

print("CTMC SIR summary")
print(f"Reproduction-number reference: beta/gamma = {beta/gamma:.2f}")
print(f"Epidemic duration: {time:.3f} days")
print(f"Number of events: {len(event_log)}")
print(f"Peak I: {peak_I}")
print(f"First time at peak: {peak_time:.3f} days")
print(f"Final state [S, I, R]: [{S}, {I}, {R}]")
print(f"Final epidemic size: {final_size}")
print(f"Attack proportion: {attack_proportion:.1%}")

print()
print("Event counts")
print(event_log["event"].value_counts().to_string())
print()
print("Last 8 events")
print(event_log.tail(8).round(4).to_string(index=False))

# Figure 1: all compartments
plt.figure(figsize=(9.5, 5.2))
plt.step(times, S_values, where="post",
         label="Susceptible", color="#2563eb")
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 SIR epidemic")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

# Figure 2: cumulative number ever infected
ever_infected = N - trajectory["S"]

plt.figure(figsize=(9.5, 4.8))
plt.step(trajectory["time"], ever_infected,
         where="post", color="#6d3fc0")
plt.scatter([time], [final_size],
            color="#174c34", zorder=3,
            label=f"Final size = {final_size}")
plt.xlabel("Continuous time (days)")
plt.ylabel("Cumulative number ever infected")
plt.title("Growth of the epidemic's final size")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Understand the stopping rule

while I > 0:

The loop runs while at least one infectious individual can cause infection or recovery. Once \(I=0\), no infectious individual remains, both rates are zero, and the outbreak is absorbed. Unlike the preceding general Gillespie lesson, no arbitrary time horizon is needed to complete this closed SIR outbreak.

Understand the outcome code

Code	Biological meaning
`final_size = N - S`	Everyone who is no longer susceptible has experienced infection.
`attack_proportion = final_size / N`	Proportion of the initial population infected during this outbreak.
`trajectory["I"].idxmax()`	Returns the row at which the largest recorded infectious count first occurs.
`event_log.tail(8)`	Displays the final eight jumps leading to extinction.

Biological interpretation

At first, many susceptible people are available, so infection may outcompete recovery. Later, susceptible depletion reduces \(\beta SI/N\) even if infectious people remain. Eventually recoveries remove the last infectious individuals.

One trajectory does not provide the probability distribution of peak size, duration or final size. It is one possible outbreak. A later lesson will simulate many independent trajectories to estimate those distributions.

Try these experiments

Set I = 1 and S = 199. Try several seeds and look for early extinction.
Reduce beta below gamma. Compare peak and final size.
Start with 50 people already removed while keeping \(N=200\). Examine how fewer susceptible people change the infection rate.
Check the identity: number of infection events plus the initial infectious count equals final \(R\).

What is unique to this lesson?

You have programmed permanent removal, susceptible depletion, certain eventual epidemic termination, and final epidemic size. The next SEIR lesson adds a latent exposed stage: infection and infectiousness no longer begin at the same event.