Exact on-lattice stochastic reaction-diffusion simulations using partial-propensity methods

Stochastic reaction-diffusion systems frequently exhibit behavior that is not predicted by deterministic simulation models. Stochastic simulation methods, however, are computationally expensive. We present a more efficient stochastic reaction-diffusion simulation algorithm that samples realizations from the exact solution of the reaction-diffusion master equation. The present algorithm, called partial-propensity stochastic reaction-diffusion (PSRD) method, uses an on-lattice discretization of the reaction-diffusion system and relies on partial-propensity methods for computational efficiency. We describe the algorithm in detail, provide a theoretical analysis of its computational cost, and demonstrate its computational performance in benchmarks. We then illustrate the application of PSRD to two- and three-dimensional pattern-forming Gray-Scott systems, highlighting the role of intrinsic noise in these systems.     

Implications of the Turing completeness of reaction-diffusion models, informed by GPGPU simulations on an XBox 360: Cardiac arrhythmias, re-entry and the Halting problem

In the arsenal of tools that a computational modeller can bring to bare on the study of cardiac arrhythmias, the most widely used and arguably the most successful is that of an excitable medium, a special case of a reaction-diffusion model. These are used to simulate the internal chemical reactions of a cardiac cell and the diffusion of their membrane voltages. Via a number of different methodologies it has previously been shown that reaction-diffusion systems are at multiple levels Turing complete. That is, they are capable of computation in the same manner as a universal Turing machine. However, all such computational systems are subject to a limitation known as the Halting problem. By constructing a universal logic gate using a cardiac cell model, we highlight how the Halting problem therefore could limit what it is possible to predict about cardiac tissue, arrhythmias and re-entry.All simulations for this work were carried out on the GPU of an XBox 360 development console, and we also highlight the great gains in computational power and efficiency produced by such general purpose processing on a GPU for cardiac simulations.     

Efficient calculations of coulombic interactions in biomolecular simulations with periodic boundary conditions

To make improved treatments of electrostatic interactions in biomacromolecular simulations, two possibilities are considered. The first is the famous particle–particle and particle–mesh (PPPM) method developed by Hockney and Eastwood, and the second is a new one developed here in their spirit but by the use of the multipole expansion technique suggested by Ladd. It is then numerically found that the new PPPM method gives more accurate results for a two-particle system at small separation of particles. Preliminary numerical examination of the various computational methods for a single configuration of a model BPTI–water system containing about 24,000 particles indicates that both of the PPPM methods give far more accurate values with reasonable computational cost than do the conventional truncation methods. It is concluded the two PPPM methods are nearly comparable in overall performance for the many-particle systems, although the first method has the drawback that the accuracy in the total electrostatic energy is not high for configurations of charged particles randomly generated. © 1993 John Wiley & Sons, Inc.     

Multiscale spatial Monte Carlo simulations: multigriding, computational singular perturbation, and hierarchical stochastic closures

Monte Carlo (MC) simulation of most spatially distributed systems is plagued by several problems, namely, execution of one process at a time, large separation of time scales of various processes, and large length scales. Recently, a coarse-grained Monte Carlo (CGMC) method was introduced that can capture large length scales at reasonable computational times. An inherent assumption in this CGMC method revolves around a mean-field closure invoked in each coarse cell that is inaccurate for short-ranged interactions. Two new approaches are explored to improve upon this closure. The first employs the local quasichemical approximation, which is applicable to first nearest-neighbor interactions. The second, termed multiscale CGMC method, employs singular perturbation ideas on multiple grids to capture the entire cluster probability distribution function via short microscopic MC simulations on small, fine-grid lattices by taking advantage of the time scale separation of multiple processes. Computational strategies for coupling the fast process at small length scales (fine grid) with the slow processes at large length scales (coarse grid) are discussed. Finally, the binomial tau-leap method is combined with the multiscale CGMC method to execute multiple processes over the entire lattice and provide additional computational acceleration. Numerical simulations demonstrate that in the presence of fast diffusion and slow adsorption and desorption processes the two new approaches provide more accurate solutions in comparison to the previously introduced CGMC method.     

Diffusion-influenced excited-state reversible transfer reactions, A* + Bright harpoon over left harpoonC* + D, with two different lifetimes: theories and simulations

