 Solution

An exact Gillespie algorithm for one particle is:

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_k{k}_{+}$ and $n_k{k}_{-}$ 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.