On the formalization of Asynchronous First Passage Algorithms

TL;DR

This paper revisits the formalization of First Passage algorithms, identifies a conceptual contradiction, and proposes a new formal framework that clarifies the principles without affecting numerical results.

Contribution

It introduces a revised formalization of First Passage algorithms that resolves a conceptual contradiction and strengthens the theoretical foundation.

Findings

Numerical results remain unaffected despite the contradiction.

A new formal framework clarifies the algorithm's principles.

Numerical evidence supports the revised formalization.

Abstract

The formalization of First Passage schemes is revisited and the emerging of a conceptual contradiction is underlined. We then show why, despite such a contradiction, the numerical results are not explicitly affected. Through a different formalization of the problem, we recast the current principles of the algorithm in a more solid conceptual framework, and numerical evidence gives further justification to our claims.

Equations

$$q_{\Omega }\left(t|{\vec{r}}_0,t_0\right)=\frac{d\left(1-S_{\Omega }\left(t|{\vec{r}}_0,t_0\right)\right)}{dt}.$$

$$g_{\Omega }(\vec{r},t|\tau ;0,0)=\frac{q_{\Omega }(\tau |\vec{r},t)p_{\Omega }(\vec{r},t|0,0)}{q_{\Omega }(\tau |0,0)}.$$

$$g_{\Omega }(\vec{r},t|0,0)=\frac{p_{\Omega }(\vec{r},t|0,0)}{S_{\Omega }(t|0,0)}.$$

$$f(\vec{r},t|0,0)=g(\vec{r},t|0,0;\vec{r}\in \Omega \forall t^{\prime }\le t)h(\vec{r}\in \Omega \forall t^{\prime }\le t)+(1-h(\vec{r}\in \Omega \forall t^{\prime }\le t))\mathcal{P}(\vec{r},t|\partial \Omega ),$$

$$g_{\Omega }(\vec{r},t|0,0)={\int }_t^{\infty }d\tau g_{\Omega }(\vec{r},t|\tau ;0,0)\frac{q_{\Omega }(\tau |0,0)}{S_{\Omega }(t|0,0)},$$

$$\frac{\partial }{\partial t}f(\vec{r},t)=D\Delta f(\vec{r},t),$$

$$f(\vec{r},t)=constp(r,t).$$

$$p_{\Omega }(r=b,t)=0;\forall t.$$

$$S_{\Omega }(t)={\int }_0^{b}dr4\pi r^2{p}_{\Omega }(r,t).$$

$$q_{\Omega }\left(t|r_0,t_0\right)=\frac{d\left(1-S_{\Omega }\left(t|r_0,t_0\right)\right)}{dt}.$$

$$\frac{\partial p(r,t)}{\partial t}=D\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial p(r,t)}{\partial r}\right).$$

$$\frac{1}{Dh(t)}\frac{\partial h(t)}{\partial t}=\frac{1}{g(r)r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial g(r,t)}{\partial t}\right).$$

$$\frac{1}{Dh(t)}\frac{\partial h(t)}{\partial t}=\frac{1}{g(r)r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial g(r,t)}{\partial t}\right)=-{\lambda }^2.$$

$$h(t)=A{e}^{-{\lambda }^{2}Dt}.$$

$$g(r)=\frac{1}{r}\left(B\sin (\lambda r)+C\cos (\lambda r)\right).$$

$$p_{\Omega }(r,t)=\frac{1}{r}\sum _m{A}_m{e}^{-\frac{m^2{\pi }^2}{b^2}Dt}\sin \left(\frac{m\pi }{b}r\right).$$

$$\begin{array}{rl}{p}_{\Omega }& (r,t|r_0,t_0)=\frac{1}{2\pi b}\frac{1}{r{r}_0}\cdot \\ & \cdot \sum _m{e}^{-\frac{m^2{\pi }^2}{b^2}D(t-t_0)}\sin \left(\frac{m\pi }{b}r\right)\sin \left(\frac{m\pi }{b}{r}_0\right).\end{array}$$

$$\begin{array}{rl}{S}_{\Omega }& (t|r_0,t_0)=-\frac{2b}{\pi r_0}\cdot \\ & \cdot \sum _m\frac{(-1{)}^m}{m}{e}^{-\frac{m^2{\pi }^2}{b^2}D(t-t_0)}\sin \left(\frac{m\pi }{b}{r}_0\right).\end{array}$$

