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

Varying the resolution of the Rouse model on temporal and spatial scales: application to multiscale modelling of DNA dynamics

A multi-resolution bead-spring model for polymer dynamics is developed as a generalization of the Rouse model. A polymer chain is described using beads of variable sizes connected by springs with variable spring constants. A numerical scheme which can use different timesteps to advance the positions of different beads is presented and analyzed. The position of a particular bead is only updated at integer multiples of the timesteps associated with its connecting springs. This approach extends the Rouse model to a multiscale model on both spatial and temporal scales, allowing simulations of localized regions of a polymer chain with high spatial and temporal resolution, while using a coarser modelling approach to describe the rest of the polymer chain. A method for changing the model resolution on-the-fly is developed using the Metropolis-Hastings algorithm. It is shown that this approach maintains key statistics of the end-to-end distance and diffusion of the polymer filament and makes computational savings when applied to a model for the binding of a protein to the DNA filament.Comment: Submitted to Multiscale Modeling and Simulatio

Dynamics of protein-protein encounter: a Langevin equation approach with reaction patches

We study the formation of protein-protein encounter complexes with a Langevin equation approach that considers direct, steric and thermal forces. As three model systems with distinctly different properties we consider the pairs barnase:barstar, cytochrome c:cytochrome c peroxidase and p53:MDM2. In each case, proteins are modeled either as spherical particles, as dipolar spheres or as collection of several small beads with one dipole. Spherical reaction patches are placed on the model proteins according to the known experimental structures of the protein complexes. In the computer simulations, concentration is varied by changing box size. Encounter is defined as overlap of the reaction patches and the corresponding first passage times are recorded together with the number of unsuccessful contacts before encounter. We find that encounter frequency scales linearly with protein concentration, thus proving that our microscopic model results in a well-defined macroscopic encounter rate. The number of unsuccessful contacts before encounter decreases with increasing encounter rate and ranges from 20-9000. For all three models, encounter rates are obtained within one order of magnitude of the experimentally measured association rates. Electrostatic steering enhances association up to 50-fold. If diffusional encounter is dominant (p53:MDM2) or similarly important as electrostatic steering (barnase:barstar), then encounter rate decreases with decreasing patch radius. More detailed modeling of protein shapes decreases encounter rates by 5-95 percent. Our study shows how generic principles of protein-protein association are modulated by molecular features of the systems under consideration. Moreover it allows us to assess different coarse-graining strategies for the future modelling of the dynamics of large protein complexes

Reaction-diffusion kinetics on lattice at the microscopic scale

Lattice-based stochastic simulators are commonly used to study biological reaction-diffusion processes. Some of these schemes that are based on the reaction-diffusion master equation (RDME), can simulate for extended spatial and temporal scales but cannot directly account for the microscopic effects in the cell such as volume exclusion and diffusion-influenced reactions. Nonetheless, schemes based on the high-resolution microscopic lattice method (MLM) can directly simulate these effects by representing each finite-sized molecule explicitly as a random walker on fine lattice voxels. The theory and consistency of MLM in simulating diffusion-influenced reactions have not been clarified in detail. Here, we examine MLM in solving diffusion-influenced reactions in 3D space by employing the Spatiocyte simulation scheme. Applying the random walk theory, we construct the general theoretical framework underlying the method and obtain analytical expressions for the total rebinding probability and the effective reaction rate. By matching Collins-Kimball and lattice-based rate constants, we obtained the exact expressions to determine the reaction acceptance probability and voxel size. We found that the size of voxel should be about 2% larger than the molecule. MLM is validated by numerical simulations, showing good agreement with the off-lattice particle-based method, eGFRD. MLM run time is more than an order of magnitude faster than eGFRD when diffusing macromolecules with typical concentrations in the cell. MLM also showed good agreements with eGFRD and mean-field models in case studies of two basic motifs of intracellular signaling, the protein production-degradation process and the dual phosphorylation cycle. Moreover, when a reaction compartment is populated with volume-excluding obstacles, MLM captures the non-classical reaction kinetics caused by anomalous diffusion of reacting molecules

Reactive Boundary Conditions as Limits of Interaction Potentials for Brownian and Langevin Dynamics

A popular approach to modeling bimolecular reactions between diffusing molecules is through the use of reactive boundary conditions. One common model is the Smoluchowski partial absorption condition, which uses a Robin boundary condition in the separation coordinate between two possible reactants. This boundary condition can be interpreted as an idealization of a reactive interaction potential model, in which a potential barrier must be surmounted before reactions can occur. In this work we show how the reactive boundary condition arises as the limit of an interaction potential encoding a steep barrier within a shrinking region in the particle separation, where molecules react instantly upon reaching the peak of the barrier. The limiting boundary condition is derived by the method of matched asymptotic expansions, and shown to depend critically on the relative rate of increase of the barrier height as the width of the potential is decreased. Limiting boundary conditions for the same interaction potential in both the overdamped Fokker-Planck equation (Brownian Dynamics), and the Kramers equation (Langevin Dynamics) are investigated. It is shown that different scalings are required in the two models to recover reactive boundary conditions that are consistent in the high friction limit where the Kramers equation solution converges to the solution of the Fokker-Planck equation.Comment: 23 pages, 2 figure

Reaction Brownian Dynamics and the effect of spatial fluctuations on the gain of a push-pull network

Brownian Dynamics algorithms are widely used for simulating soft-matter and biochemical systems. In recent times, their application has been extended to the simulation of coarse-grained models of cellular networks in simple organisms. In these models, components move by diffusion, and can react with one another upon contact. However, when reactions are incorporated into a Brownian Dynamics algorithm, attention must be paid to avoid violations of the detailed-balance rule, and therefore introducing systematic errors in the simulation. We present a Brownian Dynamics algorithm for reaction-diffusion systems that rigorously obeys detailed balance for equilibrium reactions. By comparing the simulation results to exact analytical results for a bimolecular reaction, we show that the algorithm correctly reproduces both equilibrium and dynamical quantities. We apply our scheme to a ``push-pull'' network in which two antagonistic enzymes covalently modify a substrate. Our results highlight that the diffusive behaviour of the reacting species can reduce the gain of the response curve of this network.Comment: 25 pages, 7 figures, submitted to Journal of Chemical Physic

Jump-Diffusion Approximation of Stochastic Reaction Dynamics: Error bounds and Algorithms

Biochemical reactions can happen on different time scales and also the abundance of species in these reactions can be very different from each other. Classical approaches, such as deterministic or stochastic approach, fail to account for or to exploit this multi-scale nature, respectively. In this paper, we propose a jump-diffusion approximation for multi-scale Markov jump processes that couples the two modeling approaches. An error bound of the proposed approximation is derived and used to partition the reactions into fast and slow sets, where the fast set is simulated by a stochastic differential equation and the slow set is modeled by a discrete chain. The error bound leads to a very efficient dynamic partitioning algorithm which has been implemented for several multi-scale reaction systems. The gain in computational efficiency is illustrated by a realistically sized model of a signal transduction cascade coupled to a gene expression dynamics.Comment: 32 pages, 7 figure