We report accurate Brownian simulation results for the kinetics of the pseudo-first-order diffusion-influenced excited-state reversible transfer reaction A(*) + Bright harpoon over left harpoonC(*) + D with two different lifetimes using two different propagation algorithms. The results are used to test approximate solutions for this many-particle problem. Available theories fail when one of the two reactions or (decay) rate constants is large. To remedy this situation, we develop two uniform approximations, which are based on introducing a generalized Smoluchowski term into the relaxation-time approximation. The best of these is the extended unified theory of reversible target reactions, which reduces correctly in all limits and exhibits superior agreement with simulations.     

Sample Duplication Method for Monte Carlo Simulation of Large Reaction-Diffusion System

The sample duplication method for the Monte Carlo simulation of large reaction-diffusion system is proposed in this paper. It is proved that the sample duplication method will effectively raise the efficiency and statistical precision of the simulation without changing the kinetic behaviour of the reaction-diffusion system and the critical condition for the bifurcation of the steady-states. The method has been applied to the simulation of spatial and time dissipative structure of Brusselator under the Dirichlet boundary condition. The results presented in this paper definitely show that the sample duplication method provides a very efficient way to sol-'e the master equation of large reaction-diffusion system. For the case of two-dimensional system, it is found that the computation time is reduced at least by a factor of two orders of magnitude compared to the algorithm reported in literature.     

SCF equations for pure spin states with many-particle interactions

SCF equations for any pure spin state are given for a spin-free system with many-particle interactions. The equations are very simple and explicit. Due to the use of different antisymmetric requirement our equations are different from some of the other methods. In our method, the abstract group theory formalism is converted into some explicit and straightforward equations which makes the many-particle interaction problem easier to handle.     

Comparative analysis of nucleotide translocation through protein nanopores using steered molecular dynamics and an adaptive biasing force

The translocation of nucleotide molecules across biological and synthetic nanopores has attracted attention as a next generation technique for sequencing DNA. Computer simulations have the ability to provide atomistic‐level insight into important states and processes, delivering a means to develop a fundamental understanding of the translocation event, for example, by extracting the free energy of the process. Even with current supercomputing facilities, the simulation of many‐atom systems in fine detail is limited to shorter timescales than the real events they attempt to recreate. This imposes the need for enhanced simulation techniques that expand the scope of investigation in a given timeframe. There are numerous free energy calculation and translocation methodologies available, and it is by no means clear which method is best applied to a particular problem. This article explores the use of two popular free energy calculation methodologies in a nucleotide‐nanopore translocation system, using the α‐hemolysin nanopore. The first uses constant velocity‐steered molecular dynamics (cv‐SMD) in conjunction with Jarzynski's equality. The second applies an adaptive biasing force (ABF), which has not previously been applied to the nucleotide‐nanpore system. The purpose of this study is to provide a comprehensive comparison of these methodologies, allowing for a detailed comparative assessment of the scientific merits, the computational cost, and the statistical quality of the data obtained from each technique. We find that the ABF method produces results that are closer to experimental measurements than those from cv‐SMD, whereas the net errors are smaller for the same computational cost. © 2014 The Authors Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

An accelerated algorithm for discrete stochastic simulation of reaction-diffusion systems using gradient-based diffusion and tau-leaping

Stochastic simulation of reaction-diffusion systems enables the investigation of stochastic events arising from the small numbers and heterogeneous distribution of molecular species in biological cells. Stochastic variations in intracellular microdomains and in diffusional gradients play a significant part in the spatiotemporal activity and behavior of cells. Although an exact stochastic simulation that simulates every individual reaction and diffusion event gives a most accurate trajectory of the system's state over time, it can be too slow for many practical applications. We present an accelerated algorithm for discrete stochastic simulation of reaction-diffusion systems designed to improve the speed of simulation by reducing the number of time-steps required to complete a simulation run. This method is unique in that it employs two strategies that have not been incorporated in existing spatial stochastic simulation algorithms. First, diffusive transfers between neighboring subvolumes are based on concentration gradients. This treatment necessitates sampling of only the net or observed diffusion events from higher to lower concentration gradients rather than sampling all diffusion events regardless of local concentration gradients. Second, we extend the non-negative Poisson tau-leaping method that was originally developed for speeding up nonspatial or homogeneous stochastic simulation algorithms. This method calculates each leap time in a unified step for both reaction and diffusion processes while satisfying the leap condition that the propensities do not change appreciably during the leap and ensuring that leaping does not cause molecular populations to become negative. Numerical results are presented that illustrate the improvement in simulation speed achieved by incorporating these two new strategies.     

