The Master Equation

The master equation governs the probability $P_n(t)$ of a Markovian system occupying discrete state $n$ , balancing the rate at which probability flows in against the rate it leaks out. Its gain–loss structure, the spectrum of its rate matrix, detailed balance, and its Kramers–Moyal reduction to the Fokker–Planck equation are established in the master equation article. This companion takes the equation as given and asks a different, deeply practical question: how do you actually compute with it? The answer runs in two opposite directions, and both are the beating heart of computational statistical mechanics. Read forward, the master equation generates stochastic trajectories. Because a Markov jump process waits an exponentially distributed time in each state — total rate $\alpha_n=\sum_j W_{j\leftarrow n}$ — and then jumps to state $j$ with probability $W_{j\leftarrow n}/\alpha_n$ , one can sample individual realizations exactly. This is the Gillespie stochastic simulation algorithm (SSA), the workhorse of chemical kinetics and single-cell biology, where averaging many such trajectories reproduces the smooth $P_n(t)$ the master equation predicts. Read backward, the logic inverts. Instead of deriving the rates from physics and computing the stationary distribution, one chooses rates that satisfy detailed balance with respect to a prescribed target — the Boltzmann distribution $P_n^{\mathrm{eq}}\propto e^{-\beta E_n}$ — so that the master-equation dynamics is guaranteed to sample it. The Metropolis acceptance rule $\min(1,e^{-\beta\Delta E})$ is exactly such an engineered rate, and Markov chain Monte Carlo is the resulting industry. In both readings the transition rate is the atomic object: forward it is dictated by nature, backward it is a design variable, and the spectral gap of the rate matrix decides — as convergence in one case and mixing time in the other — how long you must wait.

Key takeaways

The master equation (gain–loss balance, rate matrix, detailed balance, Fokker–Planck limit) is derived and owned in the master-equation article; this article is the computational companion.
Forward use: the Gillespie SSA samples exact trajectories — draw an exponential waiting time with rate $\alpha_n=\sum_j W_{j\leftarrow n}$ , then jump to $j$ with probability $W_{j\leftarrow n}/\alpha_n$ .
Averaging many SSA trajectories reproduces the master equation's $P_n(t)$ — the ensemble of jump paths is the ODE.
Backward use: Markov chain Monte Carlo engineers rates satisfying detailed balance for a target $P^{\mathrm{eq}}$ , so the dynamics samples it. Metropolis: accept with $\min(1,P^{\mathrm{eq}}_m/P^{\mathrm{eq}}_n)=\min(1,e^{-\beta\Delta E})$ .
Metropolis satisfies detailed balance for any symmetric proposal, guaranteeing the correct stationary distribution.
The spectral gap sets the SSA relaxation and the MCMC autocorrelation/mixing time; it closing is critical slowing down.
Key idea: the transition rate is the atomic object — read forward it dictates trajectories (Gillespie), read backward it is a design knob tuned by detailed balance to sample a target (Metropolis). Same master equation, two directions.

The core insight

The transition rate is the atomic object. Read forward, the master equation dictates exact stochastic trajectories (the Gillespie SSA); read backward, the rate becomes a design knob, engineered via detailed balance so the dynamics samples a prescribed target (Metropolis Monte Carlo). The master equation is both a law of nature and an algorithm.

Where this lives today

The two readings of the master equation power much of computational science. The Gillespie SSA is the standard tool of stochastic chemical kinetics and systems biology, where small molecule numbers make deterministic rate equations meaningless — gene-expression bursts, single-cell variability, and viral dynamics are all simulated trajectory-by-trajectory, with modern accelerations (tau-leaping, next-reaction methods) trading exactness for speed on stiff networks. In the opposite direction, Markov chain Monte Carlo — Metropolis–Hastings and its descendants (heat-bath, Wolff and Swendsen–Wang cluster algorithms, Hamiltonian Monte Carlo, parallel tempering) — is the engine of equilibrium statistical mechanics, lattice field theory, and Bayesian inference, sampling Gibbs measures over spaces far too large to enumerate. In both worlds the practical bottleneck is the spectral gap: critical slowing down near phase transitions is the master equation's slowest eigenvalue creeping toward zero, and an entire craft of algorithm design exists to engineer rates — the same rates that appear in the equation — that relax faster.

