Reactive Boundary Conditions as Limits of Interaction Potentials for Brownian and Langevin Dynamics

A popular approach to modeling bimolecular reactions between diffusing molecules is through the use of reactive boundary conditions. One common model is the Smoluchowski partial absorption condition, which uses a Robin boundary condition in the separation coordinate between two possible reactants. This boundary condition can be interpreted as an idealization of a reactive interaction potential model, in which a potential barrier must be surmounted before reactions can occur. In this work we show how the reactive boundary condition arises as the limit of an interaction potential encoding a steep barrier within a shrinking region in the particle separation, where molecules react instantly upon reaching the peak of the barrier. The limiting boundary condition is derived by the method of matched asymptotic expansions, and shown to depend critically on the relative rate of increase of the barrier height as the width of the potential is decreased. Limiting boundary conditions for the same interaction potential in both the overdamped Fokker-Planck equation (Brownian Dynamics), and the Kramers equation (Langevin Dynamics) are investigated. It is shown that different scalings are required in the two models to recover reactive boundary conditions that are consistent in the high friction limit where the Kramers equation solution converges to the solution of the Fokker-Planck equation.

(1. 2) In the special case that ϕ ≡ 0 the molecule simply moves by Brownian motion, 3) is a popular description for the movement of molecules within cells. It has been used to model spatial transport in a number of computational packages for simulating intracellular processes, including Smoldyn [2,3,30], GFRD [34,35] and FPKMC [27,28]. A common approach for then coupling molecular interactions (diffusion-limited reactions) to (1.3) in these packages is to postulate that reactions can occur by one of several possible mechanisms if the corresponding reactants are sufficiently close [1,16]. The LD model (1.1) with ϕ ≡ 0 provides a more microscopic description of diffusion than the BD model (1.3). It computes both position and velocity of a molecule by assuming that the molecule is subject to a normally distributed random force during each time increment. In particular, LD can be considered as an intermediate description between detailed molecular dynamics (MD) simulations and BD simulators [12]. Typical full-atom MD simulations use time steps of the order of 1 fs = 10 −15 s [24], while Smoldyn discretizes the equation (1.3) with times steps ranging from nanoseconds to milliseconds, depending on a particular application [30]. Stochastic descriptions which compute both position and velocity of diffusing particles, including LD, are applicable on intermediate time scales [12,13].
One advantage of Smoldyn or similar BD packages is that they can simulate wholecell dynamics. BD models based on equation (1.3) have been applied to a number of biological systems including signal transduction in E. coli [26], actin dynamics in filopodia [18], the MAPK pathway [34] and intracellular calcium dynamics [10]. In these applications, the positions of all diffusing molecules are updated according to (1.3) and the distances between each pair of possible reactants (for bimolecular reactions) are calculated. Each reaction then occurs (with a given probability) when the computed distance is smaller than a specified reaction radius (as in Smoldyn), or alternatively occurs when the computed distance exactly equals a specified reactionradius (as in GFRD and FPKMC). To our knowledge, there is no established spatiotemporal simulator of intracellular processes based on the LD model. In order to develop one, one has to first investigate how bimolecular reactions might be described in the LD context.
One possible way to implement bimolecular reactions in LD is to adopt the same approach as in BD. That is, the positions and velocities of a diffusing molecule would evolve according to the LD model (1.1) (with ϕ ≡ 0) and each bimolecular reaction would occur (with a given probability) if the distance between two reactants is smaller than the reaction radius. However, as normally formulated this description of bimolecular interactions would not make use of a molecule's velocity (as is available in LD). That is, the LD bimolecular reaction model would not provide any more physical detail than a BD model.
In this work we step back from the normal BD bimolecular reaction model to the more microscopic reaction mechanism of a molecular interaction potential. The general potential forms we consider represent an irreversible bimolecular reaction as a molecule surmounting a steep potential barrier in the separation coordinate from a stationary target molecule, after which it enters an infinitely deep well (the "bound" state). We show that the popular Smoluchowski partial-adsorption BD reaction model [22] can be derived in the simultaneous limit that the width of the barrier approaches zero and the height of the barrier becomes infinite. Using the same potential model, we then examine what limiting reactive mechanism arises in the corresponding LD model, obtaining a specular reflection boundary condition. We conclude by showing that in the high friction limit, β → ∞, the LD model with the specular reflection bimolecular reaction model converges back to the Smoluchowski partial-adsorption BD reaction model, consistent with the kinetic boundary layer studies of [7,8,23]. Our results demonstrate how an interaction potential model for a bimolecular reaction can be parametrized in either LD or BD models to be consistent with a BD model based on a partial-adsorption reaction mechanism.
2. Problem Setup. We consider the movement of a molecule that can undergo an irreversible bimolecular reaction with a second stationary molecule, hereafter called the reactive target, which is modelled as a sphere of radius r b . Let B r (0) ⊂ R d denote the d-dimensional ball of radius r centered at the origin. Then the reactive target is given as ∂B r b (0). We assume that the diffusing molecule moves within the ddimensional domain Ω which satisfies defines Ω as the domain which is the exterior to the reactive target (sphere ∂B r b (0)) and which includes all points which have a distance less than L 1 − r b > 0 from the reactive target. A simple example of Ω satisfying (2.1), for which we can find some explicit solutions in Section 3.1, is given as This is a standard ansatz for deriving formulae for the reaction radius in BD descriptions, e.g. in the Smoluchowski or Doi models [1,16,32]. In most cases, we will further assume that Ω is bounded with smooth boundary, i.e. equation (2.1) is altered by a bound from above A spherically symmetric example of domain Ω satisfying condition (2.3) is given as We assume that the molecule is adsorbed instantly upon reaching the surface of the reactive target, so that In order to reach the surface, the molecule has to overcome a potential barrier, which, denoting r = |x|, is given by where ε > 0 and ϕ > 0 are positive constants, and ψ : [0, 1] → [0, 1] is a smooth function with a (unique) global maximum at ψ(0) = 1, and ψ(1) = 0. These imply We assume that ε in (2.6) is chosen sufficiently small so that ε (L 1 − r b ) for r b and L 1 given in assumption (2.1).
Instead of studying the Langevin equation (1.1) we shall work with the corresponding equation for the probability density that X(t) = x and V (t) = v, denoted p(x, v, t), which satisfies the Kramers equation (also called the phase-space Fokker-Planck equation) where ∇ x (resp. ∇ v ) denotes the gradient in x (resp. v) variable. Considering the overdamped limit (1.2), the corresponding equation for the probability distribution p(x, t) is given as the Fokker-Planck equation In this paper, we show that equations (2.7) and (2.8) with fully adsorbing boundary condition (2.5) on ∂B r b (0) are in suitable limits equivalent to a diffusion process with partially adsorbing (reactive, Robin) boundary condition on ∂B r b (0), namely to the problem where K is a suitable Robin boundary constant (reactivity of the boundary) and x/r b is a unit normal vector to sphere ∂B r b (0) at point x. Considering spherically symmetric Ω given by (2.2) or (2.4), and a spherically symmetric initial condition, we can rewrite the Fokker-Planck equation (2.8) as where p(r, t) is the radial distribution function. Then the limiting Robin boundary problem (2.9)-(2.10) is given as The Robin boundary constant K in equations (2.10) and (2.12) is (together with D) an experimentally determinable (macroscopic) parameter. In this way, we are able to parametrize both BD and LD models using experimental data.
The rest of this paper is organized as follows. In Section 3, we investigate the BD description (1.2). To get some insights into this problem, we first consider the specific case of a linear interaction potential and spherically symmetric domain (2.2) for d = 3. In this case, we can explicitly solve the corresponding Fokker-Planck equation as shown in Section 3.1, and prove the convergence of this solution, as ε → 0 and ϕ → ∞, to the solution of a model involving a reactive Robin boundary condition. We then continue in Section 3.2 with an asymptotic analysis, in the same dual limit, of the full BD model (1.2) in general domain (2.3). In Section 4, we investigate the LD model (1.1) and derive a boundary condition in the limit ε → 0. This boundary condition is then used in Section 5 to connect the LD model to a BD model with Robin boundary condition in the dual limits of high friction, β → ∞, and large potential barrier, ϕ → ∞. Numerical examples supporting our analysis are provided in Sections 4.1 and 5.1. We conclude with discussion in Section 6.
3. Brownian Dynamics. In this section we consider the overdamped problem in which the particle moves by BD, i.e. its position X(t) evolves according to equation (1.2) with the interaction potential given by (2.6).
3.1. Simple three-dimensional example with explicit solution. Before investigating the general problem in R d , as a warm-up we first consider a simplified special case to gain insight into how to recover the Robin boundary condition (2.12) from an interaction potential. We consider the spherically symmetric domain (2.2) for d = 3. The steady-state spherically-symmetric Fokker-Planck equation (2.11) is then given by where the fully adsorbing boundary condition (2.5) implies a Dirichlet boundary condition at r = r b , i.e. p(r b ) = 0. We consider constant concentration p ∞ > 0 far away from the reactive surface, i.e. 1 lim r→∞ p(r) = p ∞ .
We also assume that the potential ϕ(r) is given by (2.6) where ψ : [0, 1] → [0, 1] is the linear function Substituting into (3.1), and making use of the boundary conditions at r = r b and r = ∞, we can solve piecewise to obtain a solution on (r b , r b + ε) and a solution on (r b + ε, ∞). On each interval the solution is defined up to an unknown constant. By enforcing continuity of the flux at r = r b we can eliminate one constant to obtain where the last constant, A, is determined by the continuity of p(r, t) at r = r b + ε. This gives To recover the Robin condition (2.12) as ε → 0 and ϕ → ∞, we need lim ε→0 ϕ→∞ which is equivalent to We note that Therefore, β(ε, ϕ) will have a finite limit if ε exp[ ϕ ]/ϕ has a finite limit. This motivates the choice which converges to zero as ε → 0. We therefore conclude that lim ε→0, ϕ→∞, where ε and ϕ are related by (3.8) β(ε, ϕ) = lim ε→0, ϕ→∞, where ε and ϕ are related by (3.8) , so that we get (3.6). In particular, we recover the Robin boundary condition with Robin constant K.
The steady-state solution to the limiting Robin boundary condition problem (2.12) with (3.2) is given by We now examine the error between this limit and the solution (3.4) for two cases: For the latter we have We can estimate that the denominator is greater than r. Consequently, Using (3.7) and the definition of ε (equation (3.8)), we have β(ε, ϕ) ∼ O(1) as ε → 0. Using (3.9), we conclude that the second term on the right hand side is O(ε/r) as ε → 0. Thus, we conclude On the other hand, we have The maximum error between the two models is therefore O(1) as ε → 0 (where ϕ satisfies (3.8)), illustrating the non-uniformity of the error within the region where the potential interaction is non-zero. Note, while the maximum norm (i.e. L ∞ norm) of the difference between solutions does not converge on (r b , r b + ε), both p(r) and p(r) are uniformly bounded on As such, using that the maximum norm error is O(ε/r) on (r b + ε, ∞) by (3.10), we find that the q-norm converges for any 3 < q < ∞, Here, the lower bound q > 3 is an artifact of our working on an unbounded domain, which requires integrability of |p(r) − p(r)| q r 2 on (r b , ∞). If our domain was bounded, i.e. given by equation (2.4), and we used the Dirichlet boundary condition that p(L) = p ∞ , we would expect a similar estimate to hold for all 1 ≤ q < ∞.
In summary, we have found that in the special case of a linear potential barrier, by taking the height of the barrier ϕ → ∞ and then choosing the width of the barrier, ε, to decrease exponentially in the height (i.e. according to equation (3.8)), we can recover the solution to the diffusion equation with Robin boundary condition. We may therefore interpret the Robin boundary condition model for bimolecular reactions between two molecules as approximating an underlying interaction potential, in which the two molecules must surmount a steep potential barrier before entering a bound state represented by a deep well.
Remark. In [31] it was shown how the Robin constant can be chosen to give the same diffusion-limited reaction rate in two steady-state spherically-symmetric models both including the same fixed interaction potential, with one using a zero Dirichlet boundary condition at r b (as in (3.1)), and the other using a Robin condition at r b + ε for any fixed ε > 0. The argument in [31] can be modified to match diffusion-limited reaction rates in the two models considered in this section (i.e. (3.1) with p(r b ) = 0, and the steady-state diffusion equation with Robin boundary condition satisfied bȳ p(r)). It is satisfying to note that choosing K to match the diffusion limited rates in this manner also gives the relationship (3.8) between K, ε and ϕ. We therefore conclude that choosing K to match the diffusion-limited reaction rates of the two models of this section is sufficient to give convergence of p(r) top(r) as ε → 0 and ϕ → ∞.

