Numerical calculation of protein-ligand binding rates through solution of the Smoluchowski equation using smoothed particle hydrodynamics

Abstract

Background

The calculation of diffusion-controlled ligand binding rates is important for understanding enzyme mechanisms as well as designing enzyme inhibitors.

Methods

We demonstrate the accuracy and effectiveness of a Lagrangian particle-based method, smoothed particle hydrodynamics (SPH), to study diffusion in biomolecular systems by numerically solving the time-dependent Smoluchowski equation for continuum diffusion. Unlike previous studies, a reactive Robin boundary condition (BC), rather than the absolute absorbing (Dirichlet) BC, is considered on the reactive boundaries. This new BC treatment allows for the analysis of enzymes with “imperfect” reaction rates.

Results

The numerical method is first verified in simple systems and then applied to the calculation of ligand binding to a mouse acetylcholinesterase (mAChE) monomer. Rates for inhibitor binding to mAChE are calculated at various ionic strengths and compared with experiment and other numerical methods. We find that imposition of the Robin BC improves agreement between calculated and experimental reaction rates.

Conclusions

Although this initial application focuses on a single monomer system, our new method provides a framework to explore broader applications of SPH in larger-scale biomolecular complexes by taking advantage of its Lagrangian particle-based nature.

Background

In the “perfect” enzyme [1] acetylcholinesterase (AChE), the rate-limiting step for catalysis is diffusional encounter [2,3]. Specifically, the active site lies at the bottom of a 20 Å-deep gorge, and the diffusion of substrate into it is accelerated by electrostatic steering [4,5]. Its diffusion-limited behavior, complex geometry, and strong electrostatic influence has made AChE a useful target for both experimental and computational studies of biomolecular diffusion [4-11].

Two major classes of methods have been used to estimate diffusion rates in biomolecular systems. Mesoscopic coarse-grained methods like Monte Carlo [12-14], Brownian dynamics (BD) [8,9,15], and Langevin dynamics [16,17] simulations trace the trajectories of individual coarse-grained particles driven by Brownian motion. Such simulations typically consider dilute ligand concentrations so that electrostatic protein-ligand interactions can be modeled by the Poisson-Boltzmann equation [18,19] with a few notable exceptions [20]. Alternatively, continuum models can be used to treat the diffusion of ligand concentration in space around a biomolecule by the Smoluchowski equation [6,21-25]. In particular, an adaptive finite element approach [26] has been used to numerically solve the Smoluchowski equation, and it shows higher accuracy in predicting experimental data about the ligand binding rates than the coarse-grained BD modeling [6]. For dilute ligand concentrations, electrostatic interactions can also be modeled with the Poisson-Boltzmann equation like the mesoscale approach [6,7]. However, for more concentrated ligand solutions, continuum models can also model the electrostatic potential near the biomolecular surface using a regularized Poisson-Nernst-Planck formulation [24,25], allowing screening of the ligand-receptor interactions by its time-dependent distribution around the protein.

Here, we follow the continuum approach but solve the Smoluchowski equation using a new smoothed particle hydrodynamics (SPH) method [27,28]. Unlike Eulerian grid-based methods such as finite element method (FEM), SPH is a Lagrangian particle-based method. SPH has been used with good accuracy for numerically solving partial differential equations (PDEs) describing momentum, mass and energy conservation laws [27]. In SPH, the domain is discretized into a set of “particles” that serve as interpolation points to numerically solve the governing PDEs. The SPH discretization of PDEs is based on a meshless interpolation scheme, which allows the PDEs to be written in the form of a system of ordinary differential equations (ODEs). SPH has a straightforward discretization without the need for time-consuming FEM mesh construction around complicated geometries such as biomolecules. Due to its Lagrangian nature, SPH has many advantages for modeling physical phenomena involving moving boundaries, large deformation of materials, multiphases, and advection-dominated diffusive transport [28-30]. Specifically, in SPH, free surfaces and interfaces between fluids move with particles, and hence, there is no need for front tracking schemes. And the non-linear advection term is embedded in the material derivative in the Lagrangian coordinate system, and hence, SPH models advection exactly. In addition, the similarity of SPH to molecular dynamics and mesoscopic coarse-grained particle methods (e.g., dissipative particle dynamics, BD, and Langevin dynamics), allows coupling of simulations across scales to build a multiscale modeling framework. This is our primary goal with the current work: to enable the multiscale and multiphysics description of biomolecular dynamics and ligand recognition. To the best of our knowledge, SPH has not been widely used in modeling biomolecular systems. Thus, in the present work, we aim to take the first step to introduce SPH into this field through the development of a SPH model for biomolecular diffusion with AChE as a test case.

In the SPH model, the Smoluchowski equation is numerically solved and the ligand binding rates are calculated from flux across the reactive boundary as in the previous studies using FEM [6,21-25]. However, in the previous FEM studies, active sites were modeled using the absolute absorbing (Dirichlet) boundary condition (BC). This BC has a simple description on the reactive boundaries but assumes infinitely fast chemical reactions between the enzyme and the ligand; i.e., a “perfect enzyme”. In our model, we take into account imperfect and non-instantaneous reactivity and thus solve the equation subject to a reactive (Robin) BC.

