A First-Passage Kinetic Monte Carlo Algorithm for Complex Diffusion-Reaction Systems

We develop an asynchronous event-driven First-Passage Kinetic Monte Carlo (FPKMC) algorithm for continuous time and space systems involving multiple diffusing and reacting species of spherical particles in two and three dimensions. The FPKMC algorithm presented here is based on the method introduced in [Phys. Rev. Lett., 97:230602, 2006] and is implemented in a robust and flexible framework. Unlike standard KMC algorithms such as the n-fold algorithm, FPKMC is most efficient at low densities where it replaces the many small hops needed for reactants to find each other with large first-passage hops sampled from exact time-dependent Green's functions, without sacrificing accuracy. We describe in detail the key components of the algorithm, including the event-loop and the sampling of first-passage probability distributions, and demonstrate the accuracy of the new method. We apply the FPKMC algorithm to the challenging problem of simulation of long-term irradiation of metals, relevant to the performance and aging of nuclear materials in current and future nuclear power plants. The problem of radiation damage spans many decades of time-scales, from picosecond spikes caused by primary cascades, to years of slow damage annealing and microstructure evolution. Our implementation of the FPKMC algorithm has been able to simulate the irradiation of a metal sample for durations that are orders of magnitude longer than any previous simulations using the standard Object KMC or more recent asynchronous algorithms.Comment: See also arXiv:0905.357

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.Comment: 13 pages, 5 figure

Efficient Reactive Brownian Dynamics

We develop a Split Reactive Brownian Dynamics (SRBD) algorithm for particle simulations of reaction-diffusion systems based on the Doi or volume reactivity model, in which pairs of particles react with a specified Poisson rate if they are closer than a chosen reactive distance. In our Doi model, we ensure that the microscopic reaction rules for various association and disassociation reactions are consistent with detailed balance (time reversibility) at thermodynamic equilibrium. The SRBD algorithm uses Strang splitting in time to separate reaction and diffusion, and solves both the diffusion-only and reaction-only subproblems exactly, even at high packing densities. To efficiently process reactions without uncontrolled approximations, SRBD employs an event-driven algorithm that processes reactions in a time-ordered sequence over the duration of the time step. A grid of cells with size larger than all of the reactive distances is used to schedule and process the reactions, but unlike traditional grid-based methods such as Reaction-Diffusion Master Equation (RDME) algorithms, the results of SRBD are statistically independent of the size of the grid used to accelerate the processing of reactions. We use the SRBD algorithm to compute the effective macroscopic reaction rate for both reaction- and diffusion-limited irreversible association in three dimensions. We also study long-time tails in the time correlation functions for reversible association at thermodynamic equilibrium. Finally, we compare different particle and continuum methods on a model exhibiting a Turing-like instability and pattern formation. We find that for models in which particles diffuse off lattice, such as the Doi model, reactions lead to a spurious enhancement of the effective diffusion coefficients.Comment: To appear in J. Chem. Phy

Gene regulatory networks: a coarse-grained, equation-free approach to multiscale computation

We present computer-assisted methods for analyzing stochastic models of gene regulatory networks. The main idea that underlies this equation-free analysis is the design and execution of appropriately-initialized short bursts of stochastic simulations; the results of these are processed to estimate coarse-grained quantities of interest, such as mesoscopic transport coefficients. In particular, using a simple model of a genetic toggle switch, we illustrate the computation of an effective free energy and of a state-dependent effective diffusion coefficient that characterize an unavailable effective Fokker-Planck equation. Additionally we illustrate the linking of equation-free techniques with continuation methods for performing a form of stochastic "bifurcation analysis"; estimation of mean switching times in the case of a bistable switch is also implemented in this equation-free context. The accuracy of our methods is tested by direct comparison with long-time stochastic simulations. This type of equation-free analysis appears to be a promising approach to computing features of the long-time, coarse-grained behavior of certain classes of complex stochastic models of gene regulatory networks, circumventing the need for long Monte Carlo simulations.Comment: 33 pages, submitted to The Journal of Chemical Physic

Diffusive Dynamics of the Reaction Coordinate for Protein Folding Funnels

The quantitative description of model protein folding kinetics using a diffusive collective reaction coordinate is examined. Direct folding kinetics, diffusional coefficients and free energy profiles are determined from Monte Carlo simulations of a 27-mer, 3 letter code lattice model, which corresponds roughly to a small helical protein. Analytic folding calculations, using simple diffusive rate theory, agree extremely well with the full simulation results. Folding in this system is best seen as a diffusive, funnel-like process.Comment: LaTeX 12 pages, figures include

MSM/RD: Coupling Markov state models of molecular kinetics with reaction-diffusion simulations

Molecular dynamics (MD) simulations can model the interactions between macromolecules with high spatiotemporal resolution but at a high computational cost. By combining high-throughput MD with Markov state models (MSMs), it is now possible to obtain long-timescale behavior of small to intermediate biomolecules and complexes. To model the interactions of many molecules at large lengthscales, particle-based reaction-diffusion (RD) simulations are more suitable but lack molecular detail. Thus, coupling MSMs and RD simulations (MSM/RD) would be highly desirable, as they could efficiently produce simulations at large time- and lengthscales, while still conserving the characteristic features of the interactions observed at atomic detail. While such a coupling seems straightforward, fundamental questions are still open: Which definition of MSM states is suitable? Which protocol to merge and split RD particles in an association/dissociation reaction will conserve the correct bimolecular kinetics and thermodynamics? In this paper, we make the first step towards MSM/RD by laying out a general theory of coupling and proposing a first implementation for association/dissociation of a protein with a small ligand (A + B C). Applications on a toy model and CO diffusion into the heme cavity of myoglobin are reported