$$\begin{array}{rl}{q}_{\Omega }& (\tau |r_0,t_0)=-\frac{2\pi D}{b{r}_0}\cdot \\ & \cdot \sum _{m}m(-1{)}^m{e}^{-\frac{m^2{\pi }^2}{b^2}D(t-t_0)}\sin \left(\frac{m\pi }{b}{r}_0\right).\end{array}$$

Full text

††preprint: AIP/123-QED

Introduction

First passage algorithms have been widely used to efficiently simulate the motion of random walkers in a lattice or Brownian particles in a continuum. Applications, across different subjects and disciplines, span from the determination of the ground state of bosonic particles KLV to the clustering of the oncoprotein Ras eGFRDras (further representative examples can be found in Refs.FPKMC ; FPKMC2 ; FPKMC3 ; AccKMC ; driftFPKMC ; nucleationFPKMC ; spatAnnFPKMC ; eGFRD ; eGFRDalldim ; eGFRDtf ; MD-GFRDanis ). Although the algorithm has been successfully used for multiple numerical applications, a proof of conceptual correctness with a solid formalism has not been treated in full yet.

The first passage algorithms allow to efficiently simulate the free diffusion of particles assuming that the motion of free particles is propagated with large spatio-temporal jumps, sampled by appropriate first-passage Green’s functions. The volume of the system is decomposed into non overlapping domains, each containing one or two particles and in every domain the underlying stochastic dynamics is tractable analytically. The relevant consequence is that the first exit-time of a single particle or of the center of mass of a particles pair can be treated exactly. Particularly effective at low molar concentrations, this approach can improve the computational performance up to several orders of magnitude with respect to a time-driven scheme MD-GFRD ; MD-GFRDanis ; MD-GFRDnew ; FPKMC ; FPKMC2 ; FPKMC3 ; eGFRDalldim . The algorithm was originally conceived in Ref.KLV as a synchronous scheme: (i) protective domains are constructed around each particle and the corresponding exit-time is sampled; (ii) the shortest exit-time is determined and particles are propagated to this time; (iii) finally, new protective domains are constructed for all particles. Differently, the first passage kinetic Monte Carlo algorithm (FPKMC) FPKMC is asynchronous (see Ref. asynch for a more comprehensive discussion about asynchronous algorithms). In this scheme, the position of the particle involved in the next exit-time event is updated locally, while the next events of all other particles are maintained in a time ordered event-list, which is updated every time a new event is sampled. The event-list must be updated also for some unscheduled events, as, indeed, when one particle is in proximity of another domain. In this case, in fact, it is impossible to construct a domain of reasonable size without overlapping with other domains. Then, the domain is burst, which means that the particle position is updated inside the domain and the exit-time previously sampled is removed from the event-list. More recently, the Molecular Dynamics-Green’s Function Reaction Dynamics (MD-GFRD) algorithm MD-GFRD ; MD-GFRDanis ; MD-GFRDnew has been developed. This is an inherently multi-scale algorithm that pairs the event-driven first exit-time sampling to a time-driven integration of the particles motion. The particle interactions are locally integrated with an inter-particle potential, while the free diffusion is updated with a first-passage sampling. For simplicity, the above-mentioned schemes are called first passage algorithms.

The event-based particle propagation is performed either by placing the particle randomly on the domain borders when the pre-sampled exit-time is reached, or by prematurely updating the particle position inside the domain when this is burst before the exit-time. The probability distribution that has been used when a domain is burst, is simply given by the diffusion propagator with the imposed condition that the particle has never reached the domain boundaries before the bursting-time FPKMC ; MD-GFRD . The position sampling is performed after the extraction of the first exit-time from the domain, but the information of the exact moment when the particle hits the boundary is missing as a condition for the probability distribution, although this would be required Redner2 . Indeed, when an observable is sampled at a time between the initial time and the time of a determined event, the probability distribution should be conditional until realization of the determined event. For instance, if the particle position is sampled immediately before the exit-time, the particle is expected to be in proximity of the domain borders. Let us assume that the sampling-time $t$ is equal to the exit-time $\tau$ minus an infinitesimal time $t=\tau-dt$, if the particle position $\vec{{r}}$ is sampled at a finite distance $\Delta\vec{R}$ from the domain border $\vec{\Omega}$ then $\vec{r}=\vec{\Omega}-\Delta\vec{R}$, which leads to the evidence that the particle travels a finite distance $\Delta\vec{R}$ in an infinitesimal time $dt$. This is a contradiction, which might occur when the condition on the exit-time is missing. Despite the conceptual contradiction underlined above, numerical applications of first passage schemes have always been shown to reproduce correctly the expected results. The discussion above naturally leads to the need of a formal proof of correctness of the algorithm which can clarify the issue about the reported conceptual contradiction. In this paper, we investigate the formal construction of the algorithm and clarify the pending aspects of the previous discussion. Specifically, we show how the position probability distribution used in previous first passage schemes can be derived from the conditional probability distribution: the two distributions result statistically identical when the bursting-time is chosen randomly, as it happens in multi-particle simulations.

