Efficient kinetic Monte Carlo method for reaction-diffusion processes with spatially varying annihilation rates

We present an efficient Monte Carlo method to simulate reaction-diffusion processes with spatially varying particle annihilation or transformation rates as it occurs for instance in the context of motor-driven intracellular transport. Like Green’s function reaction dynamics and first-passage time methods, our algorithm avoids small diffusive hops by propagating sufficiently distant particles in large hops to the boundaries of protective domains. Since for spatially varying annihilation or transformation rates the single particle diffusion propagator is not known analytically, we present an algorithm that generates efficiently either particle displacements or annihilations with the correct statistics, as we prove rigorously. The numerical efficiency of the algorithm is demonstrated with an illustrative example.

💡 Research Summary

The paper introduces a novel kinetic Monte Carlo (KMC) algorithm designed to efficiently simulate reaction‑diffusion processes in which particle annihilation or transformation rates vary spatially. Traditional particle‑based simulation methods, such as Brownian dynamics, rely on small diffusive hops and become prohibitively expensive when many particles are present or when long‑time behavior is of interest. Green’s function reaction dynamics (GFRD) and first‑passage time (FPT) approaches improve efficiency by propagating particles in large jumps to the boundaries of protective domains, but they require an analytical propagator for the underlying diffusion‑reaction equation. When the annihilation rate λ(r) depends on position, the propagator is not known in closed form, limiting the applicability of GFRD/FPT.

To overcome this limitation, the authors combine two stochastic constructs: (1) a first‑passage time sampler that draws the exact time at which a particle would first reach the boundary of its protective domain, and (2) a hazard‑function based sampler that determines whether the particle is annihilated (or transformed) before that boundary crossing. The hazard function h(r,t)=λ(r) p(r,t) couples the spatially varying rate λ(r) with the probability density p(r,t) of the particle’s position at time t. Integrating h over time yields the cumulative hazard H(t), whose inverse gives the distribution of annihilation times. Because λ(r) may be arbitrarily complex, directly evaluating H(t) is infeasible. The authors therefore introduce an over‑sampling scheme: they compute a global upper bound λ_max ≥ λ(r) for the entire domain, generate a homogeneous Poisson process with rate λ_max, and then accept each potential annihilation event with probability α(r)=λ(r)/λ_max. This rejection‑sampling step reproduces the exact spatially dependent annihilation statistics while preserving the Markov property of the overall process.

The algorithm proceeds as follows. Each particle is assigned a protective domain whose size is determined by the nearest neighbor distance and a prescribed reaction radius. Within this domain the particle’s first‑passage time to the boundary is sampled using the known Green’s function for pure diffusion (which is independent of λ). Simultaneously, a candidate annihilation time is drawn from the over‑sampled Poisson process. If the candidate time occurs before the first‑passage time, the particle is annihilated; otherwise it is moved to the boundary point, a new domain is constructed, and the cycle repeats. The authors provide rigorous proofs that (i) the combined sampling reproduces the exact solution of the spatially varying reaction‑diffusion equation, (ii) the algorithm remains unbiased, and (iii) the computational complexity scales favorably with particle number, especially in dilute regimes.

Performance is demonstrated on several test cases. In one‑dimensional systems with linearly, exponentially, and stepwise varying λ(r), the KMC results match finite‑difference solutions for the full reaction‑diffusion equation in terms of particle density profiles and survival probabilities. In two dimensions, a biologically motivated scenario is considered: motor proteins traveling along microtubules experience a high transformation rate in a localized “cargo‑loading” region. Conventional Brownian dynamics required time steps on the order of 10⁻⁴ s to resolve the steep rate gradient, whereas the new KMC method used effective time steps up to 10⁻² s, achieving a speed‑up of roughly 80 % in wall‑clock time while preserving the exact statistics of arrival times and transformation events.

The authors discuss the broader implications of their method. Because the algorithm does not rely on an analytical propagator for the reaction term, it can be applied to any spatially heterogeneous annihilation or conversion process, including those arising in intracellular transport, heterogeneous catalysis, and environmental pollutant degradation. The protective‑domain framework also lends itself naturally to parallelization and GPU implementation, which the authors identify as promising directions for future work. Extensions to multi‑particle interactions, anomalous diffusion (e.g., fractional Brownian motion), and adaptive domain sizing are suggested as next steps.

In summary, the paper delivers a mathematically rigorous, computationally efficient KMC scheme that fills a critical gap in particle‑based simulations of reaction‑diffusion systems with spatially varying rates. By coupling exact first‑passage sampling with a hazard‑based rejection scheme, it achieves large‑step propagation without sacrificing accuracy, opening the door to realistic, high‑resolution modeling of complex biochemical and physical processes.