Flooding in GROMACS: accelerated barrier crossings in molecular dynamics

The major bottleneck of today's atomistic molecular dynamics (MD) simulations is that because of the enormous computational effort involved, only processes at nanoseconds to microseconds time scales or faster can be studied directly. Unfortunately, apart from a few exceptions, relevant processes, such as chemical reactions or many large scale conformational transitions in proteins, occur at slower time scales and therefore are currently far out of reach for conventional MD. The flooding technique addresses this problem by inclusion of a flooding potential into the force field. This flooding potential locally destabilizes the educt state and thereby significantly accelerates the escape from the initial energy well without affecting the reaction pathway. Here, we summarize the theory and method for the computational chemistry community and detail the implementation within the official version 3.3 of the freely available MD program package GROMACS. Two examples shall demonstrate the application of flooding to accelerate conformational transitions and chemical reactions. The second example was carried out within a QM/MM framework.     

The distributed diagonal force decomposition method for parallelizing molecular dynamics simulations

Parallelization is an effective way to reduce the computational time needed for molecular dynamics simulations. We describe a new parallelization method, the distributed-diagonal force decomposition method, with which we extend and improve the existing force decomposition methods. Our new method requires less data communication during molecular dynamics simulations than replicated data and current force decomposition methods, increasing the parallel efficiency. It also dynamically load-balances the processors' computational load throughout the simulation. The method is readily implemented in existing molecular dynamics codes and it has been incorporated into the CHARMM program, allowing its immediate use in conjunction with the many molecular dynamics simulation techniques that are already present in the program. We also present the design of the Force Decomposition Machine, a cluster of personal computers and networks that is tailored to running molecular dynamics simulations using the distributed diagonal force decomposition method. The design is expandable and provides various degrees of fault resilience. This approach is easily adaptable to computers with Graphics Processing Units because it is independent of the processor type being used.     

Toward realistic modeling of dynamic processes in cell signaling: quantification of macromolecular crowding effects

One of the major factors distinguishing molecular processes in vivo from biochemical experiments in vitro is the effect of the environment produced by macromolecular crowding in the cell. To achieve a realistic modeling of processes in the living cell based on biochemical data, it becomes necessary, therefore, to consider such effects. We describe a protocol based on Brownian dynamics simulation to characterize and quantify the effect of various forms of crowding on diffusion and bimolecular association in a simple model of interacting hard spheres. We show that by combining the elastic collision method for hard spheres and the mean field approach for hydrodynamic interaction (HI), our simulations capture the correct dynamics of a monodisperse system. The contributions from excluded volume effect and HI to the crowding effect are thus quantified. The dependence of the results on size distribution of each component in the system is illustrated, and the approach is applied as well to the crowding effect on electrostatic-driven association in both neutral and charged environments; values for effective diffusion constants and association rates are obtained for the specific conditions. The results from our simulation approach can be used to improve the modeling of cell signaling processes without additional computational burdens.     

Molecular dynamics of cryptophane and its complexes with tetramethylammonium and neopentane using a continuum solvent model

Time scales currently obtainable in explicit–solvent molecular dynamics simulations are inadequate for the study of many biologically important processes. This has led to increased interest in the use of continuum solvent models. For such models to be used effectively, it is important that their behavior relative to explicit simulation be clearly understood. Accordingly, 5 ns stochastic dynamics simulations of a derivative of cryptophane-E alone, and complexed with tetramethylammonium and neopentane were carried out. Solvation electrostatics were accounted for via solutions to the Poisson equation. Nonelectrostatic aspects of solvation were incorporated using a surface area-dependent energy term. Comparison of the trajectories to those from previously reported 25 ns explicit–solvent simulations shows that use of a continuum solvent model results in enhanced sampling. Use of the continuum solvent model also results in a considerable increase in computational efficiency. The continuum solvent model is found to predict qualitative structural characteristics that are similar to those observed in explicit solvent. However, some differences are significant, and optimization of the continuum parameterization will be required for this method to become an efficient alternative to explicit–solvent simulation. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 956–970, 1999     

