Stochastic modelling of reaction-diffusion processes:\ud algorithms for bimolecular reactions

Several stochastic simulation algorithms (SSAs) have been recently proposed for modelling reaction-diffusion processes in cellular and molecular biology. In this paper, two commonly used SSAs are studied. The first SSA is an on-lattice model described by the reaction-diffusion master equation. The second SSA is an off-lattice model based on the simulation of Brownian motion of individual molecules and their reactive collisions. In both cases, it is shown that the commonly used implementation of bimolecular reactions (i.e. the reactions of the form A+B → C, or A+A → C) might lead to incorrect results. Improvements of both SSAs are suggested which overcome the difficulties highlighted. In particular, a formula is presented for the smallest possible compartment size (lattice spacing) which can be correctly implemented in the first model. This implementation uses a new formula for the rate of bimolecular reactions per compartment (lattice site)

Time scale of random sequential adsorption

A simple multiscale approach to the diffusion-driven adsorption from a solution to a solid surface is presented. The model combines two important features of the adsorption process: (i) the kinetics of the chemical reaction between adsorbing molecules and the surface; and (ii) geometrical constraints on the surface made by molecules which are already adsorbed. The process (i) is modelled in a diffusion-driven context, i.e. the conditional probability of adsorbing a molecule provided that the molecule hits the surface is related to the macroscopic surface reaction rate. The geometrical constraint (ii) is modelled using random sequential adsorption (RSA), which is the sequential addition of molecules at random positions on a surface; one attempt to attach a molecule is made per one RSA simulation time step. By coupling RSA with the diffusion of molecules in the solution above the surface the RSA simulation time step is related to the real physical time. The method is illustrated on a model of chemisorption of reactive polymers to a virus surface

From Brownian Dynamics to Markov Chain: an Ion Channel Example

A discrete rate theory for general multi-ion channels is presented, in which the continuous dynamics of ion diffusion is reduced to transitions between Markovian discrete states. In an open channel, the ion permeation process involves three types of events: an ion entering the channel, an ion escaping from the channel, or an ion hopping between different energy minima in the channel. The continuous dynamics leads to a hierarchy of Fokker-Planck equations, indexed by channel occupancy. From these the mean escape times and splitting probabilities (denoting from which side an ion has escaped) can be calculated. By equating these with the corresponding expressions from the Markov model the Markovian transition rates can be determined. The theory is illustrated with a two-ion one-well channel. The stationary probability of states is compared with that from both Brownian dynamics simulation and the hierarchical Fokker-Planck equations. The conductivity of the channel is also studied, and the optimal geometry maximizing ion flux is computed.Comment: submitted to SIAM Journal on Applied Mathematic

Realistic boundary conditions for stochastic simulations of reaction-diffusion processes

Many cellular and subcellular biological processes can be described in terms of diffusing and chemically reacting species (e.g. enzymes). Such reaction-diffusion processes can be mathematically modelled using either deterministic partial-differential equations or stochastic simulation algorithms. The latter provide a more detailed and precise picture, and several stochastic simulation algorithms have been proposed in recent years. Such models typically give the same description of the reaction-diffusion processes far from the boundary of the simulated domain, but the behaviour close to a reactive boundary (e.g. a membrane with receptors) is unfortunately model-dependent. In this paper, we study four different approaches to stochastic modelling of reaction-diffusion problems and show the correct choice of the boundary condition for each model. The reactive boundary is treated as partially reflective, which means that some molecules hitting the boundary are adsorbed (e.g. bound to the receptor) and some molecules are reflected. The probability that the molecule is adsorbed rather than reflected depends on the reactivity of the boundary (e.g. on the rate constant of the adsorbing chemical reaction and on the number of available receptors), and on the stochastic model used. This dependence is derived for each model.Comment: 24 pages, submitted to Physical Biolog

The Two Regime method for optimizing stochastic reaction-diffusion simulations