The physical picture

Start from the master equation as a settled fact: probability is a conserved fluid on a network of states, each link carrying a one-way current $W_{m\to n}P_m$ , and the whole evolution is fixed by the rate matrix. That picture and its spectral consequences belong to the master-equation article. The question here is operational — a single realization of the process, not the probability cloud. Watch one system. It sits in state $n$ for a random time, then jumps. Markovianity forces that residence time to be exponentially distributed with rate $\alpha_n$ equal to the sum of all rates out of $n$ — the total escape rate. When it jumps, the destination is chosen among the open channels in proportion to their individual rates. That is the entire content of the Gillespie algorithm: sample when from an exponential, sample where from the normalized rates. Run it many times and the histogram of occupancies at each instant reconstructs the smooth $P_n(t)$ the master equation would have integrated — the punchline figure shows the jumpy sample paths converging, on average, onto the deterministic curve. Now invert. Suppose you do not know the rates but you do know the distribution you want to sample — say the Boltzmann weight over a huge configuration space you cannot enumerate. Detailed balance, $W_{n\to m}P^{\mathrm{eq}}_n=W_{m\to n}P^{\mathrm{eq}}_m$ , is a sufficient condition guaranteeing that $P^{\mathrm{eq}}$ is stationary. So engineer rates that satisfy it. The Metropolis choice — propose a move symmetrically, then accept it with probability $\min(1,e^{-\beta\Delta E})$ — does exactly that, for any energy landscape, using only energy differences. The master equation you have built then relaxes to the Boltzmann distribution, and sampling its trajectory is sampling the ensemble. How long you must run is set, once again, by the rate matrix's spectral gap.

Figure

Punchline: Gillespie SSA jump trajectories scatter around the smooth master-equation solution; their ensemble mean converges to the ODE curve

Figure

The Metropolis rule: accept with min(1, e^{−βΔE}) — a detailed-balance-engineered transition rate that samples the Boltzmann target

The core equations

For orientation the master equation is quoted here in the form owned and derived by the master-equation article; this article does not re-derive it but uses it as the object to simulate and to engineer.

\frac{\mathrm{d}P_n(t)}{\mathrm{d}t}=\sum_{m\,(\neq n)}\big[\,W_{n\leftarrow m}\,P_m(t)-W_{m\leftarrow n}\,P_n(t)\,\big] \qquad (1)

where: $P_n(t)$ — probability of state $n$ ; $W_{n\leftarrow m}$ — transition rate from $m$ to $n$ (probability per unit time, $\ge0$ ). The gain–loss structure, probability conservation, the rate-matrix spectrum, and detailed balance are established in the master-equation article. The forward computational content is the residence-time law of a Markov jump process. The total escape rate from state $n$ sets an exponential holding time, and the destination is drawn from the normalized outgoing rates.

\alpha_n=\sum_{j\neq n}W_{j\leftarrow n},\qquad \tau\sim\mathrm{Exp}(\alpha_n),\qquad \Pr(n\to j)=\frac{W_{j\leftarrow n}}{\alpha_n} \qquad (2)

where: $\alpha_n$ — total exit rate from $n$ ; $\tau$ — random waiting time, exponentially distributed with mean $1/\alpha_n$ ; $\Pr(n\to j)$ — probability that the jump lands in $j$ . Equations (2) are the exact Gillespie SSA update. The backward content is detailed balance imposed as a design constraint: choose rates so a prescribed $P^{\mathrm{eq}}$ is the stationary distribution. Detailed balance is sufficient for stationarity (proved in the master-equation article) and, crucially, depends only on rate ratios.

\frac{W_{m\leftarrow n}}{W_{n\leftarrow m}}=\frac{P^{\mathrm{eq}}_m}{P^{\mathrm{eq}}_n}=e^{-\beta(E_m-E_n)} \qquad (3)

