Source-linked AI summary
Simulation of stochastic reaction-diffusion processes on unstructured meshes
Stefan Engblom, Lars Ferm, Andreas Hellander, Per Lötstedt
TL;DR
The paper addresses stochastic reaction-diffusion simulation when molecular discreteness and spatial variation matter. It combines an RDME with FEM-derived diffusion on unstructured meshes and evaluates a hybrid mesoscopic–macroscopic algorithm, which is reported to be up to three orders of magnitude faster than SSA in one example.
Problem
Small molecule copy numbers and spatial inhomogeneity require stochastic, spatially resolved biochemical models.
Method
The paper simulates the RDME on unstructured meshes, deriving mesoscopic diffusion coefficients from FEM discretization of the macroscopic diffusion operator.
Results
The hybrid algorithm is up to three orders of magnitude faster than SSA in the four-species comparison.
Takeaways & Limitations
The method provides a flexible choice between mesoscopic and macroscopic diffusion treatment while retaining reported convergence and sampling accuracy.
Takeaways & Limitations
The reported examples are limited to two-dimensional simulations and currently use a uniform diffusion constant.
Abstract
from arXiv · showhide
Stochastic chemical systems with diffusion are modeled with a reaction-diffusion master equation. On a macroscopic level, the governing equation is a reaction-diffusion equation for the averages of the chemical species. On a mesoscopic level, the master equation for a well stirred chemical system is combined with Brownian motion in space to obtain the reaction-diffusion master equation. The space is covered by an unstructured mesh and the diffusion coefficients on the mesoscale are obtained from a finite element discretization of the Laplace operator on the macroscale. The resulting method is a flexible hybrid algorithm in that the diffusion can be handled either on the meso- or on the macroscale level. The accuracy and the efficiency of the method are illustrated in three numerical examples inspired by molecular biology.
Abbreviations
The paper lists abbreviations for its stochastic, reaction-diffusion, numerical, and continuum-modeling terminology.
- The abbreviation table presents expanded forms for terms including CME, RDME, SSA, RRE, FEM, PDE, and RDE.
1 Introduction
Biochemical systems require stochastic and spatially resolved descriptions when molecule numbers are small or species are spatially inhomogeneous. The paper develops an RDME method on unstructured meshes, combining finite-element diffusion coefficients with hybrid mesoscopic–macroscopic simulation.
- Small molecule copy numbers and reaction probabilities make intrinsic noise important for biochemical-cell behavior.
- Reaction rate equations provide deterministic approximations, but their validity is limited when molecule numbers are small.
- Spatial inhomogeneity motivates the reaction-diffusion master equation, which tracks stochastic species distributions across space.
- The RDME has KN-dimensional state space for N species and K cells, making direct solution impractical and motivating Monte Carlo sampling.
- The method uses unstructured meshes, FEM-derived mesoscopic diffusion coefficients, and Strang-split hybrid integration for geometric flexibility and efficiency.
- The paper reports convergence to macroscopic diffusion and accurate RDME sampling at a fraction of the time required for full simulation.
2 Reaction-Diffusion Master Equation
The RDME represents reactions within mesh cells and diffusion between neighboring cells, while its macroscopic limit is a reaction-diffusion equation. Diffusion coefficients are linked to cell geometry, Brownian motion, and FEM or finite-difference discretizations.
- The computational domain is partitioned into nonoverlapping cells containing copy numbers for each chemical species.
- Chemical reactions change species counts within cells, whereas diffusion moves molecules between adjacent cells; the RDME evolves their joint probability density.
- Diffusion propensities connect neighboring cells, with coefficients depending on the macroscopic diffusion constant, local length scale, and cell shapes.
- The diffusion coefficient q_kj is interpreted through the expected first-exit time of a molecule from cell C_k to cell C_j.
- For Cartesian meshes, q = γ/h^2 yields convergence of the discrete diffusion solution to the corresponding PDE as h approaches zero.
- The macroscopic reaction-diffusion equation combines reaction terms with γ∆φ_i, while the discrete diffusion matrix follows mesh connectivity.
- FEM-derived coefficients recover uniform-grid jump rates and agree with Brownian first-exit probabilities, although the exit time is not generally exponentially distributed.
3 Diffusion coefficients
The paper constructs mesoscopic diffusion coefficients on an unstructured dual mesh from a finite element discretization, preserving the connectivity and relevant sign properties of the macroscopic operator. The resulting approximation is second-order accurate in 2D under the stated conditions.
- Mesh construction: The dual mesh uses polygonal cells centered at primal-mesh vertices, with boundaries formed by triangle bisectors, edge midpoints, and domain boundaries.In 2D, primal cells are triangles and dual cells are polygons; analogous primal and dual constructions are described in 1D.
- Finite element discretization: Mass lumping assigns each diagonal entry A_jj the length, area, or volume of the corresponding dual cell in 1D, 2D, or 3D.The lumped system uses these diagonal entries to represent dual-cell measures.
- Finite element discretization: The finite element system yields the approximation φ̇ = ω̂(φ) + γD̂φ, with D̂ = Â^-1Ŝ, and is generally second-order accurate in 2D.The diffusion matrix inherits sparse coupling from the stiffness matrix.
- Mesh conditions: Delaunay triangulations in 2D and meshes with non-obtuse tetrahedral dihedral angles in 3D ensure correctly signed off-diagonal diffusion entries and support a discrete maximum principle.With Dirichlet conditions and nonnegative initial data, the solution remains nonnegative under these mesh conditions.
- Mesoscopic coefficients: Mesoscopic jump coefficients are obtained from the macroscopic diffusion matrix, so molecules jump only between dual cells sharing a point, edge, or facet.Positive off-diagonal entries are interpreted through the connectivity graph and inverse expected first-exit times.
- Alternative discretization: The finite volume alternative differs only in its diffusion-term approximation, and its approximation quality can depend critically on mesh quality.The paper therefore prefers finite elements for the predominantly diffusive molecular-biology systems considered.
4 Moments of the diffusion
The paper derives exact moment equations for diffusion without reactions and uses them to characterize deterministic approximation error. The variance scales with the expected molecule number, making expectation-based diffusion reliable at high copy numbers but unsuitable for cells with few molecules.
- First and second moments: Because diffusion propensities are linear, the first-moment equation is exact and does not couple to higher-order moments.The covariance equation is also exact for the same linear-propensity structure.
- Mean behavior: The mean perturbation decays at most as cδ exp(λ2t), where λ2 is the negative eigenvalue of Q with smallest magnitude.This bounds the relaxation of the expected molecule distribution toward its long-time behavior.
- Covariance behavior: If the initial covariance is zero, the covariance norm is proportional to κ_i for large time.The covariance remains bounded as t approaches infinity.
- Stochastic fluctuations: The variance is proportional to the expected molecule number over bounded time intervals, so the standard deviation grows like the square root of the expectation.At large time, the expectation is of the same order as κ_i.
- Deterministic approximation: For large copy numbers, the standard-deviation-to-mean ratio is small, whereas expectation values should not describe cells containing few molecules.This result motivates deterministic macroscopic diffusion for sufficiently populated species or cells.
5 Time integration and hybrid diffusion
The hybrid algorithm splits diffusion by copy number, treating low-copy species stochastically and high-copy species macroscopically within a Strang-split time integrator. The method provides accuracy and preservation results, while its stochastic-to-macroscopic approximation requires assumptions about diffusion noise.
- Hybrid diffusion: Diffusion is split into low-copy species handled by SSA and high-copy species advanced macroscopically, avoiding resolution-dependent stochastic diffusion cost.The SSA diffusion cost is proportional to xγ/h2, where x is the total number of diffusing molecules.
- Time integration: The macroscopic diffusion is applied for half a step, the reduced master equation is integrated for a full step, and macroscopic diffusion is applied for another half step.This is the Strang splitting procedure used to evolve trajectories and approximate the probability density.
- Time integration: The high-copy diffusion step uses the trapezoidal method for second-order temporal accuracy, while backward Euler offers first-order accuracy for stiff integration.Explicit integration is inefficient at the reaction time step, so an implicit stiff solver is required.
- Accuracy: The splitting and trapezoidal time-discretization error is O(∆t2), while the finite-element diffusion approximation error is O(h2).The stated bounds cover temporal splitting, time discretization, and FEM spatial approximation.
- Preservation properties: For diffusion-only systems, trapezoidal integration preserves non-negativity under ∆t ≤ hmin2/6γ, whereas backward Euler has no time-step restriction for non-negativity.The proposition assumes nonnegative initial trajectories and defines hmin as the minimal vertex-to-opposing-edge distance in a mesh triangle.
- Preservation and limitations: The hybrid method preserves each species’ total molecule count under diffusion, but macroscopic diffusion requires assuming that its noise has little or no system-wide impact.Averages can reduce this error because macroscopic diffusion is exact in expectation; short steps may also waste work when many SSA realizations have no events.
6 Numerical results
The numerical examples verify convergence of the unstructured-mesh diffusion model, reproduce domain-separation behavior in a bistable system, and show substantial computational savings from the hybrid method.
- 6.1 Diffusion: The diffusion experiment compares stochastic concentrations with analytical and FEM solutions on meshes of 33 and 123 vertices using 100M molecules across M trajectories.The relative differences δa and δd are evaluated in weighted ℓ2 and ℓ∞ norms.
- 6.1 Diffusion: The mesoscopic diffusion solution converges toward the macroscopic solution, with spatial error behaving as O(h^2) and stochastic differences decaying at the correct rate as molecule number increases.For hmax = 0.5 and 0.25, the ℓ2-errors are 8.7 · 10^-4 and 2.8 · 10^-4, respectively.
- 6.2 Domain separation in a bi-stable system: With γ = 2 × 10^-13m2/s, the bistable system forms spatially separated A- and B-rich domains, including patches in different phases by t = 8.32s.A molecules concentrate in the lower-right region and B molecules in the upper-left region.
- 6.2 Domain separation in a bi-stable system: With γ = 1 × 10^-12m2/s, the same system does not separate into distinct phases and ends in an A-dominated state resembling the homogeneous case.The comparison indicates different spatial behavior for the two diffusion coefficients.
- 6.3 The hybrid method – metabolites and enzymes: The hybrid example treats enzyme reactions and diffusion stochastically while solving metabolite diffusion with a deterministic PDE, and it compares hybrid trajectories with SSA solutions.The comparison uses M = 10^4 realizations and evaluates differences at t = 200 for multiple time steps.
- 6.3 The hybrid method – metabolites and enzymes: The hybrid method is up to 1000 times faster than SSA when diffusion dominates for larger γ, while retaining small differences between the solutions.At t = 200, the speedup with Δt = 5 is about 35 as enzyme concentrations increase.
7 Conclusions
The paper combines stochastic reaction modeling with finite-element-based diffusion on unstructured meshes and uses a hybrid mesoscopic–macroscopic treatment to reduce computation for abundant, rapidly diffusing species.
- Conclusions: The RDME is discretized on an unstructured mesh, with diffusion coefficients derived from a finite-element discretization of the Laplacian.The mesh provides greater geometric flexibility than a Cartesian mesh.
- Conclusions: The hybrid method deterministically integrates diffusion for species with large copy numbers, substantially reducing computation when the diffusion constant is large.This complements stochastic simulation of the mesoscopic system.
- Conclusions: Across three systems, the method shows convergence to the macroscopic solution, agreement with structured-mesh results, and a speedup of up to three orders of magnitude over SSA.The reported examples include a diffusion-only system, a bistable biochemical system, and a four-species hybrid simulation.
- Conclusions: The numerical examples are restricted to two spatial dimensions, although the authors state that extension to three dimensions is straightforward.The implementation currently uses a uniform diffusion constant, with extensions to spatially varying diffusion and convection described as uncomplicated.