PaperlayerView exact revision on arXiv

Source-linked AI summary

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

Radek Erban, S Jonathan Chapman

arXiv:0903.1298v1physics.bio-phq-bio.QM

TL;DR

The paper examines whether bimolecular reactions are correctly represented in two commonly used stochastic reaction-diffusion algorithms. It proposes improvements for both on-lattice and off-lattice models, including a critical compartment-size formula and revised reaction-rate treatment.

  • Problem

    Common implementations of bimolecular reactions in stochastic reaction-diffusion algorithms may produce incorrect results.

  • Method

    The paper studies an on-lattice reaction-diffusion master-equation model and an off-lattice model simulating Brownian motion and reactive molecular collisions.

  • Results

    The study proposes improvements to both algorithms, including a lower bound h ≥ hcrit for compartment size and revised local bimolecular reaction rates.

  • Takeaways & Limitations

    The resulting formulae can be applied to modelling bimolecular reactions within complex reaction-diffusion processes.

  • Takeaways & Limitations

    The off-lattice model may require a simulation time step significantly below a nanosecond for the illustrative reaction parameters, with stricter limits for faster-diffusing molecules.

Abstract

from arXiv · show

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).

1. Introduction

Stochastic simulation algorithms model molecular numbers and positions when deterministic concentration-based PDE models are unsuitable. This paper examines bimolecular reactions in compartment-based and Brownian-motion-based SSAs and proposes improvements.

  • Motivation: Few molecules can make concentration-based PDE modelling unsuitable because the required neighbourhood cannot contain enough molecules.The appropriate variables are then molecular numbers and positions rather than concentrations.
  • Motivation: SSAs provide a more detailed and precise molecular description than deterministic PDE models, while agreeing for simple zero-order and first-order reactions.Bimolecular reactions and reactive boundaries are more delicate cases.
  • Scope and methods: The paper studies bimolecular reactions A+B →C and A+A →C in two commonly used reaction-diffusion SSAs.The first divides space into well-mixed compartments with diffusive jumps; the second simulates individual-molecule Brownian motion and reactive collisions.
  • Challenges: The compartment-based SSA requires a suitable compartment size balancing spatial resolution against the assumptions underlying bimolecular reactions.The molecular-based SSA can require an unrealistically small reaction radius for typical rate constants and diffusion coefficients.
  • Contribution: The paper proposes improved algorithms and derives a reaction-rate relationship connecting simulation parameters with experimentally measurable reaction constants.The improvements target the difficulties identified for both SSA formulations.

2. Bimolecular reactions - two model problems

The paper introduces heteroreaction and homoreaction examples to establish stochastic reference behaviour before evaluating reaction-diffusion models. These examples distinguish reactions between different species from reactions between identical molecules and compare stochastic and deterministic descriptions.

  • Definitions: Bimolecular reactions involve two reacting molecules, with heteroreactions using different species and homoreactions using the same species.Examples include catalytic degradation or heterodimerization for heteroreactions and homodimerization for homoreactions.
  • 2.1. A heteroreaction example: In the heteroreaction example, A is produced and catalytically degraded by B, whose molecule count remains constant.The number of A molecules therefore fluctuates around an equilibrium value after an initial transient.
  • Stochastic formulation: The heteroreaction propensity is proportional to the available A-B pairs, while homoreaction propensities differ and require separate treatment.Reaction occurrence probabilities over an infinitesimal interval are given by the corresponding propensity functions.
  • 2.1. A heteroreaction example: Figure 1 compares one Gillespie SSA trajectory and its average with the stationary distribution estimated by simulation and calculated from formula (3).For k1/ν = 0.2 sec−1, k2ν = 1 sec−1, A(0) = 5, and one B molecule, the mean A value is 5.
  • 2.1. A heteroreaction example: The heteroreaction stationary distribution is uniquely determined by the rate constants and the constant number of B molecules, and is the Poisson distribution.The simulation histogram confirms the distribution given by formula (3).
  • 2.2. A homoreaction example: In the homoreaction example, A undergoes dimerization while another reaction produces A, so dimerization removes two A molecules and production adds one.The stochastic stationary distribution involves a modified Bessel function, whereas the deterministic ODE stationary value is k2/2k1.
  • 2.2. A homoreaction example: For k1/ν = 0.2 sec−1 and k2ν = 10 sec−1, the deterministic stationary value is As = 5 while the stochastic mean is Ms ≈ 5.13.The difference is 2.5%, so the ODE provides a reasonable approximation in this example, although exact equality does not hold.
  • Compartment discretization: The computational domain is divided into K3 cubic compartments of volume h3 = (L/K)3 for the compartment-based model.Diffusion is represented by jumps between neighbouring compartments, and bimolecular reactions occur only within the same compartment.

