Programming a CTMC SIS model

This lesson applies Gillespie’s algorithm to an infection in which recovery gives no lasting immunity. Its new biological idea is reinfection: a recovered individual immediately becomes susceptible again.

Scenario: an infection with no lasting immunity

Consider 100 people in a community. Initially, 5 are infectious and 95 are susceptible. After recovery, a person can be infected again. This behaviour can represent some bacterial infections and other infections for which immunity is weak or short-lived.

Susceptible
\(S\)

infection →
← recovery

Infectious
\(I\)

Because the population is closed, \(S+I=N\). Therefore, knowing \(I\) automatically gives \(S=N-I\).

One-dimensional state: the CTMC can be programmed using the single integer \(I\), with state space \(\{0,1,2,\ldots,N\\).

Two events and their rates

Infection

\(I\) increases by 1.

\[I\to I+1,\qquad a_{\mathrm{inf}}(I)=\beta\frac{(N-I)I}{N}.\]

The rate is zero when nobody is infectious or nobody is susceptible.

Recovery

\(I\) decreases by 1.

\[I\to I-1,\qquad a_{\mathrm{rec}}(I)=\gamma I.\]

Recovery returns the individual to the susceptible compartment.

The total rate is \(a_0(I)=a_{\mathrm{inf}}(I)+a_{\mathrm{rec}}(I)\). Gillespie’s algorithm generates the next waiting time from \(\operatorname{Exp}(a_0(I))\) and then chooses infection or recovery in proportion to these rates.

Boundary states and extinction

State	Possible events	Meaning
\(I=0\)	None	Infection is extinct. This is absorbing because the model has no importation.
\(0<I<N\)	Infection and recovery	Both types of event may be possible.
\(I=N\)	Recovery only	There are no susceptible people, so the infection rate is zero.

These conditions are produced naturally by the rate formulas; separate artificial probabilities are not required.

Threshold intuition

Near the beginning, when almost everyone is susceptible, each infectious person produces infections at approximately rate \(\beta\) and recovers at rate \(\gamma\). The basic reproduction number is

\[R_0=\frac{\beta}{\gamma}.\]

If \(R_0>1\), infection initially has an upward tendency. This does not guarantee a major outbreak: the CTMC can still reach \(I=0\) through chance.

Interactive Python laboratory

Click Run code. The program simulates one SIS trajectory, prints its event log and summary, and draws the irregular-time step path. Change the seed to see a different possible epidemic.

Interactive PythonCTMC SIS model

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

# Population and biological parameters
N = 100
I = 5
beta = 0.30
gamma = 0.10

time = 0.0
max_time = 120.0
seed = 7
rng = np.random.default_rng(seed)

# Store the initial state
times = [time]
I_values = [I]
event_rows = []
stopping_reason = "time horizon reached"

while time < max_time:
    S = N - I

# Recalculate the two rates in the current state
    infection_rate = beta * S * I / N
    recovery_rate = gamma * I
    total_rate = infection_rate + recovery_rate

# I = 0 gives total_rate = 0: extinction
    if total_rate == 0:
        stopping_reason = "infection became extinct"
        break

# Generate the next continuous waiting time
    waiting_time = rng.exponential(1 / total_rate)

# Do not perform an event outside the observation window
    if time + waiting_time > max_time:
        break

# Choose infection or recovery by its rate proportion
    u = rng.random()
    if u < infection_rate / total_rate:
        I += 1
        event = "infection"
    else:
        I -= 1
        event = "recovery"

time += waiting_time

# Check that the state remains biologically valid
    if not 0 <= I <= N:
        raise ValueError("I left the valid state space.")

times.append(time)
    I_values.append(I)
    event_rows.append({
        "time": time,
        "waiting time": waiting_time,
        "event": event,
        "S": N - I,
        "I": I
    })

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

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

print("CTMC SIS summary")
print(f"R0 = {beta / gamma:.2f}")
print(f"Stopping reason: {stopping_reason}")
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 recorded I: {I}")
print(f"Final recorded S: {N - I}")

if len(event_log) > 0:
    print()
    print("Event counts")
    print(event_log["event"].value_counts().to_string())
    print()
    print("First 10 events")
    print(event_log.head(10).round(4).to_string(index=False))

# Extend the displayed state to the horizon without adding an event
plot_times = times.copy()
plot_I = I_values.copy()
if stopping_reason == "time horizon reached" and time < max_time:
    plot_times.append(max_time)
    plot_I.append(I)

plt.figure(figsize=(9, 5))
plt.step(plot_times, plot_I, where="post",
         color="#b42318", linewidth=1.7)
plt.axhline(N * (1 - gamma / beta),
            color="#39815e", linestyle="--",
            label="Deterministic endemic level")
plt.xlabel("Continuous time (days)")
plt.ylabel("Infectious people, I")
plt.title("One CTMC SIS epidemic trajectory")
plt.ylim(0, N)
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Understand the important code

Code	Meaning
`S = N - I`	Uses population conservation, so only \(I\) must be stored as the state.
`rng.exponential(1 / total_rate)`	NumPy expects the exponential mean, which is \(1/a_0\), not the rate \(a_0\).
`infection_rate / total_rate`	Converts the infection rate into the conditional probability that the next event is infection.
`0 <= I <= N`	Checks the required SIS state-space condition.
`where="post"`	Keeps the state constant after each jump until the next event.

Interpret the outcome

A rising path means infection events are temporarily occurring more often than recoveries. A falling path means recoveries dominate over that part of this particular random history.

The dashed line is the positive deterministic equilibrium \(I^*=N(1-\gamma/\beta)\) when \(R_0>1\). It is a reference level, not an absorbing CTMC state. A finite stochastic SIS process can fluctuate around it for a long time and still eventually become extinct.

Try these biological experiments

Set I = 1 and try several seeds. Record how often extinction occurs quickly.
Set beta = 0.08, giving \(R_0<1\), and observe the stronger downward tendency.
Increase N while keeping the initial infected proportion similar. Compare the relative size of fluctuations.
Remove the fixed seed only after you understand reproducibility: the result will then vary on every run.

What is unique to this lesson?

You have used one state variable, represented recovery as a return to susceptibility, and distinguished a deterministic endemic level from stochastic extinction. The next lesson changes the biology: recovered people remain removed, so reinfection is no longer possible and final epidemic size becomes meaningful.