To solve the Smoluchowski equation subject to Robin BC using SPH, we use a continuum surface reaction method [31] which we have recently adapted to solve the Navier-Stokes equations subject to slip (Robin) boundary conditions [32]. In this formulation, the Robin BC is replaced by a reflective Neumann BC and a source term added into the governing equation. The derivation of the method is based on the approximation of the sharp boundary with a diffuse interface of finite thickness by means of a color function. This method is general for any arbitrary complex geometries and thus appropriate for modeling Robin BC in biomolecular systems with complex structures.

Results and discussion

Spherical test systems

Before the numerical method was applied to a biomolecular system with complicated geometry, we verified it on simple spherical test cases. Specifically, we considered a diffusing sphere with a radius R ₁. The entire domain was confined by the outer boundary Γ _b determined as a spherical surface with the radius of R ₂=125 Å.

For the first test case, we let R ₁=0 and assume no external potential, for which the time-dependent analytical solution of the Smoluchowski equation can be easily derived. Figure 1 compares the SPH numerical solutions with the analytical solution at different times. SPH solutions are compared at different resolutions and their corresponding L ₂ errors are calculated relative to the analytical solution. Figure 1 shows that, even at the coarsest resolution (Δ x=8 Å), the SPH solution agrees well with the analytical solution with about 3% relative error. This relative error is further reduced to 1% by increasing the resolution to Δ x=2 Å.

Figure 1

Comparison of SPH solutions to the analytical solutions for the Smoluchowski equation subject to the Dirichlet BC on r=R ₂ at different times with the relative L ₂ errors for different resolutions. Specifically, L ₂=0.0326 for Δ x=8Å(green square), L ₂=0.0180 for Δ x=4Å(blue circle), and L ₂=0.0103 for Δ x=2Å(red triangle).

Next, the spherical system is assumed to have a Coulombic form of the PMF, i.e., W(r)=q/β r with +1 e charge. We set R ₁=50 Å and impose either a Dirichlet BC as specified in Eq. 6 or a Robin BC as in Eq. 5. In these two tests, the corresponding SPH solutions of concentration at steady-state are compared with the analytical solutions. The converged SPH solutions are shown for the Dirichlet BC (Figure 2) and Robin BC (Figure 3) imposed on the inner spherical boundary (r=R ₁). The reactive coefficient for the Robin BC is α=1×10³. In both tests, the SPH solutions show very good agreement with the analytical solution even at the resolution of Δ x=8 Å, which can be further improved with increasing resolution to Δ x=2 Å. Moreover, at Δ x=2Å, the calculated reaction rate is 2.83×10¹² M ⁻¹ m i n ⁻¹ for the Dirichlet BC, and is 8.24×10¹¹ M ⁻¹ m i n ⁻¹ for the Robin BC, both with L ₂ errors less than 3% relative to the analytically evaluated ones.

Figure 2

Comparison of SPH solutions to the analytical solution for the Smoluchowski equation subject to the Dirichlet BC on both r=R ₁ and r=R ₂ at steady-state with the relative L ₂ errors for different resolutions. Specifically, L ₂=0.0666 for Δ x=8Å(green square), L ₂=0.0321 for Δ x=4Å(blue circle), and L ₂=0.0153 for Δ x=2Å(red triangle).

Figure 3

Comparison of SPH solutions (symbol) to the analytical solution (line) for the Smoluchowski equation subject to the Robin BC on r=R ₁ and Dirichlet BC on r=R ₂ at steady-state with the relative L ₂ errors for different resolutions. Specifically, L ₂=0.00914 for Δ x=8Å(green square), L ₂=0.00598 for Δ x=4Å(blue circle), and L ₂=0.00377 for Δ x=2Å(red triangle).

Application to AChE-ligand binding rates

We applied the SPH method to study the ligand binding kinetics of a simple spherical cationic ligand to the mouse acetylcholinesterase (mAChE) under various ionic strength conditions. Specifically, we performed the time-dependent calculations at ionic strengths of 0.0, 0.05, 0.10, 0.15, 0.20, 0.50 and 0.67 M until the diffusion reaches the steady-state. To achieve the highest accuracy with affordable computational cost, a resolution of Δ x=2 Å was used in all the following calculations.

In previous studies by Song et al. [6], a simple but realistic set of boundaries was used inspired by Tara et al. [9], encompassing the active site as well as the gorge and the peripheral anionic site (PAS) of mAChE. We constructed these spherical active boundaries (Γ _a ) at varying distances from the active site along an axis defined by the carbonyl carbon of S203 at the origin and the gorge. Spheres 1-6 were placed at 16.6, 13.6, 10.6, 7.6, 4.6, and 1.6 Å along the this axis, respectively. The outermost spheres 1 and 2 were assigned radii of 12 and 9 Å, respectively, while all others were given radii of 6 Å. Each reactive boundary N is defined as the intersection of the (surface) union of spheres N through 6 with the mAChE structure. Figure 4A shows the discretized domain with R ₂=128 Å. Figure 4B and 4C depict the constructed reactive boundaries 1 and 4.

Figure 4

Panel A shows the discretized domain with R ₂=128 Å and the mAChE molecule in the center with the reactive boundary shown in purple. Light blue indicates the outer boundary (R ₂), blue the solvent, green the protein, and magenta the first (outermost) reactive boundary. Panels B and C show reactive boundaries 1 and 4, respectively in magenta spheres.