3. Disadvantages of current SSAs for reaction-diffusion modelling

The compartment-based and molecular-based reaction-diffusion SSAs impose modelling assumptions that can distort bimolecular reactions. The compartment model is sensitive to compartment size, while the molecular model can require unrealistically small reaction radii and time steps.

  • Compartment-based model: The compartment-based SSA divides the domain into well-mixed compartments, permits reactions only within compartments, and models diffusion as jumps between neighbours.Its formulation uses K^3 compartments of side length h = L/K, with Aijk and Bijk molecule counts and diffusion propensities proportional to D/h^2.
  • Compartment-based model: Increasing K, or decreasing h, shifts the heteroreaction stationary distribution φK(n) to the right because bimolecular removal of A becomes slower.The production rate remains unchanged as the compartment size decreases, while the bimolecular reaction is lost in the theoretical limit h →0.
  • Compartment-based model: The same rightward shift occurs for the homoreaction example when the homodimerization reaction is implemented separately in each compartment.The compartment model replaces A + A →B with K^3 local reactions Aijk + Aijk →Bijk.
  • Compartment-based model: Choosing h requires balancing spatial resolution against a lower bound: h must be small enough to resolve spatial variation but not so small that bimolecular reactions are misrepresented.The paper identifies the appropriate compartment size as a fundamental problem and proposes a critical value hcrit for correct implementation of the propensity function.
  • Molecular-based models: In molecular-based models, molecules react whenever their separation is below a reaction radius, but typical parameters can imply a radius smaller than the molecular size.For water at room temperature, the molecular radius estimate is m = 219.7 nm for D = 1 µm2 sec−1; typical protein interaction rates can require unrealistically small reaction radii.
  • Molecular-based models: The molecular-based implementation can also require a simulation time step significantly less than a nanosecond, and treating every collision as reactive neglects non-reactive collisions.The paper proposes modified algorithms that allow the reaction radius to reach the molecular radius and account for non-reactive collisions.

4. Improved SSAs for reaction-diffusion modelling

The paper modifies both compartment-based and molecular-based reaction-diffusion SSAs to address incorrect bimolecular-reaction results. The improved compartment model requires a critical lattice spacing and revised propensities, while the improved molecular model relates reaction radius and event probability to physically meaningful parameters.

  • Improved SSAs for reaction-diffusion modelling: Modified SSAs are proposed to overcome the bimolecular-reaction problems identified in the original algorithms.The modifications cover both compartment-based and molecular-based reaction-diffusion models.
  • Improved compartment-based model: hcrit is the smallest compartment size that can be correctly simulated with the modified heteroreaction propensity.Correct stationary distributions are obtained for h ≥ hcrit, whereas this is not possible for smaller h.
  • Improved compartment-based model: β∞≈0.25272 is the proposed default dimensionless parameter when system-specific information is unavailable.For large K = L/h, β is expected to converge to β∞.
  • Improved compartment-based model: The modified heteroreaction propensity reproduces the correct stationary distribution φ1(n) for K = 16, unlike the original formula.The modified result matches the correct distribution, while the original model shows a difference.
  • Improved compartment-based model: Boundary effects make the true hcrit slightly larger than the approximation in (25).The paper compensates by using β values slightly larger than β∞, with β depending weakly on K.
  • Improved molecular-based models: Reducing Pλ at fixed γ reduces κ and increases the reaction radius, enabling a more physically meaningful molecular collision model.For the heteroreaction example, Pλ = 0.77% and a 40 nm reaction radius produce excellent agreement with the exact stationary distribution.
  • Improved molecular-based models: The λ−̺ model also gives excellent agreement with the exact homoreaction stationary distribution.The homoreaction implementation adjusts the reaction-rate calculation for the removal of two A molecules and uses 2D_A in the relevant formulae.

