An Unstructured Mesh Reaction-Drift-Diffusion Master Equation with Reversible Reactions

TL;DR

The paper introduces the Convergent Reaction-Drift-Diffusion Master Equation (CRDDME), a framework that discretizes drift-diffusion with a background potential on unstructured meshes while preserving detailed balance for reversible reactions, and demonstrates convergence to the continuum Volume Reactivity (VR) model. It combines an Edge-Average Finite Element (EAFE) discretization for drift-diffusion with a finite-volume treatment of reaction terms to yield a continuous-time, discrete-space jump process whose equilibrium statistics respect Gibbs-Boltzmann distributions. A discrete detailed-balance condition is enforced to define reaction rates consistently at the mesh level, enabling accurate and mesh-refinement–stable simulations of A+B $\leftrightarrows$ C in complex geometries, including multiparticle systems. Numerical experiments, including annihilation, reversible reactions, multiparticle dynamics, and TCR-pMHC transport in immune synapse formation, validate convergence (approximately second order) and illustrate the method’s relevance to cell biology and receptor dynamics. The work provides a flexible, convergent alternative to RDME methods for drift-diffusion with potentials and lays groundwork for extensions to time-dependent potentials and two-body interactions.

Abstract

We develop a convergent reaction-drift-diffusion master equation (CRDDME) to facilitate the study of reaction processes in which spatial transport is influenced by drift due to one-body potential fields within general domain geometries. The generalized CRDDME is obtained through two steps. We first derive an unstructured grid jump process approximation for reversible diffusions, enabling the simulation of drift-diffusion processes where the drift arises due to a conservative field that biases particle motion. Leveraging the Edge-Averaged Finite Element method, our approach preserves detailed balance of drift-diffusion fluxes at equilibrium, and preserves an equilibrium Gibbs-Boltzmann distribution for particles undergoing drift-diffusion on the unstructured mesh. We next formulate a spatially-continuous volume reactivity particle-based reaction-drift-diffusion model for reversible reactions of the form $\textrm{A} + \textrm{B} \leftrightarrow \textrm{C}$. A finite volume discretization is used to generate jump process approximations to reaction terms in this model. The discretization is developed to ensure the combined reaction-drift-diffusion jump process approximation is consistent with detailed balance of reaction fluxes holding at equilibrium, along with supporting a discrete version of the continuous equilibrium state. The new CRDDME model represents a continuous-time discrete-space jump process approximation to the underlying volume reactivity model. We demonstrate the convergence and accuracy of the new CRDDME through a number of numerical examples, and illustrate its use on an idealized model for membrane protein receptor dynamics in T cell signaling.

Figures (10)

• Figure 1: Our approach to approximating the stochastic reaction-drift-diffusion process of the Volume Reactivity (VR) model (upper left) by a spatially-discrete jump process (lower left). Derivation proceeds clockwise from the former to the latter via the CRDDME we develop in this work.
• Figure 2: Survival time distributions vs. $t$ from the two-particle $\textrm{A} + \textrm{B} \to \varnothing$ reaction in a square in panel A and in a disk in panel B. In both cases each curve was estimated from 100000 simulations using the ecdf command in MATLAB. The legends give the ratio, $\varepsilon$/$h$ as the mesh is refined ($h$ is approximately successively halved). See Section \ref{['sec:annihilRx']} for parameter values.
• Figure 3: Convergence of the mean reaction time $\mathbb{E}[T_{\textrm{bind}}]$. In panel A we plot the mean reaction time $\mathbb{E}[T_{\textrm{bind}}]$ vs. $\varepsilon$/$h$ as $h$ is (approximately) successively halved. Each mean reaction time for the two cases was estimated from 100000 simulations. 95$\%$ confidence intervals are drawn on each data point, but for some points are smaller than the marker labeling the point. See Section \ref{['sec:annihilRx']} for parameter values. In panel B we demonstrate the rate of convergence when using a square and a disk domain by plotting the difference between successive points on the $\mathbb{E}[T_{\textrm{bind}}]$ vs $\varepsilon/h$ curves from panel A. The smaller of the two $h$ values is used for labeling. The effective convergence rate to zero of the FPE with both reaction mechanisms scales roughly like $O(h^2)$. Note that the search time in the square is smaller in part as the area of the circle is a factor of $\pi$ larger than the square.
• Figure 4: Probability that molecules are in bound state $P_{\textrm{bound}}(t)$ versus time for CRDDME SSA simulations and BD simulations. Dark purple curves correspond to $P_{\textrm{bound}}(t)$ from CRDDME simulations and the light blue to $P_{\textrm{bound}}(t)$ from BD simulations. The magenta solid line shows the steady-state $P_{\textrm{bound}}$ value for $K_{\textrm{d}} = 2$. $95\%$ confidence intervals for the CRDDME simulations are drawn with red error bars, and as a translucent blue ribbon for the BD simulations. Panel B shows a zoomed in version for a portion of panel A when the system reaches steady states, making clear the scale of the 95% confidence intervals. Each curve was estimated from 100000 simulations by dividing the number of simulations in which the system is in the bound state at a given time by the total number of simulations. The domain $\Omega$ is a circle centered at $(0.05, 0.05)$ with radius $r = 0.1\mu$m and was discretized into $44945$ polygons.
• Figure 5: Convergence of the probability the molecules are in the bound state, $P_{\textrm{bound}}(t)$, as $h \to 0$. In panel A we plot $P_{\textrm{bound}}(t)$ vs. $\varepsilon$/$h$ as $h$ is (approximately) successively halved. Each curve was estimated from 100000 simulations. We can observe convergence as the mesh width $h$ approaches $0$. See Section \ref{['sec:revRx']} for parameter values. In panel B we demonstrate the rate of convergence by plotting the difference between successive points on the $P_{\textrm{bound}}(t)$ vs $\varepsilon/h$ curves from panel A at $t = 0.1$s. The smaller of the two $h$ values is used for labeling. The effective convergence rate to zero of the CRDDME with both reaction mechanisms scales roughly like $O(h^2)$.