Introduction

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 102310^{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?

Watch the Molecules Jump

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">A tiny chemical system: molecules of <b>X</b> are <b>created</b> at a steady rate, and each one can <b>decay</b>. The Gillespie algorithm draws two random numbers per step — <i>how long</i> until the next reaction, and <i>which</i> one fires — so the count jumps in discrete, random steps. Press <b>Run</b> for a fresh trajectory.</p>
<div class="readout">
  <span>X = <b id="count">0</b></span>
  <span>t = <b id="time">0.00</b></span>
  <span>events: <b id="events">0</b></span>
</div>
<canvas id="plot" width="520" height="240"></canvas>
<div class="btns">
  <button id="run" type="button">Run</button>
  <button id="reset" type="button" class="ghost">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; }

// Gillespie Stochastic Simulation Algorithm (Gillespie, 1976/1977).
// System: 0 -> X  (creation, rate k1)   and   X -> 0  (decay, rate k2 per molecule).
// Propensities: a1 = k1 (constant),  a2 = k2 * X.   a0 = a1 + a2.
// Per step: tau = (1/a0) * ln(1/r1) is the exponential waiting time;
// r2 picks which reaction fires, weighted by a1/a0 vs a2/a0.
// The deterministic steady state is X* = k1/k2; the stochastic count fluctuates around it.
const k1 = 10;          // creation rate (molecules per unit time)
const k2 = 1;           // decay rate per molecule
const T_MAX = 12;       // total simulated time
const cv = document.getElementById('plot');
const ctx = cv.getContext('2d');
const W = cv.width, H = cv.height, PAD = 28, YMAX = 24;
let timer = null;

function axes(){
  ctx.clearRect(0,0,W,H);
  // steady-state reference line X* = k1/k2
  const ystar = H - PAD - (k1/k2)/YMAX * (H - 2*PAD);
  ctx.strokeStyle = '#bcd';
  ctx.setLineDash([5,4]); ctx.beginPath(); ctx.moveTo(PAD, ystar); ctx.lineTo(W-4, ystar); ctx.stroke();
  ctx.setLineDash([]);
  ctx.fillStyle = '#7a8aa0'; ctx.font = '11px system-ui';
  ctx.fillText('X* = '+(k1/k2), PAD+4, ystar-4);
  ctx.strokeStyle = '#cdd9e3'; ctx.beginPath();
  ctx.moveTo(PAD, PAD); ctx.lineTo(PAD, H-PAD); ctx.lineTo(W-4, H-PAD); ctx.stroke();
}
const px = t => PAD + (t / T_MAX) * (W - PAD - 4);
const py = x => H - PAD - Math.min(x, YMAX) / YMAX * (H - 2*PAD);

function run(){
  if (timer) return;
  axes();
  let X = 0, t = 0, events = 0;
  const traj = [[0,0]];
  // Run the full SSA up front, storing the step-trajectory.
  while (t < T_MAX) {
    const a1 = k1, a2 = k2 * X, a0 = a1 + a2;
    if (a0 <= 0) break;
    const r1 = Math.random() || 1e-12;
    t += (1 / a0) * Math.log(1 / r1);          // exponential waiting time
    if (t > T_MAX) break;
    const r2 = Math.random();
    if (r2 < a1 / a0) X += 1; else X -= 1;     // creation vs decay
    events++;
    traj.push([t, X]);
  }
  traj.push([T_MAX, X]);
  // Animate drawing the step function.
  let i = 1;
  document.getElementById('run').disabled = true;
  ctx.strokeStyle = '#e63946'; ctx.lineWidth = 2; ctx.beginPath();
  ctx.moveTo(px(traj[0][0]), py(traj[0][1]));
  timer = setInterval(() => {
    const reps = Math.max(1, Math.floor(traj.length / 60));
    for (let r = 0; r < reps && i < traj.length; r++, i++) {
      const [pt, ] = traj[i-1], [ct, cx] = traj[i];
      ctx.lineTo(px(ct), py(traj[i-1][1])); // horizontal hold
      ctx.lineTo(px(ct), py(cx));           // vertical jump
      document.getElementById('count').textContent = cx;
      document.getElementById('time').textContent = ct.toFixed(2);
      document.getElementById('events').textContent = Math.min(i, events);
    }
    ctx.stroke();
    if (i >= traj.length) { clearInterval(timer); timer = null; document.getElementById('run').disabled = false; }
  }, 30);
}
function reset(){ if (timer){ clearInterval(timer); timer = null; } document.getElementById('run').disabled = false;
  axes(); document.getElementById('count').textContent='0'; document.getElementById('time').textContent='0.00'; document.getElementById('events').textContent='0'; }
document.getElementById('run').addEventListener('click', run);
document.getElementById('reset').addEventListener('click', reset);
axes();
Result
Console

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.

How It Really Works

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 aja_j — its current probability per unit time of firing, given how many molecules are around. The total rate is $a0=∑jaja_0 = \sum_j a_j$
  • When? The waiting time until something happens is exponentially distributed with rate a0a_{0}. Draw a uniform random number r1r_{1} and set $τ=1a0ln⁡1r1\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 r2r_{2}, weighted by each reaction's share aj/a0a_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.

Where It Matters

"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.

Conclusion

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.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/gillespie/Content licensed under CC BY-NC 4.0.