5. Discussion

The paper links on-lattice and off-lattice SSAs through constraints on spatial resolution and reaction rates, and proposes improved formulations for bimolecular reactions. These results are illustrated computationally, with formulae intended for more complex reaction-diffusion systems.

  • Both on-lattice and off-lattice SSAs were studied, revealing connections between the two model types.
  • h ≥ hcrit is required for the on-lattice compartment size, and the bimolecular reaction rate per compartment becomes infinite at h = hcrit.
  • For sufficiently large ρ, the off-lattice reaction rate λ is given by the bimolecular rate constant divided by the reaction-ball volume.
  • Both models converge for large ρ and large h to local reaction rates equal to the rate constant divided by the reaction volume.
  • The results were demonstrated with illustrative examples, while the general formulae can be applied to bimolecular reactions in complex reaction-diffusion processes.

Glossary

The glossary introduces stochastic simulation terminology and the modified Bessel function used in the paper’s mathematical background.

  • The Gillespie SSA simulates the time evolution of well-stirred chemical systems, consistently with the chemical master equation.
  • A Markov chain is a stochastic process whose future depends only on the present state, not the past.
  • The modified Bessel function of the first kind is a solution of an ordinary differential equation and is available in standard mathematical software.

Appendix A. Stationary distributions, means and variances for the illustrative heteroreaction and homoreaction examples

Appendix A derives stationary distributions, stochastic means, and variances for illustrative heteroreaction and homoreaction systems, alongside deterministic approximations.

  • For the well-stirred heteroreaction example, the probabilities p_n(t) evolve according to the chemical master equation.
  • The stationary distribution φ(n) is uniquely determined by the stationary master equation and normalization condition.
  • For the heteroreaction example, the stationary stochastic variance equals the stationary stochastic mean: V_s = M_s.
  • For the homoreaction example, the stationary distribution and its mean and variance are expressed using modified Bessel functions.
  • The deterministic homoreaction description uses the concentration ODE da/dt = −2k_1a^2 + k_2.
  • A(t) = a(t)ν provides the deterministic approximation of the average molecule number in volume ν, but equation (A.5) does not describe the stochastic mean’s evolution.

Appendix B. Reaction-diffusion master equation

Appendix B formulates the reaction-diffusion master equation using molecule-count configurations, state-transition operators, and terms for reactions and diffusion jumps.

  • Molecule-count configurations n and m assign nonnegative integer populations to indexed lattice compartments.
  • The joint probability p(n, m, t) describes the compartment-wise counts of species A and B over time.
  • Operators J^e_ijk encode transferring one molecule between neighboring compartments.
  • The reaction-diffusion master equation combines chemical reaction terms with diffusion-jump terms for the lattice system.

Appendix C. Derivation of formulae (26) and (27)

The appendix derives formulas (26) and (27) by analyzing stationary reaction–diffusion behavior on a lattice and choosing a compartment reaction rate λ to match the desired mean molecule number. The derivation also accounts for the diffusive motion of both reactants through the relative diffusion constant DA + DB.

  • Stationary lattice model: The derivation first fixes DB = 0 and B0 = 1, leaving one immobile molecule of B and excluding reactions (11).The molecule of B remains in a specified compartment while the reaction–diffusion master equation is simplified accordingly.
  • Stationary lattice model: The reaction–diffusion master equation is transformed into a tensor basis satisfying an orthonormality condition.This basis change reduces the equations to relationships among the coefficients c_ijk and the reaction term at the B molecule’s compartment.
  • Derivation of formula (26): λ is selected so that the stationary average number of A molecules equals the target value Ms from equation (4).The average number is obtained by summing the stationary distribution and evaluating the contribution from the compartment containing B.
  • Derivation of formula (26): When B diffuses, the distance between A and B molecules evolves equivalently to A diffusion with diffusion constant DA + DB.The resulting calculation shows that the stationary mean of A depends on B’s compartment position.
  • Derivation of formula (27): Using the stationary mean relation, equation (26), and orthonormality, the appendix derives formula (27).The derivation obtains equation (27) by substituting the coefficient expressions into the transformed stationary equations.