In most prior studies [6,7,23], an absolute absorbing (Dirichlet) BC (Eq. 6) was assumed. However, in the present work, we demonstrated improved performance with the reactive (Robin) BC (Eq. 5) imposed on the reactive boundaries. Figure 5 shows the steady-state spatial distribution of ligand throughout the simulation domain at different ionic strengths. At zero ionic strength, there are three large ligand-attracting regions, two on either side of the active site and one on the opposite side of the protein. There is also one ligand-depleted region at the top and another one near the opening of the gorge. At non-zero ionic strengths, electrostatic screening reduces the size of the ligand-enriched and ligand-depleted regions. However, a large region around the active site remains ligand-depleted at up to 0.50 M ionic strength. Figure 6 illustrates the temporal evolution of the concentration distribution as ligand moves inward in the bulk region from the outer boundary (Γ _b ). The distribution has clearly reached steady state by 190 ns.

Figure 5

Contour of concentration distribution around mAChE (shown in dark gray) with the Robin BC (α=8×10³) on reactive boundary 1 at steady state with a range of ionic strengths. Reactive boundary 1 is shown in purple.

Figure 6

Time evolution of the concentration distribution around mAChE (shown in gray) with a Robin BC (α=8×10³) on reactive boundary 1 at 0.15 M ionic strength. Reactive boundary 1 is shown in purple.

We calculated the reaction rates from these solutions according to Eq. 22. In Figure 7, the left panel shows the time evolution of reaction rate k _on(t) on reactive boundary 1 at different ionic strengths. For this boundary, k _on(t) converges within 150 ns for all ionic strengths. The right panel shows k _on(t) on reactive boundaries 1-4, respectively, at 0.15 M ionic strength.

Figure 7

(Left) k _on as a function of t on reactive boundary 1 at different ionic strengths. Black square: 0.05M; red right-pointing triangle: 0.10M; blue asterisk: 0.15M; green circle: 0.20M; magenta diamond: 0.50M; cyan triangle: 0.67M. (Right) k _on as a function of t on reactive boundaries 1-4, respectively, at 0.15 M ionic strength. Black square: reactive boundary 1; red circle: reactive boundary 2; blue diamond: reactive boundary 3; green triangle: reactive boundary 4.

We have quantitatively compared the reaction rates calculated by SPH with experimental results [4] and previous computational studies by FEM [6,8]. Radic et al. [4] fit their experimentally measured reaction rates as a function of ionic strength using the Debye-Hückel limiting law

$$ k_{\text{on}}= \left(k_{\text{on}}^{0}-k_{\text{on}}^{\mathrm{H}}\right)10^{-1.18|Z_{E} Z_{I}|\sqrt{I}}+k_{\text{on}}^{\mathrm{H}}, $$

((1))

where I is the ionic strength, $k_{\text {on}}^{0}$ is the effective reaction rate at zero ionic strength rate, $k_{\text {on}}^{\mathrm {H}}$ is the effective limiting reaction rate at infinite ionic strength and set to the value of k _on calculated at 0.67 M ionic strength, z _E is the effective enzyme charge, and z _I is the effective inhibitor charge with a fixed value of +1 e. In SPH calculations with the Robin BC, the reaction coefficient α was varied, as shown in Figure 8, to identify the value of 8.0×10³ which optimized agreement between computational and experimental results.

Figure 8

Root mean square deviation (RMSD) of computed by SPH to experimental reaction rates (over 0-0.67 M ionic strengths) vs. α for the Robin BC.

Figure 9 and Table 1 compare the reaction rates from SPH, FEM [6,8], BD, and experimental data by Radic et al. [4]. As noted by Song et al. [6], BD simulations systematically overestimate the experimental k _on, while the FEM produces good agreement with experimental k _on at RMSD = 0.37 M ⁻¹ m i n ⁻¹. With the Dirichlet BC, SPH predicts k _on with RMSD of 0.57 M ⁻¹ m i n ⁻¹, intermediate between FEM and BD results. However, with the Robin BC, SPH predicts k _on with RMSD of 0.33 M ⁻¹ m i n ⁻¹, better than the FEM and BD results.

Figure 9

Reaction rates of mAChE on reactive boundary 1 obtained from different methods. Black: from experimental data [4] (symbol) and fitted (line) to the Debye-Hückel limiting law (Eq. 1); blue: from BD [9]; red: from FEM with Dirichlet BC [6]; green: from SPH with Dirichlet BC; magenta: from SPH with Robin BC using α=8×10³. For standardization, both computed and experimental data are fitted to the Debye-Hückel limiting law.

Table 1 Comparison of Debye-Hückel fits vs. ionic strength between experiment and simulations

We also assessed the accuracy of SPH method for describing the ligand-binding kinetics of a mAChE surface mutant. We tested the surface hexa-mutant (E84Q, E91Q, D280V, D283N, E292Q, and D372N) from Radic et al. [4], which reduces the reaction rate by about a factor of 4 across the 0 to 0.67 M ionic strengths. For the mutants, which are nearly isosteric with the wild-type protein, we used the same SPH model as the wild type, but recalculated the electrostatic potentials for the mutant charge distribution. As presented in Table 1, the Robin BC SPH model has qualitative accuracy: predicting $k_{\text {on}}^{0}$ of 2.23 compared to 1.80 from Radic et al. [4] with a k _on RMSD of 0.88 over the entire ionic strength range studied.

Conclusions

The Robin BC offers a new way to incorporate reactive surfaces into continuum diffusion models for rate calculations. This Robin-based model incorporates a new parameter α, which has units of Å/ μs and can be related to the probability of reaction within distance Δ x to the boundary and time interval Δ t by $P=1-\exp (-\alpha \frac {\Delta t}{\Delta x})$ [33]. Thus, α=0 corresponds to zero reactivity (reflective Neumann BC) while α=∞ corresponds to absolute reactivity (absorbing Dirichlet BC).