where: $P^{\mathrm{eq}}_n\propto e^{-\beta E_n}$ — the Boltzmann target; only energy differences $E_m-E_n=\Delta E$ enter, so the intractable normalization (partition function) cancels. The Metropolis–Hastings realization splits each rate into a symmetric proposal $g_{nm}=g_{mn}$ times an acceptance probability, and the acceptance is chosen to enforce (3).

W_{m\leftarrow n}=g_{nm}\,A_{n\to m},\qquad A_{n\to m}=\min\!\Big(1,\frac{P^{\mathrm{eq}}_m}{P^{\mathrm{eq}}_n}\Big)=\min\!\big(1,e^{-\beta\Delta E}\big) \qquad (4)

where: $g_{nm}$ — symmetric proposal probability; $A_{n\to m}$ — Metropolis acceptance. Downhill moves ( $\Delta E\le0$ ) are always accepted; uphill moves are accepted with $e^{-\beta\Delta E}$ . The cost of sampling is set — as for any master equation — by the spectral gap of the rate matrix, which controls both the SSA relaxation and the MCMC autocorrelation time.

\tau_{\mathrm{mix}}\sim\frac{1}{|\lambda_1|},\qquad |\lambda_1|=\text{smallest nonzero }|\mathrm{Re}\,\lambda|\ \text{of the rate matrix} \qquad (5)

where: $\lambda_1$ — slowest relaxation eigenvalue (spectral gap); $\tau_{\mathrm{mix}}$ — mixing/autocorrelation time. Its closing near a critical point is critical slowing down.

Derivation AThe Gillespie stochastic simulation algorithm reproduces the master-equation solution

The goal is to show that averaging exact stochastic trajectories generated by the Gillespie SSA reproduces the probability distribution $P_n(t)$ that the master equation integrates — that the ensemble of jump paths is the ODE. The master equation itself, its gain–loss form and rate-matrix spectrum, are owned by the master-equation article and recalled here by name only; what is derived is the sampling scheme and its consistency. The strategy is to build the SSA update from the exponential residence law and normalized jump probabilities, then argue that the empirical occupation frequency converges to $P_n(t)$ by the law of large numbers over independent trajectories. Keep in view that the SSA is exact — no time-step discretization error — because it samples the true continuous-time jump statistics.

Step 1. In state $n$ , the times to each individual channel $n\to j$ are independent exponentials with rates $W_{j\leftarrow n}$ . The first to fire determines both the waiting time and the destination; the minimum of independent exponentials is exponential with rate $\alpha_n=\sum_j W_{j\leftarrow n}$ , and channel $j$ wins with probability $W_{j\leftarrow n}/\alpha_n$ . Why: competing Poisson clocks combine into a single clock at the summed rate, and the winner is picked in proportion to its own rate — this is the exact statistics the master equation encodes for a single realization.

Step 2. The SSA update follows: draw $\tau=-\ln(u_1)/\alpha_n$ with $u_1\sim U(0,1)$ (inverse-transform sampling of the exponential), advance time by $\tau$ , then select the destination $j$ by drawing $u_2\sim U(0,1)$ against the cumulative normalized rates. Why: inverse-transform sampling turns a uniform variate into an exponential holding time exactly, and cumulative selection realizes $\Pr(n\to j)=W_{j\leftarrow n}/\alpha_n$ .

Step 3. Over an ensemble of $N$ independent trajectories, the fraction occupying state $n$ at time $t$ , $\hat P_n(t)=\tfrac1N\sum_r \mathbb{1}[X_r(t)=n]$ , is an unbiased estimator of $P_n(t)$ with variance $O(1/N)$ . Why: each trajectory is a draw from the true process, so by the law of large numbers the empirical frequency converges to the exact master-equation probability, statistical error shrinking as $1/\sqrt N$ .

