Open reaction-diffusion systems: bridging probabilistic theory and simulations across scales

TL;DR

This work develops a rigorous, cross-scale framework for open reaction-diffusion systems by extending the Chemical Diffusion Master Equation (CDME) to incorporate macroscopic reservoirs via boundary-layer fluxes modeled as creation/degradation events. By taking expectations of the CDME, the authors derive macroscopic PDEs for linear and nonlinear reactions and recover diffusion-influenced models, including Smoluchowski-type boundary conditions, establishing a principled link between particle-level stochastic dynamics and continuum descriptions. They further introduce and validate a suite of numerical schemes—tau-leap, Gillespie exact, and explicit methods—for reservoir interactions, including implicit boundary treatments—enabling consistent, efficient hybrid simulations such as open diffusion with reactive boundaries and a hybrid SIR model. The results provide a robust multiscale workflow to simulate open reaction-diffusion processes across scales, with direct implications for biochemical, epidemiological, and social-agent systems where environment exchange is essential.

Abstract

Reaction-diffusion processes are the foundational model for a diverse range of complex systems, ranging from biochemical reactions to social agent-based phenomena. The underlying dynamics of these systems occur at the individual particle/agent level, and in realistic applications, they often display interaction with their environment through energy or material exchange with a reservoir. This requires intricate mathematical considerations, especially in the case of material exchange since the varying number of particles/agents results in ``on-the-fly'' modification of the system dimension. In this work, we first overview the probabilistic description of reaction-diffusion processes at the particle level, which readily handles varying number of particles. We then extend this model to consistently incorporate interactions with macroscopic material reservoirs. Based on the resulting expressions, we bridge the probabilistic description with macroscopic concentration-based descriptions for linear and nonlinear reaction-diffusion systems, as well as for an archetypal open reaction-diffusion system. Using these mathematical bridges across scales, we finally develop numerical schemes for open reaction-diffusion systems, which we implement in two illustrative examples. This work establishes a methodological workflow to bridge particle-based probabilistic descriptions with macroscopic concentration-based descriptions of reaction-diffusion in open settings, laying the foundations for a multiscale theoretical framework upon which to construct theory and simulation schemes that are consistent across scales.

Figures (11)

• Figure 1: Structure of the phase space of the CDME. a. Phase space for a system with one chemical species $A$, consisting of discrete sets of continuous diffusion domains $\mathbb{X}$ depending on the number of particles. The transitions between sets depend on the specific reaction system. The arrows only show possible transitions between first neighbors. b. Analogous phase space but for a system with two chemical species $A$ and $B$. As particles of the same species are statistically indistinguishable, the ordering in phase space is not relevant.
• Figure 2: Diagram representing loss from and gain into the $n$ particle state due to the reaction $kA\rightarrow lA$ assuming $k>l$. The loss $\mathcal{L}_n \rho_n$ must depend on the reactant's positions, so it is a function of $x^{(n)}$, while the gain $\mathcal{G}_n\rho_{n+k-l}$ depends on the positions of the products, and thus must also be a function of $x^{(n)}$.
• Figure 3: Schematic calculation of how to obtain the loss and gain terms of the CDME for degradation and creation reactions for the $n$-particle configuration. In all cases, it follows the same procedure: 1. calculate the loss or gain due to one possible reaction in terms of the rate function. 2. calculate the total loss due to any possible reaction. Note the rate function for the degradation case is a function of the reactant's position, while for the creation is a function of the product's position. In general, regardless of the reaction, both the loss and gain are functions of the particle configuration of the current state $x^{(n)}$. For the the loss, this dependence manifests in terms of the reactant's positions, while the product's positions are integrated out (if any). Conversely, for the gain, this manifests as dependence in terms of the product's positions, while the reactant's positions are integrated out(if any).
• Figure 4: Discretizations diagrams. a. Illustration of the interaction between the system and the reservoir assuming the use of discrete rates from \ref{['eq:generalRates']}. The boundary layer $\Omega$ is shown as well as its mirror region in the reservoir domain $\hat{\Omega}$. b. Diagram of one possible discretization of the Dirac delta and its derivative as a function of a small parameter $\epsilon$. The discretization of the Dirac delta is denoted as $\delta_\epsilon(x)$, and its discretized derivative as $\delta'_\epsilon(x)$.
• Figure 5: Diagram of connections between the CDME and the reaction-diffusion PDEs (RD-PDEs) for three examples, as well as the connection between the parameters at the different scales. The PDE descriptions are recovered by taking the mean-field. Note that for linear systems the mean-field is equivalent to the large copy number limit del2018grandkostre2021coupling. a. For simple linear reactions, like creation and degradation, the parameters at the particle level match those at the concentration description. b. For nonlinear reactions involving two reactants, such as mutual annihilation and bimolecular reactions, we choose the volume reactivity model for the rate function doi1976stochastic. This depends on a microscopic rate $\lambda_0$ and reactive distance $\sigma$. Our work yields the macroscopic rate in the RD-PDE in terms of the microscopic parameters. Note in this case one also requires to assume large copy numbers to neglect the covariance arising from the nonlinear reaction kostre2021coupling. c. Although the reaction modeled is in principle nonlinear, the model focuses on the concentration of $B$ modulated by one reservoir on one edge and a reactive boundary on the other, so the particle dynamics are effectively linear, i.e. there are no pair interactions. Our approach shows how to consistently choose the microscopic rate functions to match the macroscopic description of the open system.