"All possible steps" approach to the accelerated use of Gillespie's algorithm

Many physical and biological processes are stochastic in nature. Computational models and simulations of such processes are a mathematical and computational challenge. The basic stochastic simulation algorithm was published by Gillespie about three decades ago [J. Phys. Chem. 81, 2340 (1977)]. Since then, intensive work has been done to make the algorithm more efficient in terms of running time. All accelerated versions of the algorithm are aimed at minimizing the running time required to produce a stochastic trajectory in state space. In these simulations, a necessary condition for reliable statistics is averaging over a large number of simulations. In this study the author presents a new accelerating approach which does not alter the stochastic algorithm, but reduces the number of required runs. By analysis of collected data the author demonstrates high precision levels with fewer simulations. Moreover, the suggested approach provides a good estimation of statistical error, which may serve as a tool for determining the number of required runs.     

Multiscale QM/MM molecular dynamics study on the first steps of guanine damage by free hydroxyl radicals in solution

Understanding the damage of DNA bases from hydrogen abstraction by free OH radicals is of particular importance to understanding the indirect effect of ionizing radiation. Previous studies address the problem with truncated DNA bases as ab initio quantum simulations required to study such electronic-spin-dependent processes are computationally expensive. Here, for the first time, we employ a multiscale and hybrid quantum mechanical-molecular mechanical simulation to study the interaction of OH radicals with a guanine-deoxyribose-phosphate DNA molecular unit in the presence of water, where all of the water molecules and the deoxyribose-phosphate fragment are treated with the simplistic classical molecular mechanical scheme. Our result illustrates that the presence of water strongly alters the hydrogen-abstraction reaction as the hydrogen bonding of OH radicals with water restricts the relative orientation of the OH radicals with respect to the DNA base (here, guanine). This results in an angular anisotropy in the chemical pathway and a lower efficiency in the hydrogen-abstraction mechanisms than previously anticipated for identical systems in vacuum. The method can easily be extended to single- and double-stranded DNA without any appreciable computational cost as these molecular units can be treated in the classical subsystem, as has been demonstrated here.     

On the structural convergence of biomolecular simulations by determination of the effective sample size

Although atomistic simulations of proteins and other biological systems are approaching microsecond timescales, the quality of simulation trajectories has remained difficult to assess. Such assessment is critical not only for establishing the relevance of any individual simulation but also in the extremely active field of developing computational methods. Here we map the trajectory assessment problem onto a simple statistical calculation of the "effective sample size", that is, the number of statistically independent configurations. The mapping is achieved by asking the question, "How much time must elapse between snapshots included in a sample for that sample to exhibit the statistical properties expected for independent and identically distributed configurations?" Our method is more general than autocorrelation methods in that it directly probes the configuration-space distribution without requiring a priori definition of configurational substates and without any fitting parameters. We show that the method is equally and directly applicable to toy models, peptides, and a 72-residue protein model. Variants of our approach can readily be applied to a wide range of physical and chemical systems.     

Efficient stochastic simulation of simultaneous reaction and diffusion in a gas-liquid interface

A stochastic simulation of simultaneous reaction and diffusion is proposed for the gas-liquid interface formed in the surface of a gas bubble within a liquid. The interface between a carbon dioxide bubble and an aqueous solution of calcium hydroxide was simulated as an application example, taken from the integrated production of calcium carbonate. First Gillespie’s stochastic simulation algorithm was applied in separate reaction and diffusion simulations. The results from these simulations were consistent with deterministic solutions based on differential equations. However it was observed that stochastic diffusion simulations are extremely slow. The sampling of diffusion events was accelerated applying a group molecule transfer scheme based on the binomial distribution function. Simulations of the reaction-diffusion in the gas-liquid interface based on the standard Gillespie’s stochastic algorithm were also slow. However the application of the binomial distribution function scheme allowed to compute the concentration profiles in the gas-liquid interface in a fraction of the time required with the standard Gillespie’s stochastic algorithm.