Step 4. For the two-state telegraph process this is checkable in closed form: the ensemble-mean occupation relaxes as $P_1(t)=P_1^{\mathrm{eq}}+[P_1(0)-P_1^{\mathrm{eq}}]e^{-\alpha t}$ with $\alpha=k_\uparrow+k_\downarrow$ and $P_1^{\mathrm{eq}}=k_\downarrow/(k_\uparrow+k_\downarrow)$ — precisely the master-equation solution. Why: the SSA cannot produce anything but the master equation's own dynamics, so its ensemble average must reproduce the analytic curve. *Verified: verify_aq-sm-master-equation-sm11.py runs the SSA for a two-state process over many trajectories and confirms the ensemble mean matches $P_1(t)=P_1^{\mathrm{eq}}+[\cdots]e^{-\alpha t}$ and the equilibrium occupation $k_\downarrow/(k_\uparrow+k_\downarrow)$ within Monte-Carlo error.*

Go deeper · Derivation B · 2 minMetropolis acceptance satisfies detailed balance for any symmetric proposal

Derivation BMetropolis acceptance satisfies detailed balance for any symmetric proposal

The goal is to prove that the Metropolis acceptance rule $A_{n\to m}=\min(1,P^{\mathrm{eq}}_m/P^{\mathrm{eq}}_n)$ , combined with any symmetric proposal, produces a master-equation dynamics whose stationary distribution is exactly the target $P^{\mathrm{eq}}$ — this is the inverse use of the master equation, engineering rates rather than deriving them. The fact that detailed balance implies stationarity is owned by the master-equation article and used here as a lemma; what is derived is that the Metropolis rates satisfy detailed balance. The strategy is to write the effective transition rate as proposal times acceptance and check the detailed-balance ratio directly, treating the two cases $P^{\mathrm{eq}}_m\gtrless P^{\mathrm{eq}}_n$ symmetrically. The elegance to notice: only the ratio $P^{\mathrm{eq}}_m/P^{\mathrm{eq}}_n$ appears, so the unknown normalization cancels and one samples $e^{-\beta E}$ without ever computing the partition function.

Step 1. Write the transition probability (per sweep) as $W_{m\leftarrow n}=g_{nm}A_{n\to m}$ with a symmetric proposal $g_{nm}=g_{mn}$ . The detailed-balance condition to verify is $W_{m\leftarrow n}P^{\mathrm{eq}}_n=W_{n\leftarrow m}P^{\mathrm{eq}}_m$ , i.e. $A_{n\to m}P^{\mathrm{eq}}_n=A_{m\to n}P^{\mathrm{eq}}_m$ after cancelling the symmetric $g$ . Why: symmetric proposals drop out of the balance condition, isolating the requirement on the acceptance probabilities alone.

Step 2. Suppose $P^{\mathrm{eq}}_m\ge P^{\mathrm{eq}}_n$ (downhill in free energy toward $m$ ). Then $A_{n\to m}=\min(1,P^{\mathrm{eq}}_m/P^{\mathrm{eq}}_n)=1$ and $A_{m\to n}=\min(1,P^{\mathrm{eq}}_n/P^{\mathrm{eq}}_m)=P^{\mathrm{eq}}_n/P^{\mathrm{eq}}_m$ . Why: the move into the more probable state is always accepted, and the reverse move is accepted with the ratio, exactly compensating.

Step 3. Substitute: $A_{n\to m}P^{\mathrm{eq}}_n=1\cdot P^{\mathrm{eq}}_n=P^{\mathrm{eq}}_n$ and $A_{m\to n}P^{\mathrm{eq}}_m=(P^{\mathrm{eq}}_n/P^{\mathrm{eq}}_m)P^{\mathrm{eq}}_m=P^{\mathrm{eq}}_n$ . The two are equal, so detailed balance holds; the case $P^{\mathrm{eq}}_m<P^{\mathrm{eq}}_n$ is identical by relabelling. Why: the $\min$ construction is built precisely so the accepted forward and reverse fluxes match, independent of the landscape.

