The Gillespie Algorithm

The chemistry you learned in school is smooth. Concentrations rise and fall along tidy curves, governed by differential equations that assume something quietly enormous: that there are so many molecules that randomness averages away. A beaker holds about $10^{23}$ of them, so the law of large numbers takes over and the math behaves.

But step inside a single living cell. A gene might be present in one or two copies. A signaling protein might number in the tens. Here the smooth picture collapses: a reaction either fires or it doesn't, and when it fires is a matter of chance. Counts don't glide — they jump by whole molecules, at random times.

In 1977, Daniel Gillespie gave us a way to simulate this honestly. Instead of pretending the noise away, his algorithm embraces it — asking, at every instant, two questions of pure probability: which reaction happens next, and how long until it does?

Here is a tiny chemical system with just two reactions: molecules of X are created at a steady rate, and each existing molecule can decay. Together they pull the count toward an equilibrium — but never smoothly.

<p class="hint">{{hint}}</p>
<div class="readout">
  <span>X = <b id="count">0</b></span>
  <span>t = <b id="time">0.00</b></span>
  <span>{{events_label}}: <b id="events">0</b></span>
</div>
<canvas id="plot" width="520" height="240"></canvas>
<div class="btns">
  <button id="run" type="button">{{btn_run}}</button>
  <button id="reset" type="button" class="ghost">{{btn_clear}}</button>
</div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #222; margin: 0; }
.hint { font-size: .9rem; color: #444; margin: 0 0 .7rem; line-height: 1.45; }
.readout { display: flex; gap: 1.2rem; font: 600 15px ui-monospace, monospace; color: #1d3557; margin: .2rem 0 .5rem; }
.readout b { color: #e63946; }
canvas { width: 100%; max-width: 520px; height: auto; background: #f5f8fb; border: 1px solid #cdd9e3; border-radius: 8px; }
.btns { display: flex; gap: .5rem; flex-wrap: wrap; margin-top: .6rem; }
button { font: 600 14px system-ui, sans-serif; padding: .45rem .9rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .5; cursor: default; }

// Code not found

Press Run and watch the count climb and fall in discrete jumps, each at a random moment. The more X there is, the faster decay competes with creation, so the line wobbles around a balance point instead of resting on it. Press Run again: same rules, a brand-new trajectory. That irreducible scatter is exactly what the smooth differential equation throws away — and exactly what matters inside a cell with only a handful of molecules.

The Gillespie algorithm — formally the Stochastic Simulation Algorithm (SSA), published by Daniel Gillespie in 1976–1977 — is not an approximation. It draws sample trajectories that are statistically exact solutions of the chemical master equation, the full probabilistic law governing the system. Here is the whole engine:

• Propensities. For each reaction, compute a rate $a_j$ — its current probability per unit time of firing, given how many molecules are around. The total rate is $$a_0 = \sum_j a_j$$
• When? The waiting time until something happens is exponentially distributed with rate $a_{0}$. Draw a uniform random number $r_{1}$ and set $$\tau = \frac{1}{a_0} \ln \frac{1}{r_1}$$ Advance the clock by $\tau$.
• Which? Pick the reaction that fires by a second random number $r_{2}$, weighted by each reaction's share $a_j / a_0$.
• Update. Change the molecule counts according to that reaction, recompute the propensities, and repeat.

The cost is not measured in a fixed grid of time steps but in events: the simulation leaps directly from one reaction to the next, spending no effort on the quiet stretches in between. When activity is intense, that becomes expensive — countless tiny jumps — which is why faster variants like tau-leaping and the next-reaction method exist. But for systems where every single molecule counts, the SSA is the gold standard. It is the same exact-sampling spirit behind broader Monte Carlo methods.

"Discrete things happen at random rates" describes an astonishing range of systems, and the Gillespie algorithm is the workhorse for all of them:

• Systems biology. Gene expression is famously noisy: a cell with two identical genes can produce wildly different protein levels purely by chance. The SSA is how researchers reproduce and study that noise.
• Epidemiology. In a small outbreak, whether the infection fizzles or explodes hinges on a few early random events — exactly the regime where averaged equations mislead and Gillespie shines.
• Ecology and population dynamics. Small populations can go extinct from bad luck alone; stochastic simulation captures that risk where deterministic models cannot.
• Physics and chemistry. Surface reactions and other "kinetic Monte Carlo" problems use the very same event-driven recipe.

The unifying idea is the same one behind randomized algorithms: when a process is genuinely random, the honest way to study it is to roll the dice — many times — rather than to smooth the dice away.

The Gillespie algorithm makes a deceptively simple bet: when chance truly drives a system, don't average it away — replay it, one random event at a time. Two coin flips per step, an exponential clock and a weighted choice, and out comes a trajectory that is statistically exact.

That is its quiet power. The same noise that smooth equations discard as an error term is, inside a living cell, the main event — the difference between a gene that switches on and one that stays silent. Gillespie's algorithm doesn't fight that randomness; it speaks its language. Like all good randomized methods, it wins not by being certain, but by being honestly, exactly uncertain.