3.2.
Asymptotics for a general BD model with a general interaction potential. To analyze equation (2.8) for arbitrary d ∈ N with general potential ϕ(r) in the form (2.6), we follow the approach of Pego [29] and Li [25]. We use the method of matched asymptotic expansions in a general bounded domain Ω which satisfies (2.3) and has a smooth boundary. First, we consider an inner solution in a boundary layer near ∂B r b (0). Let denote the stretched distance of a point x from the reactive boundary. We define an inner solutionp(z, x, t). The second argument should actually bex = x/ |x|, but we follow the method of [29] and instead assume that for all λ > 0 (i.e. the length of x is accounted for by z). We note the identities that where we have used (3.12) and that asp is assumed constant when the second argument is varied along thex direction. From (2.8) we find the inner problem (for 0 < z < 1) We now expand the inner and outer solutions in ε and denote the leading-order terms byp 0 and p 0 respectively. The leading-order behavior of the inner solution then satisfies Note that we will find below that in the distinguished limit ϕ depends on ε in a logarithmic manner (see equation (3.13)) so that retaining ϕ here is consistent with neglecting terms of O(ε). Solving this equation and using the Dirichlet boundary condition at z = 0, we find where A(x, t) is an unknown constant, also satisfying the dilation property The leading order term of the outer solution, p 0 , satisfies the diffusion equation (2.9) away from the reactive boundary. Our matching conditions are Combining these two equations, for where the last approximation is obtained using Laplace's method in the limit ϕ → ∞.
Thus we recover the desired Robin boundary condition (2.10) if ϕ → ∞ and ε → 0 such that In summary, we have used the method of matched asymptotic expansions to examine scaling limits for a general set of potential interactions between two particles. The potentials were assumed to be short-range, forming a high barrier as the separation between the molecules approaches a fixed "reaction-radius", and then a deep well once this barrier is surmounted. In the limit that the height of the barrier, ϕ → ∞, and the width of the barrier, ε, decreases to zero exponentially in the height, we recover the solution to the diffusion equation with Robin boundary condition. We may therefore interpret bimolecular reactions between two molecules modeled with a Robin boundary condition as an approximation to one of many possible underlying potential interactions. These interactions are characterized by the two molecules needing to surmount a steep potential barrier before entering a bound state represented by a deep well.
Remark. If ε is chosen to approach zero slower than (3.13), then we recover a zero Neumann boundary condition at the reactive surface, Likewise, if ε is chosen to approach zero faster than (3.13), then we recover the zero Dirichlet boundary condition that

