Loops for a complete epidemic trajectory

One transition gives only the next state. We now repeatedly use that one-step rule, store every state and construct one realised epidemic path through time.

The unique purpose of this lesson

A trajectory, also called a sample path or realisation, is one sequence of states produced by the stochastic model:

\[ X_0,X_1,X_2,\ldots,X_K. \]

Lessons 3–5 already explained the probabilities, random intervals and one-state update. The genuinely new programming ideas here are repetition, storage and stopping.

One trajectory is not a probability distribution and is not the average of all possible epidemics. Another random seed can produce a different valid path.

Scenario: follow one outbreak until it ends

Consider a closed SIR population of 200 people:

\[ (S_0,I_0,R_0)=(197,3,0),\qquad N=200. \]

Use \(\beta=0.3\) per day, \(\gamma=0.1\) per day and fixed step \(\Delta t=0.01\) day. We will follow the population until either:

\(I=0\), meaning that infection has become extinct; or
160 days have been simulated.

From one transition to a trajectory

1. InitialiseSet parameters and \(X_0\).

2. StoreRecord time zero and the initial state.

3. RepeatCalculate one next state during each loop.

4. StopEnd at extinction or the time horizon.

The essential relationship is:

\[ X_{n+1}=\operatorname{transition}(X_n,U_n). \]

The next state from one iteration becomes the current state for the following iteration.

Why store the initial state before the loop?

times = [0.0]
S_values = [S]
I_values = [I]
R_values = [R]

The lists begin with \(X_0\), the state before any transition occurs. The loop then appends \(X_1,X_2,\ldots\). If the initial state were omitted, the plotted path would incorrectly begin after the first transition.

The loop structure

for step in range(max_steps):
    result = one_sir_transition(...)
    S, I, R = result["next_state"]

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

    if I == 0:
        break

`for` and `range`

range(max_steps) supplies the integers from 0 to max_steps - 1. The loop body runs once for each supplied integer unless break stops it early.

State replacement

S, I, R = result["next_state"] makes the newly produced state the current state for the next iteration.

`append`

append adds one new value to the end of a list. All four lists grow together and remain aligned by time step.

`break`

When \(I=0\), break immediately leaves the loop. Further steps would only repeat the absorbing disease-free state.

Why use `step + 1` for time?

Python begins range(max_steps) at step = 0. However, the first state produced inside the loop is \(X_1\), occurring at time \(\Delta t\), not at time zero.

\[ t_{n}=n\Delta t. \]

Therefore the time stored after the first transition is:

(step + 1) * dt

This avoids an off-by-one error, where states are attached to the wrong time values.

Stopping rules and their meanings

Biological stopping rule

Stop when I == 0. In this closed SIR model, infection cannot restart because there is no importation and recovered people do not return to susceptibility.

Computational stopping rule

Stop after max_steps even if infection remains. This guarantees that the program finishes and defines the observation window.

If the time horizon is reached while \(I>0\), we cannot say that extinction occurred. We can say only that infection was still present at the end of the chosen observation period.

Choosing a valid number of steps

If the final time is \(T\) and the fixed step is \(\Delta t\), then:

\[ K=\frac{T}{\Delta t}. \]

max_steps = int(max_days / dt)

int converts the result into a whole number because range requires an integer. Here \(160/0.01=16000\) possible steps.

Interactive Python laboratory

Click Run code to generate one complete stochastic trajectory. The result includes a summary, selected rows from the stored table, event counts and a step plot.

Interactive PythonOne complete SIR DTMC trajectory

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

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

p_infection = beta * S * I / N * dt
    p_recovery = gamma * I * dt
    p_no_change = 1 - p_infection - p_recovery

if p_no_change < -1e-12:
        raise ValueError("dt is too large for valid probabilities.")

p_no_change = max(0.0, p_no_change)
    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
beta = 0.3
gamma = 0.1
dt = 0.01
max_days = 160
max_steps = int(max_days / dt)

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

# Reproducible random-number generator
seed = 7
rng = np.random.default_rng(seed)

# Store the initial time and state before the loop
times = [0.0]
S_values = [S]
I_values = [I]
R_values = [R]
events = ["initial state"]

# Repeatedly apply one transition
for step in range(max_steps):
    S, I, R, event = one_sir_transition(
        S, I, R, beta, gamma, dt, rng
    )

time = (step + 1) * dt

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

# Infection extinction is absorbing in this model
    if I == 0:
        break

# Convert aligned lists into one table
trajectory = pd.DataFrame({
    "time": times,
    "S": S_values,
    "I": I_values,
    "R": R_values,
    "event": events
})

