Coupling particle-based reaction-diffusion simulations with reservoirs mediated by reaction-diffusion PDEs

TL;DR

This paper tackles the challenge of simulating open biochemical systems by coupling particle-based reaction-diffusion (PBRD) with reservoirs described by reaction-diffusion PDEs. It develops two key mean-field results: (i) a grand-canonical diffusion mean-field that yields a consistent injection rate $\gamma = D/\delta x^2$ for coupling to constant-concentration reservoirs, and (ii) a mean-field mapping between microscopic Doi parameters $(\alpha, \sigma)$ and macroscopic RD rate $\kappa$ (e.g., $\kappa = \alpha\pi\sigma^2$ in 2D). Using these results, the authors implement a hybrid scheme that splits the domain into a particle region and a PDE reservoir region, with Strang splitting to integrate injection, reaction, and diffusion and a robust boundary treatment for bimolecular cross-boundary reactions. Numerical experiments across diffusion, first-order proliferation, and Lotka–Volterra dynamics demonstrate close agreement with PDE benchmarks and quantify convergence via Jensen–Shannon divergence as the number of simulations increases. The approach enables accurate, open, spatially varying reservoir modeling in PBRD, with potential applications in systems biology and multi-scale simulations.

Abstract

Open biochemical systems of interacting molecules are ubiquitous in life-related processes. However, established computational methodologies, like molecular dynamics, are still mostly constrained to closed systems and timescales too small to be relevant for life processes. Alternatively, particle-based reaction-diffusion models are currently the most accurate and computationally feasible approach at these scales. Their efficiency lies in modeling entire molecules as particles that can diffuse and interact with each other. In this work, we develop modeling and numerical schemes for particle-based reaction-diffusion in an open setting, where the reservoirs are mediated by reaction-diffusion PDEs. We derive two important theoretical results. The first one is the mean-field for open systems of diffusing particles; the second one is the mean-field for a particle-based reaction-diffusion system with second-order reactions. We employ these two results to develop a numerical scheme that consistently couples particle-based reaction-diffusion processes with reaction-diffusion PDEs. This allows modeling open biochemical systems in contact with reservoirs that are time-dependent and spatially inhomogeneous, as in many relevant real-world applications.

Figures (9)

• Figure 1: Some models of reaction-diffusion processes organized by their spatial scaling and the number of particles. Only the most relevant models for this work are shown here.
• Figure 2: Illustration of some of the possible reactions in a particle-based reaction-diffusion simulation. a. Diagram of the Doi model for bimolecular reactions. When the two particles are closer than a distance of $\sigma$, they react with rate $\alpha$. Note it is not the same as the macroscopic rate $\kappa$. b. An example of a unimolecular reaction (first-order) reaction, where $A$ simply transforms into $B$. c. The backward reaction of the binding given by the Doi model. The products should be placed uniformly at a distance $\delta r$ such that $\delta r \leq \sigma$ to satisfy detailed balance frohner2018reversible.
• Figure 3: Diagram of the grand canonical master equation for an open system. It allows for an arbitrary number of diffusing particles, and it is coupled to a material reservoir. This is the result of discretizing the one-dimensional particle-based diffusion processes in contact with a constant concentration reservoir.
• Figure 4: Diagram summarizing the two main results of Section \ref{['sec:coupling']}. a. Mean field limit of the particle-based diffusion open system. If the injection rate of particles from the reservoir is set to $\gamma = D/\delta x^2$ in the master equation, the mean field yields a constant concentration boundary condition. b. Mean field limit of PBRD. It relates the miscroscopic parameters of the Doi model, $\alpha$ and $\sigma$, with the macroscopic reaction rate $\kappa$.
• Figure 5: Illustrations of the boundary coupling in the hybrid scheme. The domain is divided by an interface into the particle and the concentration domain. a. Boundary coupling for one diffusing species. In the particle domain, particles diffuse freely following Brownian motion. If a particle diffuses into the concentration domain (the reservoir), it is eliminated. Along the boundary cells in the concentration domain, we convert the concentration into particle number, generally a non-integer value. The integer part is the number of virtual particles, each can jump with rate $\gamma$ into the boundary cells in the particle domain. The fractional part corresponds to a fractional virtual particle, whose jump rate is scaled by this fraction. Note virtual particles are only drawn for illustration purposes, and they do not have a specific position within the boundary cell. The same procedure applies for multiple species with up to unimolecular reactions. b. Boundary coupling for a system with three species $(A,B,C)$ involving a bimolecular reaction $A+B\rightarrow C$. The coupling is analogous to the one in Fig. \ref{['fig:HybridSetting']}a and follows the particle-based dynamics described in Section \ref{['sec:PBA']}. If a red particle (A) is close enough to a blue one (B), they can react, and the product (C) is placed in the average position between the two particles. Unlike the coupling from Fig. \ref{['fig:HybridSetting']}a, the positions of the virtual particles are sampled uniformly within the boundary cell. This allows for bimolecular reactions to occur within all boundary cells across the coupling boundary, which makes the coupling accurate and robust. If the position of the reaction product is within the concentration domain, it is eliminated. This coupling can be applied to a general system with an arbitrary number of species with up to second-order reactions.