Accurate stochastic simulation of nonlinear reactions between closest particles

TL;DR

The paper tackles the challenge of proximity-based trimolecular reactions in diffusive particle systems by deriving a nonlinear PDE for the density of the closest C molecule and a leading-order steady-state solution $P_0(r_1,r_2)$. It develops a mixed-primal finite-element framework to compute higher-order corrections ($P_h$, $\boldsymbol{s}_h$, $Q_h$) on a truncated axisymmetric domain and demonstrates how the commonly used leading-order rate $K_1(c)$ underestimates the true rate due to neglected flux in the $\boldsymbol{\hat{r}}_2$ direction. By quantifying these $O(\sigma_{\max}^2)$ corrections with FEM and then optimizing the reactive boundary $\partial \Omega_R$ via Levenberg–Marquardt to match a target rate $K_D(c)$, the authors provide a practical, analytic-looking boundary correction that improves the realism of particle-based simulations. The work clarifies the sources of discrepancy in proximity-based kinetics and offers a systematic route to more accurate trimolecular reaction modeling, with potential extensions to more general reactive boundaries and concentration regimes.

Abstract

We study a system of diffusing point particles in which any triplet of particles reacts and is removed from the system when the relative proximity of the constituent particles satisfies a predefined condition. Proximity-based reaction conditions of this kind are commonly used in particle-based simulations of chemical kinetics to mimic bimolecular reactions, those involving just two reactants, and have been extensively studied. The rate at which particles react within the system is determined by the reaction condition and particulate diffusion. In the bimolecular case, analytic relations exist between the reaction rate and the distance at which particles react allowing modellers to tune the rate of the reaction within their simulations by simply altering the reaction condition. However, generalising proximity-based reaction conditions to trimolecular reactions, those involving three particles, is more complicated because it requires understanding the distribution of the closest diffusing particle to a point in the vicinity of a spatially dependent absorbing boundary condition. We find that in this case the evolution of the system is described by a nonlinear partial integro-differential equation with no known analytic solution, which makes it difficult to relate the reaction rate to the reaction condition. To resolve this, we use singular perturbation theory to obtain a leading-order solution and show how to derive an approximate expression for the reaction rate. We then use finite element methods to quantify the higher-order corrections to this solution and the reaction rate, which are difficult to obtain analytically. Leveraging the insights gathered from this analysis, we demonstrate how to correct for the errors that arise from adopting the approximate expression for the reaction rate, enabling for the construction of more accurate particle-based simulations than previously possible.

Figures (9)

• Figure 2.1: The diffusive Jacobi coordinates for a system containing three molecules. The coordinate $\boldsymbol{\eta}_{1}$ describes the separation between the first two molecules, labelled $A$ and $B$, while $\boldsymbol{\eta}_{2}$ describes the separation of the third molecule $C$ from the centre of diffusion of $A$ and $B$ which we denote $\bar{\boldsymbol{x}}_{1}$ and is defined in Equation \ref{['eq:x_bari_def']}. Finally, $\boldsymbol{\eta}_{0}$ is the centre of diffusion of the three molecules ($\bar{\boldsymbol{x}}_{2}$) and changes in this coordinate correspond to translations of the state.
• Figure 2.2: The state space for a system that contains a single molecule of $A$ (the blue point), a single molecule of $B$ (the green point) and five molecules of $C$ (the orange points). Each combination of $A$ and $B$ with a particular molecule of $C$ gives a point within the state space on the right, shown by the orange points. The position of the $B$ molecule relative to the $A$ molecule defines a manifold of constant $\boldsymbol{\eta}_{1}$ that all states lie on. This manifold diffuses in $\boldsymbol{\eta}_{1}$ in accordance with Equation \ref{['eq:diffusion_eta1']} and the states on the manifold diffuse independently in $\boldsymbol{\eta}_{2}$ according to Equation \ref{['eq:diffusion_eta2']}. The first state that crosses the inner boundary, depicted by the blue curve and defined in Equation \ref{['eq:generalised_inner_BC']}, will be absorbed.
• Figure 3.1: Approximate solutions for the convergence test, computed with the lowest-order finite element scheme and shown on a coarse mesh.
• Figure 4.1: The (scaled) difference $\Delta P\left(r_1,r_2\right)$, defined in Equation \ref{['eq:scaled_difference']}, between the steady-state finite element solution $P_h\left(r_1,r_2\right)$ and the leading-order solution $P_0(r_1,r_2)$ for $\hat{D}_1 = 2$, $\hat{D}_2 = 1.5$, $\sigma_{\text{max}} = 0.1$, $\Gamma = 1$, and $c=10$. The leading-order solution underestimates the probability, resulting in the yellow region where $\Delta P\left(r_1,r_2\right) > 0$, that the closest molecule of $C$ is near to the origin and instead favours larger values of $r_2$, which results in the dark blue region where $\Delta P\left(r_1,r_2\right) < 0$. This discrepancy arises because $P_0$ neglects the advection towards the origin in the $\boldsymbol{\hat{r}}_{2}$ (horizontal) direction. The difference appears most pronounced near the bottom boundary $\partial \Omega_R$, defined in Equation \ref{['eq:exp_boundary']}, but this is because the solutions decay with increasing $r_1$. The surface plot in (a) shows the finite element domain from Equation \ref{['eq:finite_ele_domain']}, while the plot in (b) shows the subdomain $\left\{(r_1, r_2) : \sigma_{\text{max}}\text{exp}\left(\frac{-4\pi \Gamma r_2^3}{3}\right) \leq r_1 \leq 1,0 \leq r_2 \leq 1\right\}$.
• Figure 4.2: The corrections to the flux of the probability density that are neglected by the approximation $K_1(c)$ to the steady-state reaction rate defined in Equation \ref{['eq:closest_rate_old']}, for $\hat{D}_1 = 2$, $\hat{D}_2 = 1.5$, $\sigma_{\text{max}} = 0.1$, $\Gamma = 1$, and $c=10$. The surface plot in (a) shows the (scaled) correction $\Delta s_1(r_1, r_2)$ to the diffusive flux in the $\boldsymbol{\hat{r}}_{1}$ (vertical) direction defined by Equation \ref{['eq:r1_flux_correction']}. The diffusive flux of the leading-order solution $P_0(r_1,r_2)$ underestimates the magnitude of the $\boldsymbol{\hat{r}}_{1}$ component of the finite element solution $\boldsymbol{s}_h\left(r_1,r_2\right)$ for values of $r_2$ less than approximately $0.25$, resulting in the dark blue region where $\Delta s_1(r_1, r_2) < 0$. In contrast, the magnitude of this flux is overestimated for larger values of $r_2$, resulting in the yellow region where $\Delta s_1(r_1, r_2) > 0$. Similarly, the surface plot in (b) shows the (scaled) correction $\Delta s_2(r_1, r_2)$ to the flux accounted for by $K_1(c)$ in the $\boldsymbol{\hat{r}}_{2}$ (horizontal) direction defined by Equation \ref{['eq:r2_flux_correction']}. Since $K_1(c)$ entirely neglects the flux in this direction over $\partial \Omega_R$, as it is an $O(\sigma_{\text{max}}^2)$ contribution, this correction is simply the $\boldsymbol{\hat{r}}_{2}$ component of the finite element solution $\boldsymbol{s}_h\left(r_1,r_2\right)$. Both plots display the subdomain $\left\{(r_1, r_2) : 0 \leq r_1 \leq 1,0 \leq r_2 \leq 1\right\}$.