Step 4. By the lemma (detailed balance ⟹ stationarity, from the master-equation article), $P^{\mathrm{eq}}$ is the stationary distribution of the engineered chain; with an irreducible aperiodic proposal it is the unique attractor, so the dynamics samples $P^{\mathrm{eq}}\propto e^{-\beta E}$ . Only $\Delta E=E_m-E_n$ enters, so the partition function never appears. Why: this is what makes Monte Carlo work at all — sampling a Gibbs measure over a space too large to enumerate, using only local energy differences. *Verified: verify_aq-sm-master-equation-sm11.py checks $A_{n\to m}P^{\mathrm{eq}}_n=A_{m\to n}P^{\mathrm{eq}}_m$ symbolically for both orderings, and runs a Metropolis chain on a two-level system confirming $p_{\mathrm{up}}/p_{\mathrm{down}}=e^{-\beta\Delta E}$ within Monte-Carlo error.*

Limiting cases

The computational readings of the master equation across regimes:

Case	Result	Meaning
Single trajectory (forward)	$\tau\sim\mathrm{Exp}(\alpha_n),\ \Pr(n\to j)=W_{j\leftarrow n}/\alpha_n$	Exact Gillespie SSA update; no time-step error, samples the true continuous-time jump process.
Ensemble of trajectories	$\hat P_n(t)\to P_n(t)$ as $N\to\infty$	The empirical occupation frequency converges to the master-equation solution, error $\sim1/\sqrt N$ .
Engineered rates (backward)	$A_{n\to m}=\min(1,e^{-\beta\Delta E})$	Metropolis: detailed balance for any symmetric proposal ⟹ samples $P^{\mathrm{eq}}\propto e^{-\beta E}$ using only $\Delta E$ .
Large-copy-number limit	$\dot{\mathbf x}=\mathbf f(\mathbf x)$ (rate equations)	Averaging the chemical master equation gives deterministic kinetics; SSA fluctuations $\sim1/\sqrt{\Omega}$ (recall: Kramers–Moyal, owned by the master-equation article).
Critical / near-degenerate	$\tau_{\mathrm{mix}}\sim1/\|\lambda_1\|\to\infty$	Spectral gap closes: critical slowing down; cluster algorithms enlarge the gap.

Worked examples

Example 1 — Gillespie SSA for a two-state process reproduces the equilibrium occupation. Take a telegraph process with $k_\uparrow=3.0\,\mathrm{s^{-1}}$ (down→up) and $k_\downarrow=1.0\,\mathrm{s^{-1}}$ (up→down). The master-equation equilibrium is $P^{\mathrm{eq}}_{\mathrm{up}}=k_\uparrow/(k_\uparrow+k_\downarrow)$ .

Step 1. Compute the predicted occupation: $P^{\mathrm{eq}}_{\mathrm{up}}=3.0/(3.0+1.0)=\mathbf{0.75}$ , and the relaxation rate $\alpha=k_\uparrow+k_\downarrow=\mathbf{4.0\ s^{-1}}$ , so the relaxation time is $1/\alpha=\mathbf{0.25\ s}$ .

Step 2. Run the SSA: in the down state the holding time is $\mathrm{Exp}(k_\uparrow)$ with mean $1/3\,\mathrm s$ ; in the up state it is $\mathrm{Exp}(k_\downarrow)$ with mean $1\,\mathrm s$ . The long-run fraction of time spent up equals the ratio of mean up-time to total cycle time: $\frac{1/k_\downarrow}{1/k_\downarrow+1/k_\uparrow}=\frac{1}{1+k_\downarrow/k_\uparrow}=\frac{1}{1+1/3}=\mathbf{0.75}$ . The trajectory-time fraction matches the master-equation equilibrium exactly.

Step 3. Averaging $N=20000$ SSA trajectories at $t=2\,\mathrm s$ (eight relaxation times, fully equilibrated) gives $\hat P_{\mathrm{up}}\approx0.75$ within Monte-Carlo error $\sim1/\sqrt N\approx0.007$ . The ensemble of jump paths reproduces the smooth $P_n(t)$ .

Example 2 — Metropolis sampling of a two-level Boltzmann target. Two states with energy gap $\Delta E=E_1-E_0>0$ at inverse temperature $\beta$ . The exact Boltzmann ratio is $p_1/p_0=e^{-\beta\Delta E}$ ; verify that Metropolis reproduces it.

