Comparing deterministic, CTMC and SDE results

A fair comparison holds the biological rates, initial state and observation times fixed while recognising that the three approaches represent randomness differently.

Scenario and comparison question

Use an SIS infection in a population of \(N=200\), initially with \(I(0)=3\). All three models use \(\beta=0.25\) and \(\gamma=0.15\) per day and are compared over the same 60-day period.

The comparison asks three separate questions:

How does the deterministic prevalence path compare with stochastic ensemble means?
How wide are the CTMC and SDE outcome distributions?
How similarly do the stochastic models represent the disease-free boundary?

One biology, three mathematical descriptions

Approach	State and path	Main information
Deterministic ODE	Continuous, one smooth path	Average-direction reference
CTMC	Integer, piecewise-constant jumps at random times	Discrete events and extinction
SDE	Continuous, irregular diffusion path	Approximate continuous fluctuations

The parameter values and biological rates must match before outcomes are compared.

Shared SIS drift

\[f(I)=\beta\frac{(N-I)I}{N}-\gamma I.\]

The deterministic model uses \(dI=f(I)dt\). The CTMC uses infection and recovery events separately. The SDE uses the same drift and diffusion \(g(I)=\sqrt{a_1(I)+a_2(I)}\).

Requirements for a fair comparison

Must match	Why
Population and initial state	Different starting risks produce different outcomes.
Infection and recovery rate formulas	The models must represent the same biological mechanisms.
Parameter values and units	A per-day rate cannot be compared with a per-week rate.
Observation times	Outcomes must refer to the same dates.
Boundary and stopping definitions	Extinction must mean the same biological event.

Numerical choices such as Euler–Maruyama \(\Delta t\) must also be reported, even though the CTMC does not use that fixed step.

Interactive Python laboratory

Interactive PythonDeterministic, CTMC and SDE comparison

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

N = 200
I0 = 3
beta = 0.25
gamma = 0.15
T = 60
dt = 0.02
M = 300
seed = 107
seed_sequence = np.random.SeedSequence(seed)
ctmc_seed, sde_seed = seed_sequence.spawn(2)
ctmc_rng = np.random.default_rng(ctmc_seed)
sde_rng = np.random.default_rng(sde_seed)

observation_times = np.arange(0, T + 1)

def drift(I):
    return beta * (N - I) * I / N - gamma * I

# Deterministic Euler solution
fine_times = np.linspace(0, T, int(T / dt) + 1)
I_ode = np.zeros(len(fine_times))
I_ode[0] = I0

for n in range(len(fine_times) - 1):
    I_ode[n + 1] = I_ode[n] + drift(I_ode[n]) * dt

ode_daily = I_ode[(observation_times / dt).astype(int)]

# CTMC paths observed on the common daily grid
I_ctmc = np.zeros((M, len(observation_times)))

for run in range(M):
    time = 0.0
    I = I0
    event_times = [0.0]
    values = [I]

while I > 0 and time < T:
        infection_rate = beta * (N - I) * I / N
        recovery_rate = gamma * I
        total_rate = infection_rate + recovery_rate
        wait = ctmc_rng.exponential(1 / total_rate)

if time + wait > T:
            break

time += wait
        if ctmc_rng.random() < infection_rate / total_rate:
            I += 1
        else:
            I -= 1

event_times.append(time)
        values.append(I)

positions = np.searchsorted(
        event_times, observation_times, side="right"
    ) - 1
    I_ctmc[run] = np.array(values)[positions]

# SDE Euler–Maruyama paths
steps = int(T / dt)
I_sde_fine = np.zeros((M, steps + 1))
I_sde_fine[:, 0] = I0

for n in range(steps):
    current = I_sde_fine[:, n]
    infection_rate = beta * (N - current) * current / N
    recovery_rate = gamma * current
    diffusion = np.sqrt(np.maximum(
        infection_rate + recovery_rate, 0
    ))
    dW = np.sqrt(dt) * sde_rng.normal(0, 1, size=M)
    I_sde_fine[:, n + 1] = np.clip(
        current + drift(current) * dt + diffusion * dW,
        0, N
    )

