Kinetic Monte Carlo against the master equation

The exact solution in the reference, posts/refs/kinetic-monte-carlo-against-the-master-equation.py, is the matrix exponential $p(t)=e^{Wt}p(0)$ from scipy.linalg.expm. The program solves the same equation a different way, by Runge-Kutta steps, and the first test compares the two. The same file runs a separate Gillespie code 200 times with numpy's random numbers, and its run-to-run spread matches the error bars used here.

Rates and the master equation

Many processes are jumps between a few states at known rates: a molecule changing shape, an atom hopping between sites on a surface, a gene switching on and off. A rate $k_{i\to j}$ means that in a short time $dt$ a system in state $i$ jumps to $j$ with probability $k_{i\to j}dt$, independently of how long it has already been in $i$.

Here the states are A, B and C in a chain,

and every molecule starts in A. The probability $p_i(t)$ of being in state $i$ gains from jumps into $i$ and loses from jumps out of it. That balance is the master equation written above, a linear system $dp/dt=Wp$ with a 3 by 3 matrix $W$. Its long-time answer needs no solving: each forward and backward pair balances when $p_B/p_A=2/0.5=4$ and $p_C/p_B=1/0.25=4$, so $p\to (1,4,16)/21$.

The program integrates $dp/dt=Wp$ with the classical fourth-order Runge-Kutta method at step $10^{-3}$. The fastest decay rate of $W$ is 2.82, so a step of $10^{-3}$ makes the method's error at $t=1$ between $3\times 10^{-14}$ and $3\times 10^{-13}$.

Kinetic Monte Carlo

The master equation tracks probabilities. Kinetic Monte Carlo (KMC) follows one molecule instead, and gets the probabilities by counting many molecules. The version here is Gillespie's algorithm (1977). In state $i$ the total rate of leaving is $a={\sum }_j{k}_{i\to j}$. Since each jump happens at a constant rate, the waiting time until the next one is exponential with mean $1/a$, drawn as $-\ln (u_1)/a$ from a uniform random number $u_1$. Which jump it is follows from a second draw: jump to $j$ with probability $k_{i\to j}/a$. No time step is chosen, and none of the jumps is approximated.

The cell runs 10000 molecules to $t=5$, records which state each is in at the 25 sample times, and notes when it first reaches C. Three trajectories are kept for the picture.

How close should they be?

At a fixed time each molecule is in state $i$ or not, independently of the others, so the count in $i$ is binomial and the fraction has standard error $\sqrt{p(1-p)/N}$, at most 0.005 for $N=10000$. The reference's 200 repeated runs spread by that amount. There are 75 comparisons (3 states at 25 times), so the test allows four standard errors: one comparison exceeds that with probability $6.3\times 10^{-5}$, and any of the 75 with at most 0.5%.

The mean first-passage time from A to C is the average time until a molecule starting in A first reaches C. Treating C as a trap, it solves ${\tau }_A=\frac{1}{2}+{\tau }_B$ and ${\tau }_B=\frac{1}{1.5}+\frac{0.5}{1.5}{\tau }_A$ (wait in the state, then take the next jump), so ${\tau }_A=1.75$. The measured mean has a standard error from the sample spread of the 10000 first-passage times.

The bars on the points are four standard errors, the test's tolerance; the largest deviation of the 75 is printed as max_z, in standard errors: 2.6. The measured first-passage mean, 1.721, is 2.05 standard errors below 1.75; a normal variable lands that far out about one time in 25. B is an intermediate: it fills as A empties and drains into C, and its peak of 0.460 at $t=0.70$ is above its final share of $4/21=0.19$. By $t=5$ the slowest part of the approach, which decays as $e^{-0.931t}$ (the smaller nonzero eigenvalue of $W$), has not died out, so $p_C$ is still 0.011 short of $16/21$.