The computer simulation of stochastic reaction-diffusion processes in biology is often done using either compartment-based (spatially discretized) simulations or molecular-based (Brownian dynamics) approaches. Compartment-based approaches can yield quick and accurate mesoscopic results but lack the level of detail that is characteristic of the more computationally intensive molecular-based models. Often microscopic detail is only required in a small region but currently the best way to achieve this detail is to use a resource intensive model over the whole domain. We introduce the Two Regime Method (TRM) in which a molecular-based algorithm is used in part of the computational domain and a compartment-based approach is used elsewhere in the computational domain. We apply the TRM to two test problems including a model from developmental biology. We thereby show that the TRM is accurate and subsequently may be used to inspect both mesoscopic and microscopic detail of reaction diffusion simulations according to the demands of the modeller

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

Stochastic modelling of reaction-diffusion processes: algorithms for bimolecular reactions

Several stochastic simulation algorithms (SSAs) have been recently proposed for modelling reaction-diffusion processes in cellular and molecular biology. In this paper, two commonly used SSAs are studied. The first SSA is an on-lattice model described by the reaction-diffusion master equation. The second SSA is an off-lattice model based on the simulation of Brownian motion of individual molecules and their reactive collisions. In both cases, it is shown that the commonly used implementation of bimolecular reactions (i.e. the reactions of the form A + B -> C, or A + A -> C) might lead to incorrect results. Improvements of both SSAs are suggested which overcome the difficulties highlighted. In particular, a formula is presented for the smallest possible compartment size (lattice spacing) which can be correctly implemented in the first model. This implementation uses a new formula for the rate of bimolecular reactions per compartment (lattice site).Comment: 33 pages, submitted to Physical Biolog

Dynamics of polydisperse irreversible adsorption: a pharmacological example

Many drug delivery systems suffer from undesirable interactions with the host immune system. It has been experimentally established that covalent attachment (irreversible adsorption) of suitable macromolecules to the surface of the drug carrier can reduce such undesirable interactions. A fundamental understanding of the adsorption process is still lacking. In this paper, the classical random irreversible adsorption model is generalized to capture certain essential processes involved in pharmacological applications, allowing for macromolecules of different sizes, partial overlapping of the tails of macromolecules, and the influence of reactions with the solvent on the adsorption process. Working in one dimension, an integro-differential evolution equation for the adsorption process is derived, and the asymptotic behavior of the surface area covered and the number of molecules attached to the surface are studied. Finally, equation-free dynamic renormalization tools are applied to study the asymptotically self-similar behavior of the adsorption statistics

Multiscale reaction-diffusion algorithms: PDE-assisted Brownian dynamics

Two algorithms that combine Brownian dynamics (BD) simulations with mean-field partial differential equations (PDEs) are presented. This PDE-assisted Brownian dynamics (PBD) methodology provides exact particle tracking data in parts of the domain, whilst making use of a mean-field reaction-diffusion PDE description elsewhere. The first PBD algorithm couples BD simulations with PDEs by randomly creating new particles close to the interface which partitions the domain and by reincorporating particles into the continuum PDE-description when they cross the interface. The second PBD algorithm introduces an overlap region, where both descriptions exist in parallel. It is shown that to accurately compute variances using the PBD simulation requires the overlap region. Advantages of both PBD approaches are discussed and illustrative numerical examples are presented.Comment: submitted to SIAM Journal on Applied Mathematic

Analysis of a stochastic chemical system close to a sniper bifurcation of its mean field model

A framework for the analysis of stochastic models of chemical systems for which the deterministic mean-field description is undergoing a saddle-node infinite period (SNIPER) bifurcation is presented. Such a bifurcation occurs for example in the modelling of cell-cycle regulation. It is shown that the stochastic system possesses oscillatory solutions even for parameter values for which the mean-field model does not oscillate. The dependence of the mean period of these oscillations on the parameters of the model (kinetic rate constants) and the size of the system (number of molecules present) is studied. Our approach is based on the chemical Fokker Planck equation. To get some insights into advantages and disadvantages of the method, a simple one-dimensional chemical switch is first analyzed, before the chemical SNIPER problem is studied in detail. First, results obtained by solving the Fokker-Planck equation numerically are presented. Then an asymptotic analysis of the Fokker-Planck equation is used to derive explicit formulae for the period of oscillation as a function of the rate constants and as a function of the system size