daily_indices = (observation_times / dt).astype(int)
I_sde = I_sde_fine[:, daily_indices]

ctmc_mean = I_ctmc.mean(axis=0)
sde_mean = I_sde.mean(axis=0)
ctmc_low, ctmc_high = np.quantile(
    I_ctmc, [0.05, 0.95], axis=0
)
sde_low, sde_high = np.quantile(
    I_sde, [0.05, 0.95], axis=0
)

selected_days = [0, 10, 20, 40, 60]
selected_indices = [int(day) for day in selected_days]

time_comparison = pd.DataFrame({
    "day": selected_days,
    "deterministic I": ode_daily[selected_indices],
    "CTMC mean I": ctmc_mean[selected_indices],
    "SDE mean I": sde_mean[selected_indices],
    "CTMC fraction zero": (
        I_ctmc[:, selected_indices] == 0
    ).mean(axis=0),
    "SDE fraction zero": (
        I_sde[:, selected_indices] == 0
    ).mean(axis=0)
})

print("Matched-model comparison at selected times")
print(time_comparison.round(4).to_string(index=False))

plt.figure(figsize=(9.5, 5.4))
plt.fill_between(
    observation_times, ctmc_low, ctmc_high,
    color="#bdd7f5", alpha=0.45,
    label="CTMC 5th–95th percentiles"
)
plt.fill_between(
    observation_times, sde_low, sde_high,
    color="#f2c7c2", alpha=0.42,
    label="SDE 5th–95th percentiles"
)
plt.plot(observation_times, ode_daily,
         color="black", linewidth=2,
         label="Deterministic")
plt.plot(observation_times, ctmc_mean,
         color="#2563eb", label="CTMC mean")
plt.plot(observation_times, sde_mean,
         color="#b42318", label="SDE mean")
plt.xlabel("Time (days)")
plt.ylabel("Infectious population")
plt.title("Same SIS biology, different mathematical descriptions")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()

plt.figure(figsize=(9, 4.9))
bins = np.arange(0, N + 5, 5)
plt.hist(I_ctmc[:, -1], bins=bins, alpha=0.55,
         label="CTMC", color="#2563eb")
plt.hist(I_sde[:, -1], bins=bins, alpha=0.50,
         label="SDE", color="#b42318")
plt.axvline(ode_daily[-1], color="black",
            linestyle="--", label="Deterministic value")
plt.xlabel("Infectious population at day 60")
plt.ylabel("Number of paths")
plt.title("Final-time outcomes")
plt.legend()
plt.grid(axis="y", alpha=0.25)
plt.tight_layout()

Output

Run the code to see the result.

Why the stochastic mean need not equal the deterministic path

The SIS drift is nonlinear because it contains \(I^2\). In general,

\[\mathbb E[f(I)]\ne f(\mathbb E[I]).\]

The deterministic equation evaluates the drift at one deterministic state. The stochastic mean depends on the full distribution, including variance and extinction. Therefore, matching biological rates does not require the deterministic path and stochastic ensemble mean to be identical.

How to interpret agreement

The deterministic path is not expected to equal every stochastic path. Agreement is assessed through ensemble summaries. For sufficiently large populations away from boundaries, the SDE distribution may approximate the CTMC distribution. Near extinction and small counts, discrepancies can be important.

When is the SDE approximation useful?

The diffusion approximation is generally most plausible when event rates and population counts are sufficiently large, changes of one individual are small relative to the population scale, and trajectories remain away from absorbing boundaries.

It becomes less reliable when infectious counts are small, extinction probability is central, discrete superspreading events matter, or the numerical method frequently corrects invalid states. In those cases, the CTMC provides a more direct event-level description.

Do not claim that one model is universally better

Use the ODE for a fast average-direction baseline.
Use the CTMC when integer events, exact absorption or small populations matter.
Use the SDE when continuous fluctuation is a useful and efficient approximation.
Model choice depends on the biological question, scale and required outcomes.

What this lesson adds

You can now compare three models under matched assumptions, align their outputs on common times, distinguish path agreement from distributional agreement, and explain why extinction boundaries challenge diffusion approximations.