= Solution
An exact <Gillespie algorithm> for one particle is:
* Set the current well $k$ and time $t=0$.
* If $k$ is absorbing, stop.
* Set $a=k_++k_-$ and draw $\Delta t=-\log U_1/a$.
* Move to $k+1$ if $U_2<k_+/a$, and otherwise to $k-1$.
* Set $t\leftarrow t+\Delta t$ and repeat, stopping if the jump crosses an absorbing end.
For $N\ll K$, simulate particle identities independently and maintain a priority queue of their next event times. For $N\gg K$, store occupation numbers $n_k$ and use aggregate event rates $n_kk_+$ and $n_kk_-$ for each well; one population-level Gillespie event then decrements one $n_k$ and increments its neighbor. This replaces work proportional to particle number by work proportional to the number of occupied wells.
Back to article page