Step 1. Take $\beta\Delta E=2.0$ , so the target ratio is $p_1/p_0=e^{-2.0}=\mathbf{0.1353}$ , giving $p_1=e^{-2}/(1+e^{-2})=\mathbf{0.1192}$ and $p_0=\mathbf{0.8808}$ .

Step 2. Metropolis with a symmetric proposal (flip the state): from $0\to1$ (uphill) accept with $\min(1,e^{-\beta\Delta E})=e^{-2}=0.1353$ ; from $1\to0$ (downhill) accept with $\min(1,e^{+\beta\Delta E})=1$ . Check detailed balance: $A_{0\to1}p_0=0.1353\times0.8808=0.1192$ and $A_{1\to0}p_1=1\times0.1192=0.1192$ — equal, detailed balance holds.

Step 3. Running the chain for $N=200000$ sweeps yields an occupation ratio $\hat p_1/\hat p_0\approx e^{-2}=\mathbf{0.135}$ within Monte-Carlo error. The engineered master-equation dynamics samples the Boltzmann distribution using only the energy difference $\Delta E$ .

Origins · 1 mintransition rate and how this was discovered

transition rate

the atomic object of jump-process kinetics, 1940s–1970s · A transition rate

W_{n\leftarrow m}

is the probability per unit time that a Markov system jumps from state

m

to state

n

. It is the fundamental input of the master equation and the pivot of both its computational uses: read forward, the rates dictate the exponential holding times and jump probabilities of the Gillespie stochastic simulation algorithm; read backward, they become design variables tuned by detailed balance so the dynamics samples a target Boltzmann distribution, as in the 1953 Metropolis algorithm. Whether dictated by scattering cross sections and reaction propensities or engineered for Monte Carlo, the transition rate — carrying units of inverse time and constrained only in its ratios by detailed balance — is where physics and computation meet.

History

1906

— Andrey Markov defines the memoryless chains whose continuous-time limit is the master equation

1928

— Sydney Chapman and 1931 Andrei Kolmogorov give the composition law and analytic theory behind the rate description

1940

— Hendrik Kramers casts escape kinetics as a rate process, sharpening the notion of a transition rate

1953

— Metropolis, Rosenbluth, Rosenbluth, Teller and Teller devise the Monte Carlo acceptance rule that engineers detailed-balanced rates (J. Chem. Phys. 21, 1087)

1970

— W. Keith Hastings generalizes Metropolis to asymmetric proposals, founding Markov chain Monte Carlo as a general method

1976

— Daniel Gillespie formulates the exact stochastic simulation algorithm for chemical master equations (J. Comput. Phys. 22, 403)

Going further

References. Pathria & Beale, Statistical Mechanics, 3rd ed., Ch. 15 (Fluctuations and stochastic processes) · Kardar, Statistical Physics of Particles, Ch. 3 · Sethna, Statistical Mechanics: Entropy, Order Parameters, and Complexity, Ch. 8 · van Kampen, Stochastic Processes in Physics and Chemistry, Ch. V · Gardiner, Handbook of Stochastic Methods, Ch. 3 · Newman & Barkema, Monte Carlo Methods in Statistical Physics · Gillespie, J. Comput. Phys. 22, 403 (1976) and J. Phys. Chem. 81, 2340 (1977) · Metropolis et al., J. Chem. Phys. 21, 1087 (1953)

Provenance. The master equation, its gain–loss form, rate-matrix spectrum, detailed balance, and Kramers–Moyal reduction to Fokker–Planck are derived and owned by the master-equation article and are recalled here by name, not re-derived. A reviewer should re-derive the computational content: (1) the Gillespie update from the exponential residence law and normalized jump probabilities, confirming the ensemble mean reproduces the master-equation $P_n(t)$ for the two-state process; (2) that the Metropolis acceptance $\min(1,e^{-\beta\Delta E})$ satisfies detailed balance for any symmetric proposal, so the chain samples $P^{\mathrm{eq}}\propto e^{-\beta E}$ ; (3) that only energy differences enter. Cross-references: Gardiner Ch. 3, van Kampen Ch. V, Newman & Barkema Chs. 2–4.