The distinction between the conditional and unconditional probability becomes relevant when a systematic bursting of the domains is performed. This might be the case of MD-GFRD, a scheme where the particles interaction is simulated with an inter-particle potential in a time-discretization of the overdamped Langevin equation. When this approach is used, the particles clock must be synchronized during their interaction MD-GFRDnew . This synchronization can be obtained by projecting the exit-time onto a discrete temporal grid. The projection can be performed after the exit-event with a fractional Brownian motion propagation or before the exit-event by sampling the particle position in the last discrete time before the exit-time. For the latter case, the choice of the probability distribution is crucial, for the reason that a systematic use of an unconditional probability distribution under these circumstances would bring the simulation to inaccurate results.

Results and Discussion

The first exit-time from a protective domain $\Omega$ is derived by imposing absorbing boundary conditions on the domain boundary $\partial\Omega$Redner . The Green’s function $p_{\Omega}(\vec{r},t|\vec{r}_{0},t_{0})$ is the solution to the diffusion equation when absorbing boundary conditions on $\partial\Omega$ and initial conditions in $\vec{r}_{0}$ at $t_{0}$ are imposed, and it represents the probability that the particle is located in the position $\vec{r}$ at time $t$ without having previously crossed the domain borders. The integration of the probability distribution $p_{\Omega}(\vec{r},t|\vec{r}_{0},t_{0})$ over the whole domain gives the survival probability $S_{\Omega}(t|\vec{r}_{0},t_{0})$, i.e. the probability that the particle has not reached the domain borders at $t$ yet. The first exit-time from the domain $q_{\Omega}(t|\vec{r}_{0},t_{0})$ is derived from the survival probability

[TABLE]

It is assumed that the particle is initially located in the origin, hence the exit-time is sampled from Eq. (1) with $\vec{r}_{0}\rightarrow$0 at $t_{0}=0$.

After the extraction of the particle escape at time $\tau$, the conditional probability distribution of the particle position at an intermediate time $0<t<\tau$ is obtained from the Bayes theorem

[TABLE]

In first passage schemes the unconditional probability is used, renormalized by the survival probability

[TABLE]

For a detailed derivation of the above probability distributions, we refer to the appendix A.

Derivation of the unconditional probability

Eq. (2) and Eq. (3) are clearly different (see Fig. 1) and the assumption that these two distributions can be used interchangeably might result surprising. In particular in asynchronous schemes as FPKMC and MD-GFRD, the correctness of the use of the unconditional probability is not evident, because one may ask how the position probability of the particle inside the domain is characterized when the position of other particles is updated and known. This can be justified by considering the probability distribution conditional with respect to the exit-time. Indeed, when the exit-time is extracted, it is ensured that the particle is inside the domain and that the particle position is probabilistically described by Eq. (2) until realization of the exit-event.

It is possible to sample from Eq. (3), only when the sampling time is independent of the exit-time extraction. The problem can also be reformulated by introducing $f(\vec{r},t|0,0)$, which represents the unconditional probability that the particle is in $\overrightarrow{r}$ at $t$. The particle is supposed to be initially placed in a protective domain $\Omega$. According to the Bayes theorem, the probability distribution $f(\vec{r},t|0,0)$ can be decomposed in two conditional probabilities with respect to the event that the particle never crossed the domain borders $\partial\Omega$ until $t$

[TABLE]