There are two possible origins for the differences between the current SPH model results and past FEM calculations using the Dirichlet BC. First, the current SPH work uses a more recent mAChE structure (4B82) while the previous FEM calculations used an older structure (1MAH). Second, our SPH model uses a fixed resolution uniformly on both solution domain and boundaries, while the FEM adaptively meshes the reactive boundary with higher resolution.

This work has provided an initial demonstration that the Lagrangian (particle-based) SPH method out-performs the Eulerian (grid-based) FEM [6] in accurately predicting ligand binding rates in AChE. This result is important because while both methods can be used to study molecules of the size of AChE, SPH is more scalable to larger systems such as the synapse geometry where AChE operates. Additionally, due to its Lagrangian nature, SPH can easily incorporate other physical phenomena such as fluid flow or protein flexibility.

We have demonstrated that superior performance can be achieved using a probabilistic reactive (Robin) BC rather than a simple Dirichlet BC. In fact, the Robin BC is likely more biologically relevant than the Dirichlet BC. While the AChE enzyme is considered nearly “perfect” with a diffusion-limited reaction rate, there is experimental evidence that a very small fraction of substrates entering the active site gorge do not react. Specifically, recent kinetic experiments suggest that through unknown mechanisms, the PAS limits the rate of progression of non-substrates of any size to the catalytic site [34]. In addition, molecular dynamics simulations suggest that the PAS provides a selective gating function, for example by fluctuations in the gorge width that are likely to let acetylcholine but not let larger molecule pass through [35,36].

Methods

Governing equation and boundary conditions

The time-dependent Smoluchowski equation can be written as:

$$ \frac{d p(\mathbf{x},t)}{dt}=\nabla\cdot \mathbf{J} (\mathbf{x},t), ~~~~\mathbf{x}\in\Omega, $$

((2))

where p(x,t) is the concentration distribution of the reactants, and the concentration flux J(x,t) is defined as:

$$ \mathbf{J} (\mathbf{x},t) = D (\mathbf{x})[\nabla p(\mathbf{x},t) + \beta p(\mathbf{x},t) \nabla W(\mathbf{x})], $$

((3))

where D(x) is the diffusion coefficient; for simplicity, it is assumed to be constant. β=1/k _B T is the inverse Boltzmann energy with the Boltzmann constant k _B and kinetic temperature T. W(x) is the potential mean force (PMF) for the diffusing particle due to solvent-mediated interactions with the target molecule. The equation is solved in a three-dimensional domain Ω, subject to the following boundary conditions. First,

$$ p(\mathbf{x},t)= p_{\text{bulk}} ~\text{for} ~\mathbf{x}\in\Gamma_{b}, $$

((4))

specifying a Dirichlet BC on the outer boundary Γ _b where the concentration is equal to a bulk concentration p _bulk. The outer boundary is often a spherical surface with a radius chosen to ensure that the ligand-protein potential is spherically symmetric and/or can be approximated analytically [6]. For the current study with mAChE, this outer boundary has radius R ₂≈128 as determined following a procedure similar to Song et al. and Chen et al. [6,23]. Also following Song et al. and Chen et al., p is normalized such that p _bulk=1.

The active site boundary Γ _a was modeled using either reactive Robin or absolute absorbing Dirichlet BC:

$$ \mathbf{n}(\mathbf{x}) \cdot \mathbf{J} (\mathbf{x},t) = \alpha p(\mathbf{x},t) ~ \text{for} ~ \mathbf{x}\in\Gamma_{a}, $$

((5))

or

$$ p(\mathbf{x})= 0 ~ \text{for} ~ \mathbf{x}\in\Gamma_{a}, $$

((6))

respectively. The coefficient α is chosen to model an intrinsic reaction rate for the active site. Finally, a reflective Neumann BC is defined on the non-reactive boundary of molecule

$$ \mathbf{n}(\mathbf{x}) \cdot \mathbf{J} (\mathbf{x},t) = 0 ~ \text{for} ~ \mathbf{x}\in\Gamma_{m}. $$

((7))

Figure 10 shows the simulation domain along with all boundaries.

Figure 10

Illustration of the simulation domain and all boundaries: Γ _b indicates the outer boundary, Γ _m the molecular surface, and Γ _a the reactive boundary 1; Ω indicates the problem domain between Γ _b and Γ _a ∪Γ _m .

Given a solution to Eq. 2, the reaction rate is calculated from the integral of the flux across the reactive boundary [37]:

$$ k_{\text{on}}= p_{\text{bulk}}^{-1} \iint\limits_{\Gamma_{a}} \mathbf{n}(\mathbf{x}_{s}) \cdot \mathbf{J} (\mathbf{x}_{s},t) d\mathbf{x}^{\prime}_{s}. $$

((8))

In order to solve the Smoluchowski equation (Eq. 2) subject to the reactive Robin BC (Eq. 5), the simulation domain is extended to include a sub-domain Ω _a that is separated from Ω by Γ _a , and we then reformulate Eq. 2 as:

$$ {\fontsize{8.1pt}{9.6pt}\selectfont{ \begin{aligned} {} \frac{dp^{r}(\mathbf{x},t)}{dt}= &\nabla\cdot \left(D (\mathbf{x})[\nabla p^{r}(\mathbf{x},t) + \beta p^{r}(\mathbf{x},t) \nabla W(\mathbf{x})]\right) \\ - &\alpha p^{r} (\mathbf{x}, t) \iiint\limits_{\Omega_{a}} [\mathbf{n}(\mathbf{x}) + \mathbf{n}(\mathbf{x}')]\cdot \nabla_{\mathbf{x}} w(\mathbf{x}-\mathbf{x}', h_{r}) d\mathbf{x}', ~~~\mathbf{x}\in\Omega, \end{aligned}}} $$

((9))

subject to the reflective Neumann BC:

$$ \mathbf{n}(\mathbf{x}) \cdot \mathbf{J}^{r} (\mathbf{x},t) = 0 ~ \text{for} ~ \mathbf{x}\in\Gamma_{a}. $$

((10))

The derivation of Eq. 9 is detailed in Appendix A, which demonstrates

$$ {\lim}_{h_{r} \to 0} p^{r}(\mathbf{x},t) = p(\mathbf{x},t). $$

((11))

In Eq. 9, the normalized kernel function, w(x), is a positive bell-shaped function with at least first continuous derivative and compact support κ h _r such that w(|r|>κ h _r )=0. The value of κ depends on the specific functional form of w(x), which is specified in Section ‘SPH discretization of equations and boundary conditions’. In particular, w(x) satisfies the following conditions:

$$ \iiint\limits_{\Omega\cup\Omega_{a}} w(\mathbf{x}-\mathbf{x}', h_{r}) d\mathbf{x}' =1 $$

((12))

and

$$ {\lim}_{h_{r}\to 0}w(\mathbf{x}-\mathbf{x}', h_{r})=\delta(\mathbf{x}-\mathbf{x}'). $$

((13))

The normal unit vector n in Eq. 9 can be found in terms of a smoothed color function $\tilde {\phi }$ as defined in Appendix A:

$$ \mathbf{n}(\mathbf{x}) = \frac{\nabla \tilde{\phi}(\mathbf{x})}{|\nabla \tilde{\phi}(\mathbf{x})|}, ~~~ \mathbf{x} \in \Omega\cup\Omega_{a}. $$

((14))

SPH discretization of equations and boundary conditions

In this section, we present SPH discretization of the Smoluchowski equation, using Eq. 2 if the Dirichlet BC is used and Eq. 9 if the Robin BC is assumed. To simplify notation, we omit superscript r for the variables in Eq. 9 in the subsequent derivations.

The domain Ω, and the boundaries Γ _a and Γ _b (extended as domains Ω _a and Ω _b respectively), are discretized with a set of N points with positions denoted by a vector r _i (i=1,...,N). The points (which are commonly referred to as particles in SPH) are used to discretize and solve the governing equation. Initially, the particles are distributed uniformly (e.g., placed on a regular cubic lattice) with d _i as the prescribed number density at r _i . The discretization is based on a meshless interpolation scheme:

$$ A_{i} \approx \sum\limits_{j} \frac{A_{j}}{d_{j}}w(\mathbf{r}_{ij},h), $$

((15))

where, A _i =A(r _i ) is a function defined at particle i, A _j =A(r _j ) is the function defined at particle i’s neighboring particles j with distances r _ij =r _i −r _j , and w(r _ij ,h) is the weighting kernel function. The interpolation scheme assumes a summation over all neighboring SPH particles but, due to the compact support of w, only particles within distance κ h from r _i have a non-zero contribution to the summation. Spatial derivatives of A can be calculated as

$$ \nabla_{i} A_{i} \approx \sum\limits_{j} \frac{A_{j}}{d_{j}} \nabla_{i} w(\mathbf{r}_{ij},h). $$

((16))

In the present work, we use a cubic spline kernel as the weighting function

$$ w(\mathbf{r},h) = \frac{1}{\pi h^{3}} \begin{cases} 1-\frac{3}{2}q^{2}+\frac{3}{4}q^{3} & ~~~0 \leq q \leq 1 \\ \frac{1}{4}(2-q)^{3} & ~~~1 < q \leq 2 \\ 0 & ~ q>2, \end{cases} $$

((17))

where q=|r|/h. With this form of weighting function, only particles within 2h distance from particle i contribute to the summations in the SPH equations. We have chosen h=1.3Δ x where Δ x is the size of the cubic lattice.

The SPH approximation of functions and their spatial derivatives allows the Smoluchowski equation subject to the Dirichlet BC (Eq. 2) to be written as a ODE governing the evolution of concentration on particle i as:

$$\begin{array}{@{}rcl@{}} \frac{d p_{i}}{dt} &=& \sum_{j\in \Omega\cup\Omega_{b}\cup\Omega_{a}} \frac{D_{i}+D_{j}}{d_{j}} (p_{i}-p_{j}) \frac{1}{r_{ij}} \frac{d w(r_{ij}, h)}{d r_{ij}} \\ &&+ \beta \sum_{j\in \Omega} \frac{D_{i} p_{i}\,+\,D_{j} p_{j}}{d_{j}} (W_{i}\,-\,W_{j}) \frac{1}{r_{ij}} \frac{d w(r_{ij}, h)}{d r_{ij}}. \end{array} $$

((18))

The derivations of the first and second terms on the right-hand side of Equation 18 can be found in Monaghan et al. [27] where r _ij is the magnitude of the vector r _ij . If the reactive Robin BC on reactive boundary is imposed, Equation 9 is then solved instead and its corresponding SPH discretization form is:

$$\begin{array}{@{}rcl@{}} \frac{d p_{i}}{dt} &=& \sum_{j\in \Omega\cup\Omega_{b}} \frac{D_{i}+D_{j}}{d_{j}} (p_{i}-p_{j}) \frac{1}{r_{ij}} \frac{d w(r_{ij}, h)}{d r_{ij}} \\ &&+\beta \sum_{j\in \Omega} \frac{D_{i} p_{i}+D_{j} p_{j}}{d_{j}} (W_{i}-W_{j}) \frac{1}{r_{ij}} \frac{d w(r_{ij}, h)}{d r_{ij}} \\ &&- \alpha p_{i} \sum_{k\in \Omega_{a}} \frac{\mathbf{n}_{i} + \mathbf{n}_{k}}{d_{k}}\cdot \frac{\mathbf{r}_{ik}}{r_{ik}} \frac{d w(r_{ik}, h_{r})}{d r_{ik}}. \end{array} $$

((19))

To integrate the SPH Eqs. 18 and 19, an explicit Verlet scheme [38] is employed. The last term in Eq. 19 is obtained by discretizing the integral in Eq. 9 as a Riemann sum:

$$ \begin{aligned} {} \alpha p(\mathbf{x}, t) &\iiint\limits_{\Omega_{a}} [\mathbf{n}(\mathbf{x}) + \mathbf{n}(\mathbf{x}')]\cdot \nabla_{\mathbf{x}} w(\mathbf{x}-\mathbf{x}', h_{r}) d\mathbf{x}' \\ &= \alpha p(\mathbf{x}, t) \sum_{k \in \Omega_{a}} \Delta V_{k} [\mathbf{n}(\mathbf{x}) + \mathbf{n}_{k}] \cdot \nabla_{\mathbf{x}} w(\mathbf{x} - \mathbf{r}_{k}, h_{r}), \end{aligned} $$

((20))

where $\Delta V_{k} = \frac {1}{d_{k}}$ is the volume of particle k and $\sum _{k \in \Omega _{a}}$ is the summation over the reactive boundary particles within 2h _r distance from particle i. The SPH expression for calculating the normal unit vector is obtained as:

$$ \mathbf{n}_{i} = \frac{\sum\limits_{j \in \Omega\cup\Omega_{a}}\frac{1}{d_{j}}(\phi_{j}-\phi_{i})\nabla_{i} w(\mathbf{r}_{ij}, h_{r})}{\left |\sum\limits_{j \in \Omega\cup\Omega_{a}}\frac{1}{d_{j}}(\phi_{j}-\phi_{i})\nabla_{i} w(\mathbf{r}_{ij}, h_{r}) \right |}. $$

((21))

In the simulations presented below, we set h _r =h but it could be set differently in future applications.

Note that the reflective Neumann BC (Equation 7 or 10) can be simply enforced in SPH by excluding the contribution from the boundary particles in the summation. The Dirichlet BC (Equation 4 or 6) is enforced by assigning the fixed boundary value of concentration on the boundary particles. If the Robin BC is imposed, the reaction rate k _on(t) can be calculated by Equation 42 in Appendix B and its corresponding SPH discretization form is:

$$ k_{\text{on}}= \sum_{i \in \Omega} \alpha p_{i} \left[\sum_{k\in \Gamma_{a}} \frac{\mathbf{n}_{i} + \mathbf{n}_{k}}{d_{k}}\cdot \frac{\mathbf{r}_{ik}}{r_{ik}} \frac{d w(r_{ik}, h)}{d r_{ik}}\right]. $$

((22))

Otherwise, when the Dirichlet BC is enforced on the reactive boundary, the discretization of Eq. 41 in Appendix B is:

$$ k_{\text{on}}= \sum_{i \in \Omega} (\mathbf{n}_{i} \cdot \mathbf{J}_{i}) \left[\sum_{k\in \Gamma_{a}} \frac{\mathbf{n}_{i} + \mathbf{n}_{k}}{d_{k}}\cdot \frac{\mathbf{r}_{ik}}{r_{ik}} \frac{d w(r_{ik}, h)}{d r_{ik}}\right], $$

((23))

where

$$ \begin{aligned} \mathbf{J}_{i} &= D_{i} \sum_{j\in \Omega\cup\Omega_{b}\cup\Omega_{a}} \frac{p_{j}-p_{i}}{d_{j}}\frac{\mathbf{r}_{ij}}{r_{ij}} \frac{d w(r_{ij}, h)}{d r_{ij}}\\ &\quad+ \beta D_{i} p_{i} \sum_{j\in \Omega} \frac{W_{j}-W_{i}}{d_{j}} \frac{\mathbf{r}_{ij}}{r_{ij}} \frac{d w(r_{ij}, h)}{d r_{ij}}. \end{aligned} $$

((24))

Calculation of potentials of mean force

We calculated the potential of mean force W(x) using the recently published 2.1 Å resolution structure of mAChE [39]. To prepare this structure for the calculation, we assigned titration states of ionizable residues using PROPKA [40] at pH 7, and we used PDB2PQR v1.8 [41,42] to assign atomic radii and charges. APBS v1.4 was used to perform a nonlinear Poisson-Boltzmann multi-grid calculation of the electrostatic potential over the entire simulation domain [43]. The small and large domains were set to 600 Å and 400 Å, respectively, with a fine grid spacing of 0.600 Å. For APBS calculations, we used the single Debye-Hückel boundary condition, a smoothed molecular surface, and protein and solvent dielectrics of 2 and 78.54, respectively. Atomic charges were mapped onto the grids using cubic B-spline discretization. The calculated potential was mapped onto the SPH discretization points of protein and solvent via trilinear interpolation.

Appendix A: continuum surface reaction method

In this appendix, we present a detailed derivation of the continuum surface reaction method for solving the Smoluchowski equation subject to Robin BC. We start from a two-sided problem; i.e., the concentration field p(x,t) is extended into the sub-domain Ω _a that is separated from Ω by Γ _a such that Eq. 2 can be approximated as

$$ \begin{aligned} \frac{dp^{r}(\mathbf{x},t)}{dt} &= \nabla\cdot \left(D (\mathbf{x})[\nabla p^{r}(\mathbf{x},t) + \beta p^{r}(\mathbf{x},t) \nabla W(\mathbf{x})]\right)\\ &\quad- P_{\Omega}(\mathbf{x},t) ~ \text{for} ~ \mathbf{x}\in\Omega\cup\Omega_{a} \end{aligned} $$

((25))

subject to

$$ \mathbf{n}(\mathbf{x}_{s}) \cdot [\mathbf{J}^{r}(\mathbf{x}_{s},t)|_{\mathbf{x}_{s}{\in\Gamma^{F}_{a}}}-\mathbf{J}^{r}(\mathbf{x}_{s},t)|_{\mathbf{x}_{s}{\in\Gamma^{S}_{a}}}]= 0 ~ \text{for} ~ \mathbf{x}_{s}\in\Gamma_{a}, $$

((26))

where ${\Gamma ^{F}_{a}}$ and ${\Gamma ^{S}_{a}}$ are the two sides of Γ _a , respectively. The boundary condition Eq. 26 emphasizes that the extended concentration field is continuous across Γ _a . Comparison of the weak formulations of Eq. 2 subject to Eq. 5 and Eq. 25 subject to Eq. 26 yields the relationship

$$ \iiint\limits_{\Omega\cup\Omega_{a}} P_{\Omega} (\mathbf{x},t) d\mathbf{x} = \iint\limits_{\Gamma_{a}} \alpha p(\mathbf{x}'_{s},t) d \mathbf{x}'_{s}. $$

((27))

This weak formulation is obtained by integrating the equations over their respected domains and then applying Gauss’ theorem with the corresponding boundary conditions. To derive the formulation of P _Ω , we define a color function (i.e., a sharp characteristic function) as:

$$ \phi(\mathbf{x}) = \begin{cases} 0, & \mathbf{x}\in\Omega,\\ 1, & \mathbf{x}\in\Omega_{a}. \end{cases} $$

((28))

and its smooth counterpart as

$$ \tilde{\phi}(\mathbf{x}) = \iiint\limits_{\Omega\cup\Omega_{a}} \phi(\mathbf{x}') w(\mathbf{x}-\mathbf{x}', h_{r}) d\mathbf{x}^{\prime}. $$

((29))

The gradient of $\tilde {\phi }$ can then be found from Eq. 29 as

$$ \nabla \tilde{\phi}(\mathbf{x}) = \iiint\limits_{\Omega\cup\Omega_{a}} \phi(\mathbf{x}') \nabla_{\mathbf{x}} w(\mathbf{x}-\mathbf{x}', h_{r}) d\mathbf{x}^{\prime}. $$

((30))

Using the definition of the surface delta function [44]:

$$ {} \delta [\mathbf{n}(\mathbf{x}_{s}) \cdot (\mathbf{x}_{s}-\mathbf{x})] = \mathbf{n}(\mathbf{x}) \cdot \nabla {\phi}(\mathbf{x}), ~~~ \mathbf{x} \in \Omega\cup\Omega_{a}, ~\text{for}~ \mathbf{x}_{s}\in\Gamma_{a}, $$

((31))

and noting that

$$ {\lim}_{h_{r}\to 0} \tilde{\phi} = \phi, $$

((32))

we can rewrite the surface delta function in terms of $\tilde {\phi }$ as:

$$ {} \delta [\mathbf{n}(\mathbf{x}_{s}) \cdot (\mathbf{x}_{s}-\mathbf{x})] \,=\, \mathbf{n}(\mathbf{x}) \cdot {\lim}_{h_{r}\to 0} \nabla \tilde{\phi}(\mathbf{x}), \text{for}\, \mathbf{x} \!\in\! \Omega\cup\Omega_{a}, \, \mathbf{x}_{s} \!\in\! \Gamma_{a}. $$

((33))

The surface integral can then be rewritten as a volume integral through the surface delta function:

$${\fontsize{8.3pt}{9.6pt}\selectfont{ \begin{aligned} {} \iint\limits_{\Gamma_{a}} \alpha p(\mathbf{x}'_{s}, t) d\mathbf{x}'_{s} &= \iiint\limits_{\Omega\cup\Omega_{a}} \alpha p(\mathbf{x},t) \delta [\mathbf{n}(\mathbf{x}_{s}) \cdot (\mathbf{x}_{s}-\mathbf{x})] d\mathbf{x}, ~\text{for}~ \mathbf{x}_{s}\in\Gamma_{a}, \\ & = \iiint\limits_{\Omega\cup\Omega_{a}} \alpha p^{r}(\mathbf{x},t) \mathbf{n}(\mathbf{x}) \cdot \nabla \tilde{\phi}(\mathbf{x}) d\mathbf{x}. \end{aligned}}} $$

((34))

To uniquely define P _Ω (x,t), we require it to vanish at a normal distance greater than h _r from Γ _a and require that

$$ {\lim}_{h_{r}\to 0}\iiint\limits_{\Omega\cup\Omega_{a}} P_{\Omega} (\mathbf{x},t) d\mathbf{x} =\iint\limits_{\Gamma_{a}} \alpha p(\mathbf{x}'_{s},t) d\mathbf{x}'_{s}. $$

((35))

Comparing Eqs. 34 and 35 yields an expression for P _Ω (x,t) as:

$$ P_{\Omega}(\mathbf{x},t) = \alpha p^{r}(\mathbf{x},t) \mathbf{n}(\mathbf{x}) \cdot \nabla \tilde{\phi}(\mathbf{x}), ~\text{for}~ \mathbf{x} \in \Omega\cup\Omega_{a}. $$

((36))

Eq. 25 can then be rewritten by combining Eqs. 36 and 30 as:

$$ \begin{aligned} {} \frac{dp^{r}(\mathbf{x},t)}{dt} &= \nabla\cdot \left(D (\mathbf{x})[\nabla p^{r}(\mathbf{x},t) + \beta p^{r}(\mathbf{x},t) \nabla W(\mathbf{x})]\right)\\ &\quad- \!\alpha p^{r} (\mathbf{x}, t)\! \iiint\limits_{\Omega\cup\Omega_{a}}\! \mathbf{n}(\mathbf{x}) \cdot [\phi(\mathbf{x}') \nabla_{\mathbf{x}} w(\mathbf{x}\,-\,\mathbf{x}', h_{r})]\\ &\quad \times d\mathbf{x}', \text{for}~\mathbf{x}\!\in\!\Omega\cup\Omega_{a}. \end{aligned} $$

((37))

Since p ^r is not uniquely defined on Ω _a , we introduce a one-sided formulation by approximating Eq. 37 as:

$${\fontsize{8.7pt}{9.6pt}\selectfont{ \begin{aligned} {} \frac{dp^{r}(\mathbf{x},t)}{dt} &= \nabla\cdot \left(D (\mathbf{x})[\nabla p^{r}(\mathbf{x},t) + \beta p^{r}(\mathbf{x},t) \nabla W(\mathbf{x})]\right)\\ &\quad- \!\alpha p^{r} (\mathbf{x}, t)\! \iiint\limits_{\Omega\cup\Omega_{a}}\! [\!\mathbf{n}(\mathbf{x})+\mathbf{n}(\mathbf{x})] \cdot [\phi(\mathbf{x}') \nabla_{\mathbf{x}} w(\mathbf{x}-\mathbf{x}', h_{r})]\\ &\quad \times d\mathbf{x}', ~\text{for}~\mathbf{x}\in\Omega, \end{aligned}}} $$

((38))

subject to the reflective Neumann BC (Eq. 10). Note that ϕ is non-zero only in Ω _a , where it is equal to 1 as defined in Eq. 28. Thus, the modified governing equation takes its final form as Eq. 9.

Appendix B: calculation of reaction rate

Similar to the derivation in Eq. 34, using the definition of the surface delta function and given p _bulk=1, the expression for the reaction rate can be rewritten as

$$ k_{\text{on}}=\iiint\limits_{\Omega\cup\Omega_{a}} [\mathbf{n}(\mathbf{x}) \cdot \mathbf{J} (\mathbf{x},t)] [\mathbf{n}(\mathbf{x}) \cdot \nabla \tilde{\phi}(\mathbf{x})] d\mathbf{x}. $$

((39))

Substituting Eq. 30 into the above equation and using Eq. 28, a new expression of k _on can be obtained:

$$ {} k_{\text{on}}=\iiint\limits_{\Omega\cup\Omega_{a}} [\mathbf{n}(\mathbf{x}) \cdot \mathbf{J} (\mathbf{x},t)] \iiint\limits_{\Omega_{a}} \mathbf{n}(\mathbf{x}) \cdot \nabla_{\mathbf{x}} w(\mathbf{x}-\mathbf{x}', h_{r}) d\mathbf{x}' d\mathbf{x}. $$

((40))

Similar to Eq. 38, the corresponding one-sided formulation is:

$${\fontsize{8.1pt}{9.6pt}\selectfont{ \begin{aligned} {} k_{\text{on}}=\iiint\limits_{\Omega} [\mathbf{n}(\mathbf{x}) \cdot \mathbf{J} (\mathbf{x},t)] \iiint\limits_{\Omega_{a}} [\mathbf{n}(\mathbf{x}) + \mathbf{n}(\mathbf{x}')] \cdot \nabla_{\mathbf{x}} w(\mathbf{x}-\mathbf{x}', h_{r}) d\mathbf{x}' d\mathbf{x}. \end{aligned}}} $$

((41))

If the Robin BC (Eq. 5) is enforced, Eq. 41 can be reduced to

$$ {} k_{\text{on}} \,=\,\! \iiint\limits_{\Omega}\!\! \alpha p^{r}(\mathbf{x},t) \!\!\iiint\limits_{\Omega_{a}} \! [\!\mathbf{n}(\mathbf{x}) + \mathbf{n}(\mathbf{x}')]\! \cdot \nabla_{\mathbf{x}} w(\mathbf{x}-\mathbf{x}', h_{r}) d\mathbf{x}' d\mathbf{x}. $$

((42))