Builds on: Markov Processes and the Markov Property · Chapman–Kolmogorov Equation · The Master Equation · Detailed Balance and Microscopic Reversibility. Leads to: Fokker–Planck Equation · Metropolis–Hastings Algorithm · Monte Carlo Method · Kinetic Theory and the Boltzmann Equation.

Isn't this just the same as the master equation article — why a separate page?

No — the two are deliberately complementary and de-overlapped. The master-equation article owns the physics of the equation: its gain–loss derivation from Chapman–Kolmogorov, the rate-matrix spectrum (zero eigenvalue = stationary state, others = relaxation), detailed balance, the H-theorem, and the Kramers–Moyal reduction to Fokker–Planck. This article owns the computation: how to sample exact trajectories forward (Gillespie SSA) and how to engineer rates backward so the dynamics samples a prescribed target (Metropolis Monte Carlo). It recalls the equation's structure by name and never re-derives it.

Why is a single Gillespie trajectory noisy if the master equation is deterministic?

Because the master equation describes the probability distribution over states, an average object, while a Gillespie run is one realization — an actual sequence of jumps. Any single system executes a jagged, random path; the smooth $P_n(t)$ emerges only when you average over many independent realizations. This is the same relationship as between a single die roll and the probability distribution of outcomes: the distribution is deterministic, each roll is not. The SSA's job is to produce statistically exact samples, so their ensemble average reconstructs the master-equation curve with error $\sim1/\sqrt N$ .

How does the Gillespie algorithm avoid time-step error?

It does not discretize time at all. Instead of asking "what happens in the next $\Delta t$ ?", it asks "when is the next event, and which event is it?", and answers both exactly by sampling the true continuous-time statistics: the next-event time is an exponential with rate equal to the total propensity, and the event is chosen in proportion to individual rates. Because these are the exact laws of the Markov jump process, the trajectory has no discretization bias whatsoever — the only error is statistical, from finite sampling. Fixed-step methods (like naive Euler on the master ODE, or $\tau$ -leaping) trade this exactness for speed.

Why must the Metropolis proposal be symmetric, and what if it isn't?

The plain Metropolis acceptance $\min(1,P^{\mathrm{eq}}_m/P^{\mathrm{eq}}_n)$ assumes the proposal probabilities are symmetric, $g_{nm}=g_{mn}$ , so they cancel in the detailed-balance ratio, leaving a condition on the acceptance alone. If the proposal is asymmetric, one must correct for it — this is exactly Hastings' 1970 generalization, with acceptance $\min\!\big(1,\frac{P^{\mathrm{eq}}_m g_{mn}}{P^{\mathrm{eq}}_n g_{nm}}\big)$ . The extra ratio $g_{mn}/g_{nm}$ restores detailed balance for biased proposals, and Metropolis–Hastings is the resulting general algorithm.

What sets how long a Monte Carlo simulation must run?

The spectral gap of the effective rate matrix — the same object that governs relaxation in the forward reading. Writing the chain's evolution in its eigenmodes, every non-stationary mode decays at a rate given by an eigenvalue, and the slowest (the smallest nonzero $|\mathrm{Re}\,\lambda|$ ) sets the autocorrelation time $\tau_{\mathrm{mix}}\sim1/|\lambda_1|$ . Independent samples are separated by roughly $\tau_{\mathrm{mix}}$ sweeps. Near a continuous phase transition this eigenvalue approaches the stationary one, the gap closes, and autocorrelation times diverge — critical slowing down — which is why cluster algorithms, engineered to enlarge the gap, are essential near criticality.

Can the Gillespie SSA and Metropolis be seen as the same object?

Yes, and that is the point of this article. Both are master-equation dynamics; they differ only in whether the rates are given or chosen. In the Gillespie case the rates come from physics (reaction propensities, scattering cross sections) and you sample the resulting trajectories to learn the kinetics. In the Metropolis case you fix the stationary distribution you want and solve detailed balance for admissible rates, then run the dynamics to sample it. Forward: rates → distribution. Backward: distribution → rates. The transition rate is the shared atomic object, which is why it is the portrait of this article.