where $h(\vec{r}\in\Omega\thinspace\forall t^{\prime}\leq t)$ is the probability that the particle has never escaped the domain, and $\mathcal{P}(\vec{r},t|\partial\Omega)$ is the probability that the particle is in $\vec{r}$ after crossing the domain borders $\partial\Omega$. The particle position at $t$ can be sampled from the above equation as follows:

Sample the first exit-time $\tau$. 2. 2.

if $t<\tau$: then the particle position is sampled from $g(\vec{r},t|0,0;\vec{r}\in\Omega\thinspace\forall t^{\prime}\leq t)$. 3. 3.

else if $t>\tau$: then the particle position is the result of a Wiener process of time length $t-\tau$, starting in a random position on the domain borders.

For a given observation time $t$ the first-exit time $\tau$ is extracted, the comparison between these times, i.e. the $\it{if}$ statements in the list, in effect samples $h(\vec{r}\in\Omega\thinspace\forall t^{\prime}\leq t)$. The only relevant assumption to sample $h(\vec{r}\in\Omega\thinspace\forall t^{\prime}\leq t)$ is that $t$ and $\tau$ are uncorrelated.

The probability distribution $g(\vec{r},t|0,0;\vec{r}\in\Omega\thinspace\forall t^{\prime}\leq t)$ requires that the particle is known to be inside the domain at time $t$ and that it has never left it before, conditions satisfied by both Eq. (2) and Eq. (3). The connection between these probability distributions is that Eq. (3) is given by the integral of Eq. (2) over all possible exit-times $\tau>t$ weighted with their likelihood

[TABLE]

where $q_{\Omega}(\tau|0,0)/S_{\Omega}(t|0,0)$ is the probability that the particle exits at $\tau$, conditional with respect to the information that at $t$ is still inside the domain. From Eq. (5) it is possible to infer that, given the information that the particle at $t$ is still inside the domain, the position sampled from Eq. (3) is on average the same as if sampled from Eq. (2), when the exit-time is extracted from Eq. (1). An evident assumption that has been made is that $t$ and $\tau$ are uncorrelated, otherwise the integral in Eq. (5) could not be defined. This assumption is the essential reason that explains why sampling from Eq. (3) and from Eq. (2) are statistically equivalent, see Eqs. (4) and (5).

In Fig. 2, simulations have been performed using the two discussed probability distributions to update the particles position when domains are burst. Since in these simulations the domain bursts occur only in occasion of the random collision of an external particle with the protective domain of another particle, it is reasonable to assume that the bursting-time and the first exit-time are uncorrelated, hence sampling from Eq. (3) is statistically equivalent to sampling from Eq. (2). The mean squared displacement in a simple diffusive process and the kinetics of a two species annihilation reaction-diffusion system have been observed, and, as expected, the use of the conditional and unconditional distributions led to indistinguishable results.

In asynchronous schemes as FPKMC and MD-GFRD the particle position is updated serially while other particles lie in protective domains, but it is still an open question how the motion of particles in such protective domains is characterized before their exit-time. Sampling the particles position from Eq. (3) simply gives an average behavior that is justifiable when the bursting-time is given randomly. More rigorously, the particle position should be sampled from Eq. (2) until realization of the exit-time. Let us suppose that a domain is burst at a time $t$ and the particle that was inside the domain interacts with another particle on the time grid $\left\{t_{i}\right\}_{i=0}^{n}$, where $t_{i+1}=t_{i}+dt$ and $dt$ is the integration step. The particle position should be sampled from Eq. (2) at all time steps until the pre-sampled exit time, $t>\tau$. However, in Eq. $(\ref{eq:pos_cond}$) the domain has been assumed to be interaction-free until $\tau$, but the particle is also assumed to be interacting. It is still possible to use the conditional sampling because the Langevin equation is linear, and the contribution of the free displacement can be separated from the interaction term. When the particle is burst at time $t$ all particle displacements until $\tau$ are sampled and recorded. These will give the free displacements of the particle that will then be summed to the interaction term (see Fig. 3). The sampling procedure just enlightened involves many integration steps between the bursting time and the first-exit time, and in each of these steps the particle position is sampled from Eq. (2). This is clearly computationally more expensive than sampling once from Eq. (3), that can be done efficiently using the rejection method as described in Refs. FPKMC2 ; FPKMC3