# Validate every stored state
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.")

# Summarise the realised trajectory
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_row = trajectory["I"].idxmax()
peak_time = trajectory.loc[peak_row, "time"]

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

print("Trajectory summary")
print(f"Seed: {seed}")
print(f"Stored states: {len(trajectory)}")
print(f"Final time: {final_time:.2f} days")
print(f"Stopping reason: {stopping_reason}")
print(f"Peak infectious count: {peak_I}")
print(f"First time at that peak: {peak_time:.2f} days")
print(f"Final state: ({final_S}, {final_I}, {final_R})")
print()

print("First 10 stored states")
print(trajectory.head(10).to_string(index=False))
print()
print("Last 10 stored states")
print(trajectory.tail(10).to_string(index=False))
print()

event_counts = (
    trajectory["event"]
    .value_counts()
    .reindex(["infection", "recovery", "no change"], fill_value=0)
)
print("Event counts")
print(event_counts.to_string())

plt.figure(figsize=(9, 5.2))
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.xlabel("Time (days)")
plt.ylabel("Number of people")
plt.title("One stochastic SIR trajectory")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

How to read the output table

Column	Meaning
`time`	The model time attached to the stored state.
`S`, `I`, `R`	The integer compartment counts at that time.
`event`	The event that produced this row from the preceding row.

The first row is labelled initial state because no transition produced it. Every later row records infection, recovery or no change.

Why use a step plot?

A DTMC state remains at its recorded value until the next fixed observation time, then may jump to another integer state. Therefore:

plt.step(..., where="post")

is more faithful than joining the values with sloping lines. A sloping line could visually suggest that fractional people change continuously between time points.

Understand the new code

Code	Meaning
`for step in range(max_steps)`	Allows at most `max_steps` repetitions of the one-transition rule.
`times.append(time)`	Adds the new time to the end of the time list.
`break`	Stops the loop early when the biological stopping condition is met.
`pd.DataFrame({...})`	Combines the aligned lists into a labelled rectangular table.
`.iloc[-1]`	Selects the last stored value by its position.
`.idxmax()`	Returns the row label at which the infectious count first reaches its maximum.
`.value_counts()`	Counts how many times each event label occurs.
`.all()`	Checks whether a condition is true for every value being tested.

Three validations after the loop

Aligned storageEvery list receives exactly one value per stored state.

Population conservationEvery row must satisfy \(S+I+R=N\).

Non-negativityNo stored compartment count may be below zero.

The one-step function protects each transition. The checks after the loop protect the complete stored dataset.

Try these experiments

Change seed = 7. Compare the peak, final state and stopping time with the original trajectory.
Use seed = 2. Does early randomness produce a noticeably different path?
Change max_days = 20. If \(I>0\) at the end, explain why the result is not an extinction time.
Start with S, I, R = 199, 1, 0. The outbreak may disappear very early or grow.
Comment out the break statement. Observe that the disease-free state is repeatedly stored until the time horizon.

Biological interpretation

The plotted line is one possible history of the epidemic. Early infection and recovery events can change whether the outbreak grows, how high it peaks and when it becomes extinct.

The final removed count is the number of people infected during this particular realisation because everyone begins susceptible or infectious and removed people do not return to susceptibility.

We should not generalise from one path. Estimating extinction probability, expected peak size or the distribution of final epidemic size requires many independent trajectories. That Monte Carlo task is reserved for Lesson 12.

What this lesson has—and has not—done

You can now:

repeat a one-step DTMC transition with a loop;
store the initial state and every later state in aligned lists;
attach states to the correct discrete times;
stop at extinction or a fixed time horizon;
validate and tabulate a complete trajectory; and
display an integer-valued stochastic path with a step plot.

Next: Lessons 7–9 apply this general trajectory machinery to the distinctive biology and state structures of SIS, SIR and SEIR models. Lesson 7 begins with reinfection in an SIS epidemic.

Loops for a complete epidemic trajectory

The unique purpose of this lesson

Scenario: follow one outbreak until it ends

From one transition to a trajectory

Why store the initial state before the loop?

The loop structure

for and range

State replacement

append

break

Why use step + 1 for time?

Stopping rules and their meanings

Biological stopping rule

Computational stopping rule

Choosing a valid number of steps

Interactive Python laboratory

How to read the output table

Why use a step plot?

Understand the new code

Three validations after the loop

Try these experiments

Biological interpretation

What this lesson has—and has not—done

`for` and `range`

`append`

`break`

Why use `step + 1` for time?