Stochastic CTMC Models

Epimodels provides stochastic epidemic models based on Continuous-Time Markov Chains (CTMCs), solved using the Gillespie Stochastic Simulation Algorithm (SSA). These models are the stochastic counterparts of the deterministic ODE models and are especially useful when population sizes are small or when extinction risk analysis is required.

Available Models

Four classic models are available in stochastic form:

Model

Compartments

Description

SIR

S, I, R

Classic susceptible-infectious-removed

SIS

S, I

No immunity, reinfection possible

SIRS

S, I, R

Waning immunity

SEIR

S, E, I, R

Latent period (exposed)

All models inherit from CTMCModel, which extends BaseModel with stochastic simulation capabilities.

Quick Start

from epimodels.stochastic.CTMC import SIR

model = SIR()
model(
    inits=[990, 10, 0],
    trange=[0, 100],
    totpop=1000,
    params={'beta': 0.3, 'gamma': 0.1},
    reps=100,
    seed=42,
)
model.plot_traces()

Simulation Parameters

CTMC models accept all standard model parameters plus stochastic-specific ones:

Parameter

Default

Description

reps

1

Number of independent stochastic trajectories

seed

None

Random seed for reproducibility

n_jobs

1

Number of parallel processes for replicate execution

n_points

1000

Number of time points in the output grid

solver

None

Custom solver (defaults to GillespieSolver)

Parallel Execution

For large numbers of replicates, set n_jobs > 1 to use multiprocessing:

model = SIR()
model([990, 10, 0], [0, 100], 1000,
      {'beta': 0.3, 'gamma': 0.1},
      reps=1000, n_jobs=4, seed=42)

Accessing Results

Single Replicate

get_replicate() returns a single replicate’s trajectory as a dictionary:

rep_0 = model.get_replicate(0)
print(rep_0['time'])  # Time points
print(rep_0['I'])     # Infectious trajectory

Summary Statistics

Compute summary statistics across all replicates:

# Mean trajectory across replicates
mean_traces = model.get_mean()
plt.plot(mean_traces['time'], mean_traces['I'], label='Mean I')

# Variance across replicates
var_traces = model.get_variance()
plt.fill_between(
    mean_traces['time'],
    mean_traces['I'] - np.sqrt(var_traces['I']),
    mean_traces['I'] + np.sqrt(var_traces['I']),
    alpha=0.3, label='±1 std'
)

Quantile Bands

get_quantiles() returns quantile envelopes, useful for credible/confidence intervals:

# 95% credible interval
q = model.get_quantiles([0.025, 0.5, 0.975])

plt.plot(q['time'], q['I_0.5'], 'b-', label='Median')
plt.fill_between(q['time'], q['I_0.025'], q['I_0.975'],
                 alpha=0.3, label='95% CI')

Event Times

get_event_times() returns the times at which specific events occurred across replicates:

peak_times = model.get_event_times('peak_detected')
print(f"Peak occurred at t={np.mean(peak_times):.1f} ± {np.std(peak_times):.1f}")

DataFrame Export

Export a specific replicate to a pandas DataFrame:

df = model.to_dataframe(replicate=0)

The Gillespie Solver

GillespieSolver implements the Gillespie Direct Method for exact stochastic simulation. At each step:

  1. Compute propensities for all transition events

  2. Draw time to next event from an exponential distribution

  3. Select which event fires (weighted random selection)

  4. Update the state vector

from epimodels.stochastic.CTMC import SIR, GillespieSolver

solver = GillespieSolver()
model = SIR()
model([990, 10, 0], [0, 100], 1000,
      {'beta': 0.3, 'gamma': 0.1},
      reps=50, solver=solver)

When to Use CTMC Models

Scenario

Notes

Small populations (< 10,000)

Stochastic effects are significant

Extinction risk analysis

Only stochastic models can capture extinction

Confidence intervals

Run multiple replicates for uncertainty

Large populations (> 100,000)

ODE models may suffice; CTMC is slow for many replicates

API Reference

Models

Solvers