In the above derivations $t$ and $\tau$ have been assumed uncorrelated, therefore it is possible to sample from Eq. $(\ref{eq:pos_uncond}$). Although this assumption might result natural when the bursting-time comes from the random encounters between different particles in a multi-particle simulation, a systematic bursting happening at a time which is function of the exit-time could be used in a multi-scale simulation as MD-GFRD. If the particle synchronization, performed in this multi-scale algorithm, is obtained by time-projection onto the last discrete time before the exit-time of the particle, this systematic bursting would clearly insert a correlation between $t$ and $\tau$. Eq. (5) would not be valid because it is not possible to integrate $\tau$ independently from $t$ and in Eq. (4) $h(\vec{r}\in\Omega\thinspace\forall t^{\prime}\leq t)$ is not sampled.

In Fig. 4, the mean squared displacements of particles simulated as illustrated in Fig. 3 and particles that use the position sampling from Eq. (3) are shown and compared to a brute-force integration that is taken as reference. The synchronization step was performed before the exit-time during simulations, and it is clear from the figure that the use of Eq. (3) makes the particles diffuse slower. Indeed, immediately before the sampled exit-time the particle is expected to be in proximity of the domain borders, but this is not ensured without the condition on the exit-time (see Fig. 1). The correct slope of the mean square displacement is instead reproduced when the Eq. (2) is sampled. In the simulation, particles are interacting and their mean squared displacements lie below the expected free diffusion value due to crowding effects.

Conclusion

We have analyzed the conceptual consistency of First Passage Algorithms. We assess that the use of the unconditional probability is only justified under the assumption that the bursting time and the first-exit time are uncorrelated. We also show that the unconditional probability distribution can be derived from a formally correct conditional probability distribution. The use of the unconditional probability distribution routinely employed leads to formal contradictions, although this aspect does not explicitly emerge in current numerical applications. Given the increasing number of algorithms which implement asynchronous first passage schemes, we note that the enlightened contradiction might represent a pitfall in future applications. We finally give numerical evidence that confirms the solidity of our formal derivation, and thus places First Passage Algorithms on firm conceptual ground.

Appendix A Probability distributions derivation

The motion of a Brownian particle is assumed to be described stochastically by the diffusion equation

[TABLE]

where $f(\vec{r},t)$ is the probability to find the particle in position $\vec{r}$ at time $t$, and $D$ is the diffusion coefficient. We assume that the particle position is known at time $t_{0}$, and we are interested in the first exit-time of the particle from a domain $\Omega$ that is constructed around the particle. The domain $\Omega$ is supposed to be spherical and centered on the particle position. Given the spherical symmetry of the system, the angular coordinates can be averaged out from Eq. (6)

[TABLE]

The first exit-time from a protective domain $\Omega$ is derived by imposing absorbing boundary conditions on the domain boundary $\partial\Omega$Redner . The domain $\Omega$ is here assumed to be a sphere of radius $b$

[TABLE]

This condition ensures that once the particle hits the domain borders, it is removed from the system. The integral of Eq. (8) over the whole domain gives the survival probability $S_{\Omega}(t)$, i.e. the probability that, until $t$, the particle has always remained inside the domain

[TABLE]

The probability of crossing the domain borders for the first time at $t$ is derived from the survival probability

[TABLE]

In this paper, our aim is to characterize the probability distribution of the particle position before it escapes the domain. Firstly, the Green’s function for the diffusion equation is found, then exit-time boundary conditions are applied to the solution. The diffusion equation is only solved for the radial component

[TABLE]

The solution is found by assuming that it can be separated in the variables $t$ and $r$, $p(r,t)=g(r)h(t)$,

[TABLE]

Each term of the above equation has only one variable, in this case both terms can be equated to a constant value resulting in two different equations

[TABLE]

The solution of the first equation is an exponential function

[TABLE]

The second equation can be expressed as a harmonic oscillator differential equation by substituting $\tilde{g}(r)=r\thinspace g(r)$. The solution is then straightforward

[TABLE]

The capital letters $A,B,C$ used are normalization factors. Absorbing boundary conditions $p_{\Omega}(b,t)=0$ are imposed on the equation $p(r,t)=g(r)h(t)$ obtained from Eqs. (14) and (15)

[TABLE]

The particle is known to be in $r_{0}$ at $t_{0}$

[TABLE]

The survival probability is obtained as shown in Eq. (9)

[TABLE]

The first exit-time is then derived from the survival probability (see Eq. (10))

[TABLE]