Appendix D. Derivation of formula (28)

The appendix derives formula (28) by Fourier-transforming the lattice equations for the stationary mean concentration and taking the continuum limit. It also relates the lattice calculation to a relative-distance model for reacting molecules.

  • Fourier and continuum derivation: Formula (27) is identified as a Riemann sum for a definite integral.This provides the passage from the finite lattice expression toward its infinite-lattice or continuum representation.
  • Relative-distance interpretation: A direct derivation of formula (28) fixes B at the origin and lets A diffuse with diffusion constant DA + DB, tracking relative rather than absolute distance.The construction preserves the distribution of relative positions needed for bimolecular reactions.
  • Fourier and continuum derivation: The stationary lattice equations are transformed using Fourier modes, and the lattice operator simplifies to a cosine expression.The transformed quantities include the Fourier transform bµxyz and exponential factors indexed by lattice coordinates.
  • Reaction-rate matching: The reaction term λM000 represents removal of A molecules at the origin and is matched to the deterministic reaction rate k1/h3 M∞.This comparison identifies the lattice reaction-rate relationship used in the derivation.
  • Reaction-rate matching: Taking K →∞ in formula (27), using cosine periodicity, and integrating over z yields formulas (D.1) and (28).The appendix presents both the limiting Riemann-sum route and the resulting integral formula.

Appendix E. Derivation of formula (31)

The appendix derives formula (31) by solving the steady diffusion problem for A around a reactive ball. Boundary conditions determine the concentration profile and the total flux, which is matched to the bimolecular reaction rate.

  • Diffusion-to-a-ball model: The derivation models molecules of A diffusing toward a ball centered at the origin and removes them at rate λ.The equilibrium concentration c(r) is assumed to be continuous with a continuous derivative and satisfies radial diffusion equations.
  • Boundary conditions: The general solutions of the radial second-order ODEs are constrained by the boundary condition at infinity and continuity at the origin.These conditions determine a1 = c∞ and a4 = −a3 before the remaining constants are fixed.
  • Boundary conditions: Continuity of concentration and its derivative at the ball boundary determines the remaining solution constants.The resulting radial profile is then used to compute the flux through the boundary.
  • Reaction-rate identification: The total boundary flux is equated with k1c∞, and division by c∞ yields formula (31).The derivation uses diffusion to the ball and notes that the approximation can be justified using a more general many-particle evolution equation.

Appendix F. Derivation of (34)–(36) and a numerical method for solving it

The appendix formulates a fixed-point equation for the steady radial distribution function of reacting particles and derives a numerical procedure for computing the reaction rate constant κ. The equation is discretized into a linear system and evaluated across γ and Pλ values.

  • Fixed-point formulation: The concentration of A after diffusion over Δt is expressed through an integral kernel K(r, r′; γ).Particle removal within a radius around the origin modifies the diffusion update by the removal probability Pλ.
  • Fixed-point formulation: Equation (34) is the fixed-point equation for the iterative scheme and defines g(r) as a steady-state radial distribution function generalized to arbitrary Pλ.The appendix notes that the earlier RDF calculation covered only Pλ = 1.
  • Numerical results: Substituting the computed g(ri) values into equation (36) gives κ, and repeating for different γ and Pλ produces the results in Figure 5.The numerical sweep varies both the diffusion-related parameter γ and the removal probability Pλ.
  • Numerical solution: The numerical method imposes g(r) →1 as r →∞ and approximates g(r) = 1 beyond a sufficiently large radius S.A two-part mesh resolves the interval from the reaction radius to S.
  • Numerical solution: Discretizing equation (34) produces a linear system for g(ri), solvable for example by Gaussian elimination.The right-hand side can be evaluated using the error function erf.
Loading 0903.1298v1…