Skip to main content
🧠0% 🤖100%

Stochastic compartmental

Gillespie (exact, in continuous time), tau-leaping (close to exact, for a run of many events), and chain-binomial (in discrete time). One model, several ways of being random.

Recommendation

  1. Stochasticity is a run-time choice with several named solvers, not a property of the model.
  2. Prefer continuous time where the model allows it, because it eliminates rate-to-probability conversion entirely.
  3. Name the approximation. Gillespie, tau-leaping, and chain-binomial are different models with different validity conditions, and the result should say which one produced it.
  4. One sim is one trajectory; the distribution comes from many runs, and running many is one call.
  5. Quantities that are not people are states too. An environmental reservoir or a vector population is a state holding a count, with its own transitions, and disease rates may use that count directly.
  6. Extinction is the main reason to use this paradigm; make it visible (msim.plot(lines=True) shows every run, and a run with branching=True gives the share that died out, msim.branching.p_extinct).

Why

One model, several stochastic methods, already demonstrated

camdl runs the same model text on Gillespie SSA (exact continuous-time), chain-binomial (Euler-multinomial discrete-time), and ODE. epipack runs the same processes deterministically and stochastically, including “a correct inhomogeneous-Poisson Gillespie implementation” for time-varying rates, “a real capability that most frameworks skip or approximate”. odin shares initial() between deriv() and update().

PyRoss is the sharpest evidence for the ladder rather than the binary: it implements Gillespie, tau-leaping, and Gaussian/van Kampen approximations of the same model side by side. MEmilio has a whole sde_* family alongside its ode_* family, separately implemented.

Continuous time removes an error class

SimInf is continuous-time Gillespie throughout, and its review makes the point directly: “Being continuous-time is a real advantage for correctness — the exponential-conversion errors that afflict discrete-time frameworks simply do not arise — and a real limitation for scale, since Gillespie cost grows with event count.” EoN is continuous-time for the same reason.

The implication drawn in the SimInf review: “The lingua franca should treat discrete-time versus continuous-time as a backend choice with a declared cost, not as a property of the model.”

Disease state coupled to an environment

SimInf’s u/v split (integer compartment counts governed by the Markov chain, alongside continuous per-node variables updated by a post-timestep function) is used by its SISe family for a SIS disease coupled to an environmental compartment. Its review: “Environmental reservoirs, accumulated exposure, and vector abundance are all naturally continuous and coupled to a discrete disease process.” The need is real. The lingua franca meets it without a second kind of state: the reservoir is an ordinary state that holds a count instead of people, and rates are written in terms of that count.

Each stochastic solver has conditions under which it is valid

Gillespie’s method is exact, and its cost grows with the number of events. Tau-leaping is valid when the rates are roughly constant over a leap. A chain-binomial run is valid when rate × dt is small. The diffusion approximation is “valid only at moderate-to-large populations” (MEmilio), so it fails exactly where stochasticity was the reason to use it; it is not offered as a solver. An ODE is valid when counts are large and extinction is not the question. Because the model is data, these conditions can be checked rather than documented, and the run says when one is in doubt. The exact method is called 'gillespie' because that is the name the field uses for it (EoN, epipack, SimInf, camdl, and PyRoss all say Gillespie).

The specification

12-backends.md specifies the solvers ('auto', 'ode', 'binomial', 'gillespie', 'tau', and 'abm', with their options as ss.solvers.ode(dwells='exact') and the like), how a value that varies in time is followed within a step, Erlang stages for dwell times, and the notes on small counts, possible extinction, and large steps. States that hold a count (a reservoir, a vector population) are in 03-states-and-transitions.md, and many runs in 09-results.md. The docs show how to use these in Solvers and backends and Branching and rescaling.

Trade-offs

  • Exact Gillespie simulation does not scale, and the point at which tau-leaping pays off depends on the population and the rates. So the choice is the user’s: solver='gillespie' stays exact however many events a run takes, and solver='tau' takes leaps. An extinction question suits the exact solver, since it arises when counts are small and events are few.
  • A reservoir with very large counts would make exact Gillespie simulation expensive. Bacteria in water are shed and decay in numbers far larger than the number of infections, and simulating each one as an event would be wasteful. A state that is not people holds a plain number and is not stepped event by event.
  • Count-based rates are literal. beta * W.n / (W.n + K) is used exactly as written, so the library cannot check that it is a sensible force of infection or adjust it for mixing. That is the intent, and it means the user owns the formula.
  • Continuous-time methods make “results at time t” a resampling question. MEmilio’s interpolate_simulation_result() and SimInf’s tspan both handle this and most frameworks leave it to the user. Here every solver reports its results at the sim’s time points.

Rejected

  • Stochasticity as a model property (MEmilio’s sde_* classes, EpiModel’s icm versus dcm). It is a choice of solver.
  • A single stochastic solver (SimInf: Gillespie only; Epydemix: chain-binomial only). Each is right for a different population size.
  • Substituting an approximation without saying so. If tau-leaping was used, the result says so.
  • A separate kind of continuous state for environmental quantities (SimInf’s v). Ordinary states holding counts, with count-based rates, cover the same cases with no new concept.
  • A stochastic sim that is many trajectories. A sim is one run; ss.parallel(sim, n_runs=500) is many.
  • The diffusion approximation as the default stochastic solver. It fails at small counts, which is where stochasticity matters. (The 'sde' solver has since been removed altogether.)
  • Requiring the user to choose. solver='auto' chooses and prints the choice, and any solver can be asked for by name, the same policy as calibration. (A stochastic solver is always asked for by name.)

Open questions

The open questions from this document are in 19-open-questions.md.