Turning Reaction Rates Into Gillespie Propensities

From Reaction Rates to Event Rates

A reaction-rate constant is not, by itself, the rate at which a reaction fires in a stochastic simulation. Its meaning depends on the rate-law convention, the current molecule counts, and—when the rate law is written in concentrations—the volume. The propensity is the state-dependent event rate needed to describe reaction events. A second ingredient is the state-change vector, which records how one firing changes the molecule counts.

Let \((\Omega,\mathcal{F},\mathbb{P})\) be a probability space and let \(X_t\) be a continuous-time reaction network with state \(X_t=x\in S\), where \(S\) is a subset of nonnegative integer count vectors. Write \(x_i\) for the current count of species \(i\), and let \((\mathcal{F}_t)_{t\ge0}\) be the history of the process through time \(t\). Each reaction channel \(r\) has a propensity \(a_r(x)\), with units of events per unit time. Its local meaning is that, conditional on \(X_t=x\), the probability of one firing of channel \(r\) in a short interval of length \(h\) is \(a_r(x)h+o(h)\) as \(h\) decreases to zero; the probability of two or more firings is \(o(h)\). This is a local rate description, not a claim that \(a_r(x)h\) is an exact event probability for arbitrary \(h\).

Write a reaction as a list of reactants and products. If \(\alpha_{ir}\) molecules of species \(i\) are consumed and \(\beta_{ir}\) are produced in one firing of channel \(r\), its state-change vector is \(\nu_r=\beta_r-\alpha_r\). Thus a firing takes \(x\) to \(x+\nu_r\), provided the reactants are available. The propensity describes how quickly the channel fires; the vector describes what the firing does.

These formulas apply when the specified reactants are present. If a reaction requires more molecules than the state contains, its propensity is zero. In the formulas, \(k\) for a unimolecular reaction has units of inverse time; \(k\) for the stated bimolecular concentration law has units of volume per molecule per time, if counts are literal molecule counts. Concentration and volume units must be consistent. For example, a rate constant expressed per mole uses concentrations in moles per volume, with counts converted accordingly.

Why the Count Formulas Have These Factors

The count factors come from counting available reactant molecules or combinations. They also explain why identical-reactant formulas require care: the number of distinct pairs is not \(x_A^2\).

Proof. If there are \(x_A\) molecules, a unimolecular event can occur through any one of them. Adding their rates gives \(x_Ak=kx_A\). For distinct species, each of the \(x_A\) molecules of \(A\) can be paired with each of the \(x_B\) molecules of \(B\), so there are \(x_Ax_B\) available pairs. Summing the assigned rate \(k/V\) over these pairs gives \(kx_Ax_B/V\). For identical reactants, an event uses two different molecules, and the number of unordered pairs is \(\binom{x_A}{2}=x_A(x_A-1)/2\). Multiplying that count by the assigned rate \(2k/V\) per pair gives \((2k/V)x_A(x_A-1)/2=kx_A(x_A-1)/V\). If \(x_A=0\) or \(1\), the pair count and propensity are both zero. \(\square\)

For the distinct-species formula, dividing the propensity by volume gives \(a(x)/V=k(x_A/V)(x_B/V)=k[A][B]\), exactly the stated concentration rate law. For identical reactants under the no-\(1/2\) convention, the finite-count rate per volume is instead

This is an exact finite-volume identity. At fixed concentration as volume grows, the correction \(k[A]/V\) tends to zero, so the rate per volume converges to \(k[A]^2\). It is not generally equal to \(k[A]^2\) at finite volume. If a source defines the macroscopic identical-reactant law as \(v=(k/2)[A]^2\), the corresponding pair-count formula under the same counting convention is \(a(x)=k x_A(x_A-1)/(2V)\). Always check whether the factor \(1/2\) is included in the rate-law definition before using a formula.

The state-change vector is determined separately from the rate constant. If a reaction consumes one \(A\) and produces one \(B\), then its \(A,B\) changes are \(-1,+1\), regardless of how large or small its propensity is. More generally, \(x+\nu_r\) must have nonnegative entries whenever \(a_r(x)>0\).

Worked Conversions

Common Errors and a Reliable Checklist

A frequent error is to use \(kx_A^2/V\) for an identical-reactant channel with the no-\(1/2\) macroscopic convention. For \(k=2,\ V=3,\ x_A=4\), this incorrect expression gives \(2\cdot16/3=32/3\) events per minute. The correct propensity counts distinct unordered pairs:

The correct finite-count rate per volume is \(8/3\). Since \([A]=4/3\), the macroscopic expression is \(2(4/3)^2=32/9\), and the correction is \(2(4/3)/3=8/9\); subtracting gives \(32/9-8/9=24/9=8/3\). The error uses ordered/self-pairs rather than the available unordered pairs and also mistakes the macroscopic law for an exact finite-count identity.

For each reaction channel, a reliable conversion proceeds in this order:

1. Write the reactants, products, and the precise convention for the rate law, including any factor \(1/2\).
2. Check whether the law is in counts or concentrations; for a concentration law, include the appropriate volume conversion.
3. Count available reactants: one molecule for a unimolecular event, \(x_Ax_B\) pairs for distinct species, and \(x_A(x_A-1)/2\) unordered pairs for identical species.
4. Compute the propensity and verify its units are events per time and its value is nonnegative.
5. Write the state-change vector as products minus reactants, then verify that a permitted firing leaves nonnegative counts.

Propensities are inputs describing the rates of the individual reaction channels at the current state. Distinct channels remain distinct even if they have the same state-change vector: their separate rate constants and reactant requirements still matter. The subsequent simulation method uses these inputs to determine event times and channel choices; the conversion itself does not select an event or generate a waiting time.