Langevin Dynamics.
We now consider the LD model (1.1) in a bounded d-dimensional domain satisfying (2.3). The reactive target is again taken to be the surface of the d-dimensional sphere of radius r b about the origin. It is again assumed that the molecule is adsorbed instantly upon reaching the surface of the sphere, i.e. we consider the boundary condition (2.5) together with spherically symmetric interaction potential ϕ given by (2.6). We leave unspecified any boundary condition on ∂Ω \ ∂B r b (0), as it is not needed in the following analysis.
Instead of studying the Langevin equation (1.1) we work with the corresponding Kramers equation (2.7). Since x/r b is the normal to ∂B r b (0) at x, the adsorbing boundary condition (2.5) means that the Kramers equation (2.7) is coupled with the Dirichlet boundary condition and the unspecified boundary condition on ∂Ω \ ∂B r b (0). We are interested in various limits of (2.7) as ε → 0, ϕ → ∞ and β → ∞ in which the interaction potential can be approximated by an appropriate reactive boundary condition. In the BD models of Section 3 our goal was to derive the widely-used Robin boundary condition. To our knowledge, in the LD model we now investigate there is no standard reactive boundary condition for bimolecular reactions. Therefore, we wish to see what, if any, reactive boundary condition arises when considering similar limits of ε and ϕ. As in Section 3.2, to study these limits we will match an inner solution within an O(ε) boundary layer about ∂B r b (0) to an outer solution when far from the reactive boundary. Using the same notation as in Section 3.2, we denote by z the stretched distance from the boundary, given by (3.11). We also introduce a re-scaled radial velocity, wherex = x/ |x| is a unit vector in the x direction, v r = v ·x is the radial velocity into Ω from ∂B r b (0), andv r is a scaling constant in the radial velocity that will be specified later. In the inner region we denote the density byp(z, x, ν, v, t), where we assume thatp is constant whenever x is varied in the radial direction, and when the component of the velocity, v, in the radial direction is varied. That is for α ∈ R.
In addition to the identities (3.12), we note that Using (3.12) and (4.2) we find the derivative operators in (2.7) transform in the new coordinates to where we have used thatx Using (2.6), equation (2.7) can be transformed to We consider an asymptotic expansion of the inner solution,p (in the boundary layer about ∂B r b (0)) as ε → 0 with all other parameters held fixed, Similarly we also consider an expansion of the outer solution, valid far from ∂B r b (0), for which we will abuse notation and also denote it by p, To simplify this equation, we letv r = ϕ β D. This choice emphasizes large velocities, for which we expect the particle to be most likely to escape over the (reactive) potential barrier. Substituting (4.4) into (4.3), we obtain The boundary condition (4.1) implies divides the half-plane into three regions. In the region where ν > 0 and 1 − ψ(z) < ν 2 /2,p (0) is zero due to the boundary condition (4.6). This region corresponds to particles originating at the boundary and moving into the domain. Since particles are only adsorbed at the boundary, and not emitted, we find that the solution is zero throughout the region. Where ν < 0 and 1 − ψ(z) < ν 2 /2 the solution will be determined by matching to the outer solution. This region corresponds to particles entering from outside the boundary layer with sufficient velocity to escape over the potential barrier, and thereby exit through the reactive surface at z = 0. Finally, in the region 1 − ψ(z) > ν 2 /2, trajectories with ν < 0 are reflected in a symmetric manner. For points with ν < 0 the solution will again be given by matching to the outer solution, while for those with ν > 0 the solution is determined by reflection. This region corresponds to particles that enter from outside the boundary layer with insufficient velocity to escape over the potential barrier. These particles are then reflected, moving back out of the boundary layer. Note, in the original variables the curve separating these three regions is given by For x ∈ ∂B r b (0), we match the inner solution as z approaches the edge of the boundary layer to the outer solution as the particle's spatial position, y, approaches x. That is, we require Using that ψ(1) = 0, this matching condition implies the outer solution satisfies the following reactive boundary condition for v r ≥ 0 and x ∈ ∂B r b (0) (4.7) When the moving particle reaches the reactive boundary with radial velocity in the outward direction greater than √ 2 ϕ β D it leaves the domain (i.e. undergoes reaction). In contrast, when the particle reaches the boundary with a slower radial velocity in the outward direction it is reflected back into the domain along the direction of the normal at the point where it hit the reactive boundary. This "specular reflection" boundary condition is also obtained in the corresponding deterministic Newtonian mechanics model. We demonstrate this explicitly for a simple one-dimensional example in Appendix A.
A version of this boundary condition, given in terms of an arbitrary threshold velocity, is assumed in the kinetic boundary layer investigations of the one-dimensional Kramer's equation in [7,8], and the 3D spherically-symmetric steady-state Kramer's equation in [23]. It is also (briefly) mentioned in [23] that the specific threshold of √ 2ϕβD we derive is what one might impose across an interface where the potential is discontinuous with a jump of size ϕ (again, for the 3D spherically-symmetric steadystate Kramer's equation). Our asymptotic analysis shows hows this specific threshold velocity arises in the general Kramer's equation as the limit of a shrinking potential boundary layer bordering a Dirichlet boundary condition. We obtain an effective jump in potential at the reactive boundary, as opposed to an interface within the domain. Contrast this to the limit of the BD problem from Section 3.2, where in taking ε → 0 with ϕ fixed the influence of the potential is completely lost (e.g. a zero Dirichlet boundary condition is recovered).
We therefore conclude that the LD model with the interaction potential is in the limit that ε → 0 equivalent to solving the (zero-potential) Kramers equation with the specular reflection reactive boundary condition (4.7) on ∂B r b (0) and whatever boundary condition was imposed on ∂Ω. That is, in the limit that the width of the potential approaches zero, with the barrier height held fixed, we find that the potential can be approximated by a velocity threshold boundary condition. Here particles moving sufficiently fast relative to the height of the potential barrier undergo bimolecular reactions when reaching the reactive boundary, while those moving too slow are reflected back into the domain. This result should be applicable for general short-range potential interactions that form a high barrier as the separation between two molecules approaches a fixed "reaction-radius", and then a deep well once this barrier is surmounted. In contrast to the diffusive case, taking the barrier height ϕ → ∞ (with β fixed) leads to a complete loss of reaction; all particles reaching the reactive boundary are simply reflected back into the domain. We choose a small time step ∆t and compute the position X(t + ∆t) and the velocity V (t + ∆t) from the position X(t) and the velocity V (t) by where χ [0,ε] : R → {0, 1} is the characteristic function of the interval [0, ε] and ξ is a normally distributed random variable with zero mean and unit variance. We implement adsorbing boundaries at both ends x = 0 and x = L of the simulation domain [0, L] by terminating the computed trajectory whenever X(t) < 0 or X(t) > L.
In this paper, we are interested in understanding the dependence of the behavior of the LD model (1.1) on its parameters ε, ϕ and β. In particular, we choose the values of other parameters equal to 1, namely D = L = 1. (4.11) In this section, we are interested in the limit ε → 0. We want to illustrate the boundary condition (4.7). Therefore we fix the values of ϕ and β and simulate the LD model (1.1) for different values of ε.
For each value of ε, we compute many trajectories according to (4.9)-(4.10), starting from the middle of the domain, i.e. X(0) = L/2, with initial velocity V (0) sampled from the normal distribution with zero mean and variance β D. Whenever a particle enters the region [0, ε], we record its incoming velocity. Then we follow its trajectory in the region [0, ε] and record one of two possible outputs: (1) the particle leaves Ω through its left boundary (i.e. X(t) < 0); or (2) the particle returns back to the region (ε, L] of domain Ω (i.e. X(t) > ε). The fraction of particles which left the domain Ω through its left boundary as a function of the incoming velocity is plotted in Figure 4.1(b). This curve can be interpreted as the probability that the particle escapes over the potential barrier, given its incoming velocity. We also plot the vertical dashed line at − √ 2 ϕ β D in Figure 4.1(b). This threshold separates the two cases of boundary condition (4.7). We observe that the probability of escape converges to the step function as ε → 0, i.e. we have numerically confirmed boundary condition (4.7).
We will return to this example in Section 5.1 when we study the convergence of the LD model to the diffusion process with a Robin boundary condition. In particular, we use values β = 10 3 and ϕ = 2.593 in Figure 4.1(b). This choice of ϕ will be explained in the following section and then used again in one of our simulations presented in Section 5.1.

From
Langevin Dynamics to Brownian Dynamics. We now study the overdamped limit of the LD model (1.1). We wish to show that in the overdamped limit where β → ∞, taking ϕ → ∞ in a β-dependent manner will recover a Robin boundary condition for the limiting diffusion equation. To study the β → ∞ limit, we extend the asymptotic analysis of the one-dimensional Kramers equation in [14] to the d-dimensional (zero-potential) Kramers equation (4.8). The kinetic boundary layer studies in [7,8,23] previously investigated the relationship between the velocity threshold in the specular reflection boundary condition and effective adsorption rate, K, in the Robin condition. Our approach here differs from those studies, which primarily used numerical solutions of truncated moment equations or basis function expansions to estimate empirically determined formulas for the Robin constant, K. We focus on deriving an explicit formula relating how the potential barrier height, ϕ, should be chosen as β → ∞ to recover a specified Robin constant.
We begin by re-scaling velocity as v = √ β η and let f (x, η, t) = p(x, v, t). Substituting into (4.8), we obtain Substituting into (5.1), we find Implicit in these equations is the assumption that we are interested in timescales for which We therefore interpret (5.2) as implying that the velocity distribution component of f 0 relaxes to equilibrium on a faster timescale than t. Denote by τ this faster timescale. As discussed in the introduction to [9], up to a normalization constant there is a unique solution to (5.2) that corresponds to the equilibrium solution of the fast-timescale, time-dependent equation This equilibrium solution is then where (x, t) is independent of η. We similarly find that the general solution to (5.3) is given by where ξ(x, t) is also independent of η. Substituting these into (5.4) we see that Integrating both sides of this equation for all η ∈ R d , and using that f 2 (x, η, t) → 0 and |∇ η f 2 (x, η, t)| → 0 as |η| → ∞, we conclude that (x, t) satisfies the diffusion equation The probability density that the particle has position x at time t and has not reacted is given by implying that to leading order u is proportional to . As such, we expect as β → ∞ that u satisfies the diffusion equation. We now show that , and hence u, satisfies a Robin boundary condition as β → ∞ when ϕ is chosen to approach infinity in a β-dependent manner. To leading order in β, the outward flux through a point, where as in previous sectionsx = x/ |x|. Let denote the set of velocities at which particles may escape through the reactive boundary. The specular reflection boundary condition (4.7) implies Comparing with (5.5), we find that to leading order , and hence u, satisfy the reactive Robin boundary condition for β → ∞ As ϕ → ∞, the denominator of K approaches one. This suggests the scaling We then find (x, t) satisfies the desired Robin boundary condition, in the limit β → ∞. More generally, we could impose the equality √ β K(ϕ, D) = K in equation (5.6). Solving for β, we obtain the following relation between ϕ and β: (5.8) In Figure 5.1(a), we compare both formulae (5.7) and (5.8) for K = D = 1. As expected, they are equivalent in the limit β → ∞. In the following section, we will present illustrative simulations, confirming that (5.8) provides a more accurate approximation of the limiting Robin boundary condition. Further improvements to these formulas could presumably be made by incorporating a more detailed representation of the kinetic boundary layer near ∂B r b (0).

A numerical example illustrating limit β → ∞.
As in Section 4.1, we consider the LD model (1.1) for d = 1 in the interval Ω = [0, L] with the linear potential (3.3) near the left-hand boundary only. We choose a small time step ∆t and compute the position X(t + ∆t) and the velocity V (t+∆t) from the position X(t) and the velocity V (t) by (4.9)-(4.10). We implement adsorbing boundaries at both ends x = 0 and x = L of the simulation domain [0, L] by terminating the computed trajectory whenever X(t) < 0 or X(t) > L.
In this section, we want to show that both equations (5.7) and (5.8) correctly recover the Robin boundary condition in the limit β → ∞. In particular, we choose the values of other parameters equal to 1, namely (compare with (4.11)) K = D = L = 1. (5.9) We vary β and we use either equation (5.7) or equation (5.8) to calculate the corresponding value of ϕ (these values are plotted in Figure 5.1(a)). For each value of β, we compute 10 5 trajectories according to (4.9)-(4.10), starting from the middle of the domain, i.e. X(0) = L/2, with initial velocity V (0) sampled from the normal distribution with zero mean and variance β D. Each trajectory is calculated until it leaves the domain [0, L] either through the left or right boundary point. In Figure 5.1(b), we plot the probability that a trajectory leaves the domain Ω through the left boundary (the so-called splitting probability), estimated as the fraction of all trajectories which are terminated because X(t) < 0.
Our goal is to illustrate that equations (5.7) and (5.8) can be used to connect the LD model with the limiting Robin boundary problem (2.9)-(2.10). Let Π(x) be the probability that the BD particle leaves the domain Ω = [0, L] through the left boundary. Since and Π(L) = 0, we find Since all trajectories start from the middle of the domain, X(0) = L/2, we have for the parameter values given by (5.9). This value is plotted in Figure 5.1(b) as the black dashed line. We confirm that the results estimated from simulations approach (5.10) as β → ∞. We also confirm that simulations based on the higher-order approximation (5.8) converge more quickly to the limiting Robin boundary problem than the simulations based on equation (5.7).
6. Discussion. We have considered three parameters, ε (potential width), ϕ (potential height) and β (friction constant), and studied several limits of these parameters which lead to the Robin (reactive) boundary condition (2.10). Parameters ε and ϕ are shared by both the BD and LD models. In Section 3, we have shown that the BD model can recover the Robin boundary condition in the limit ε → 0 and ϕ → ∞ when these parameters are related by (3.13). For the case of the linear potential (3.3) this relation can be rewritten as The LD model has an additional parameter β. In Section 5, we have derived two formulae (5.7) and ( Considering that (experimentally determinable) parameters K and D are given constants, we can compare our BD result (6.1) with our LD result (6.2). They can be both used to specify the height of the potential barrier ϕ, which is given as a function of ε in (6.1) and as a function of β in (6.2). On the face of it, the LD result (6.2) does not depend on the parameter ε. However, the derivation of the LD result (6.2) is only valid for small ε. More precisely, our conclusion for the LD model (1.1) can be stated as follows: (1) If ε D β 1, then ϕ is independent of ε and can be written in the form (6.2); (2) If 1 ε D β , then ϕ is independent of β and is given by (6.1). In our illustrative computations in Figure 5.1(b) we have used ε = 10 −3 . In this case, the BD result (6.1) implies that ϕ = 9.118. We observe that this value is higher than the values of ϕ plotted in Figure 5.1(a) which are used for our simulations in Figure 5.1(b). We can also substitute this value ϕ = 9.118 into (6.2). We obtain β = 5.224 × 10 8 . For these values of β and ε we are in case (2), so that the LD result (6.2) would no longer be applicable. We instead recover the limiting Robin boundary condition (in the limit β → ∞ with ε = 10 −3 ) by using the BD result (6.1).
The above results (6.1)-(6.2) can be used to connect the experimentally determinable parameters K and D with parameters of computer simulations to design reaction-diffusion models based on LD, i.e. to simulate diffusion-limited bimolecular reactions between Kramers particles. There are also other situations where our analysis will be applicable. One of them is adsorption to surfaces [15,17]. Our set up includes interactions with a reactive surface of a sphere and can be used in modeling interactions of small molecules with large reactive spheres, for example, for adsorption of polymers to a surface of a virus [17] or for coating of spherical particles by reactive polymers [33]. We also expect our results should be easy to extend to general non-spherical reaction surfaces, assuming sufficient regularity. Another possible application area is modeling excluded-volume effects. We have observed that the short range repulsive interaction potential (2.6) leads to zero Neumann boundary condition (3.14) if ε is chosen to approach zero slower than (3.13) (in the limit ϕ → ∞). A similar potential mechanism has been used to enforce Neumann boundary conditions on global domain boundaries in [5], and can be used to model excluded-volume effects in models of intracellular macro-molecular crowding [4,6].
Since LD requires a smaller time step than the overdamped BD model, the numerical simulations are in general more computationally intensive for LD. However, for many biological applications the LD model is only required close to a reactive surface (where we have a non-zero potential). In particular, one could replace the computationally intensive LD model with the BD model in the part of the computational domain which is far from the reactive surface [12]. In some applications, one could further substitute the BD simulation algorithms by even coarser and more efficient simulation techniques, including lattice-based simulations [11,19] or even mean-field equations [15,20]. In this way, one could design LD simulation methods which simulate intracellular processes on comparable time scales as the BD simulation packages which are available in the literature [3,30,34]. In some applications, one could also design an efficient first-passage-time scheme by replacing discretization (4.9)-(4.10) (which uses fixed time step ∆t) by estimating the next time when a Kramers particle becomes sufficiently close to the reactive target [21]. Similar approaches have been used to accelerate BD simulations in the literature [27,28,35].
Appendix A. A simple example illustrating how specular reflection arises in a classical mechanical system. In this section we illustrate with a simple example how the specular reflection boundary condition (4.7) arises in a classical mechanical system in the absence of noise.
We consider the simple case of a particle moving in the infinite one-dimensional potential, x > 0, and experiencing friction, with friction constant β. Note, for the purposes of this section we only make use of the value of the potential for x ∈ (−∞, 0]. Newton's equation when the particle's position, X(t), satisfies X(t) ≤ 0 with negative velocity is then Using the Einstein Relation this reduces to the one-dimensional analogue of the LD model (1.1) with the noise term neglected. We consider the initial conditions and ask for what initial velocity range the molecule successfully moves a distance ε to the left before changing direction and falling down the potential gradient. This is consistent with a molecule reaching the reactive boundary in the Kramer's equation model (2.7). Given the molecule's negative initial velocity, until the molecule's direction of motion reverses and the molecule moves back across the origin we find The velocity then becomes zero at time At this time, the molecule's position is The molecule then successfully travels further than −ε, analogous to penetrating the "reactive boundary" at x = −ε, if or equivalently that Dϕ εβ (α − ln(1 − α)) < −ε.
As ε → 0 We therefore find that the molecule successfully penetrates the reactive boundary in the limit that ε → 0 if and only if v 0 < − 2ϕβD, consistent with the specular reflection boundary condition (4.7).