Skip to main content

Advertisement

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

BMC Biophysics volume 8, Article number: 7 (2015) Cite this article

Article metrics

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 1. The entire domain was confined by the outer boundary Γ b determined as a spherical surface with the radius of R 2=125 Å.

For the first test case, we let R 1=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 2 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
figure1

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

Full size image

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 1=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 1). The reactive coefficient for the Robin BC is α=1×103. 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×1012 M −1 m i n −1 for the Dirichlet BC, and is 8.24×1011 M −1 m i n −1 for the Robin BC, both with L 2 errors less than 3% relative to the analytically evaluated ones.

Figure 2
figure2

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

Full size image
Figure 3
figure3

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

Full size image

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 2=128 Å. Figure 4B and 4C depict the constructed reactive boundaries 1 and 4.

Figure 4
figure4

Panel A shows the discretized domain with R 2=128 Å and the mAChE molecule in the center with the reactive boundary shown in purple. Light blue indicates the outer boundary (R 2), 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.

Full size image

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
figure5

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

Full size image
Figure 6
figure6

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

Full size image

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
figure7

(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.

Full size image

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×103 which optimized agreement between computational and experimental results.

Figure 8
figure8

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.

Full size image

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 −1 m i n −1. With the Dirichlet BC, SPH predicts k on with RMSD of 0.57 M −1 m i n −1, intermediate between FEM and BD results. However, with the Robin BC, SPH predicts k on with RMSD of 0.33 M −1 m i n −1, better than the FEM and BD results.

Figure 9
figure9

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×103. For standardization, both computed and experimental data are fitted to the Debye-Hückel limiting law.

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

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 2≈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
figure10

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 .

Full size image

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))

References

  1. 1

    Eigen, M. 1974. Diffusion control in biochemical reactions. In Quantum statistical mechanics in the natural sciencespp, ed. Mintz SL and Widemeyer SM. 37–61. New York: Plenum.

  2. 2

    Bazelyansky, M, E Robey, and JF Kirsch. 1986. Fractional diffusion-limited component of reactions catalyzed by acetylcholinesterase. Biochem.25(1): 125–30. doi:10.1021/bi00349a019.

  3. 3

    Nolte, HJ, TL Rosenberry, and E Neumann. 1980. Effective charge on acetylcholinesterase active-sites determined from the ionic-strength dependence of association rate constants with cationic ligands. Biochem.19(16): 3705–11. doi:10.1021/bi00557a011.

  4. 4

    Radic, Z, PD Kirchhoff, DM Quinn, JA McCammon, and P Taylor. 1997. Electrostatic influence on the kinetics of ligand binding to acetylcholinesterase - Distinctions between active center ligands and fasciculin. J Biol Chem.272(37): 23265–77. doi:10.1074/jbc.272.37.23265.

  5. 5

    Tan, RC, TN Truong, JA McCammon, and JL Sussman. 1993. Acetylcholinesterase - electrostatic steering increases the rate of ligand-binding. Biochem.32(2): 401–3. doi:10.1021/bi00053a003.

  6. 6

    Song, YH, YJ Zhang, TY Shen, CL Bajaj, A McCammon, and NA Baker. 2004. Finite element solution of the steady-state Smoluchowski equation for rate constant calculations. Biophys J.86(4): 2017–29.

  7. 7

    Song, YH, YJ Zhang, CL Bajaj, and NA Baker. 2004. Continuum diffusion reaction rate calculations of wild-type and mutant mouse acetylcholinesterase: Adaptive finite element analysis. Biophys J.87(3): 1558–66. doi:10.1529/biophysj.104.041517.

  8. 8

    Wade, RC, ME Davis, BA Luty, JD Madura, and JA McCammon. 1993. Gating of the active-site of triose phosphate isomerase - Brownian dynamics simulations of flexible peptide loops in the enzyme. Biophys J.64(1): 9–15.

  9. 9

    Tara, S, AH Elcock, PD Kirchhoff, JM Briggs, Z Radic, P Taylor, and al et. 1998. Rapid binding of a cationic active site inhibitor to wild type and mutant mouse acetylcholinesterase: Brownian dynamics simulation including diffusion in the active site gorge. Biopolymers.46(7): 465–74.

  10. 10

    Davis, ME, JD Madura, BA Luty, and JA McCammon. 1991. Electrostatics and diffusion of molecules in solution - Simulations with the University-of-Houston-Brownian dynamics program. Comput Phys Commun.62(2–3): 187–97.

  11. 11

    Davis, ME, JD Madura, J Sines, BA Luty, SA Allison, and JA McCammon. 1991. Diffusion-controlled enzymatic-reactions. Methods Enzymology.202: 473–97.

  12. 12

    Genest, D. 1989. A Monte-Carlo simulation study of the influence of internal motions on the molecular-conformation deduced from two-dimensional NMR experiments. Biopolym.28(11): 1903–11. doi:10.1002/bip.360281107.

  13. 13

    Saxton, MJ. 1992. Lateral diffusion and aggregation - a Monte-Carlo study. Biophys J.61(1): 119–28.

  14. 14

    Kerr, RA, TM Bartol, B Kaminsky, M Dittrich, J-CJ Chang, SB Baden, and al et. 2008. Fast Monte Carlo simulation methods for biological reaction-diffusion systems in solution and on surfaces. SIAM J Sci Comput.30(6): 3126–49.

  15. 15

    Northrup, SH, JO Boles, and JCL Reynolds. 1988. Brownian dynamics of cytochrome-c and cytochrome-c peroxidase association. Sci.241(4861): 67–70. doi:10.1126/science.2838904.

  16. 16

    Eastman, P, and S Doniach. 1998. Multiple time step diffusive Langevin dynamics for proteins. Proteins-Structure Funct Genet.30(3): 215–27.

  17. 17

    Yeomans-Reyna, L, and M Medina-Noyola. 2001. Self-consistent generalized Langevin equation for colloid dynamics. Phys Rev E. 64(6). Article No. 066114. doi:06611410.1103/PhysRevE.64.066114.

  18. 18

    Baker, NA.2005. Biomolecular applications of Poisson-Boltzmann methods. Reviews in Computational Chemistry(KB Lipkowitz, R Larter, and TR Cundari, eds.), Vol. 21. New York: Wiley-VCH, Inc.

  19. 19

    Baker, NA. 2005. Improving implicit solvent simulations: a Poisson-centric view. Curr Opinion Struct Biol.15(2): 137–43.

  20. 20

    Senapati, S, CF Wong, and JA McCammon. 2004. Finite concentration effects on diffusion-controlled reactions. J Chem Phys.121(16): 7896–900.

  21. 21

    Kurnikova, MG, RD Coalson, P Graf, and A Nitzan. 1999. A lattice relaxation algorithm for three-dimensional Poisson-Nernst-Planck theory with application to ion transport through the gramicidin A channel. Biophys J.76(2): 642–56. doi:10.1016/s0006-3495(99)77232-2.

  22. 22

    Smart, JL, and JA McCammon. 1998. Analysis of synaptic transmission in the neuromuscular junction using a continuum finite element model. Biophys J.75(4): 1679–88.

  23. 23

    Cheng, Y, JK Suen, D Zhang, SD Bond, Y Zhang, Y Song, and et al.2007. Finite element analysis of the time-dependent Smoluchowski equation for acetylcholinesterase reaction rate calculations. Biophys J.92(10): 3397–406. doi:10.1529/biophysj.106.102533.

  24. 24

    Lu, B, MJ Holst, JA McCammon, and YC Zhou. 2010. Poisson-Nernst-Planck equations for simulating biomolecular diffusion-reaction processes I: Finite element solutions. J Comput Phys.229(19): 6979–94. doi:10.1016/j.jcp.2010.05.035.

  25. 25

    Lu, B, and JA McCammon. 2010. Kinetics of diffusion-controlled enzymatic reactions with charged substrates. PMC Biophys.3: 1. doi:10.1186/1757-5036-3-1.

  26. 26

    Holst, MJ. 2001. Adaptive numerical treatment of elliptic systems on manifolds. Adv Comput Mathematics.15(1): 139–91. doi:10.1023/a:1014246117321.

  27. 27

    Monaghan, JJ. 2005. Smoothed particle hydrodynamics. Rep Prog Phys.68: 1703–59.

  28. 28

    Monaghan, JJ. 2012. Smoothed particle hydrodynamics and its diverse applications. Ann Rev Fluid Mech.44: 323–46.

  29. 29

    Pan, W, AM Tartakovsky, and JJ Monaghan. 2013. Smoothed particle hydrodynamics non-Newtonian model for ice sheet and ice shelf dynamics. J Comput Phys.242(1): 828–42.

  30. 30

    Pan, W, D Li, AM Tartakovsky, S Ahzi, M Khraisheh, and M Khaleel. 2013. A new smoothed particle hydrodynamics non-Newtonian model for friction stir welding: Process modeling and simulation of microstructure evolution in a magnesium alloy. Int J Plasticity.48: 189–204.

  31. 31

    Ryan, EM, AM Tartakovsky, and C Amon. 2010. A novel method for modeling Neumann and Robin boundary conditions in smoothed particle hydrodynamics. Comput Phys Commun.181: 2008–23.

  32. 32

    Pan, W, J Bao, and AM Tartakovsky. 2014. Smoothed particle hydrodynamics continuous boundary force method for Navier-Stokes equations subject to a Robin boundary condition. J Comput Phys.259: 242–59.

  33. 33

    Gillespie, DT. 2007. Stochastic simulation of chemical kinetics. Annu Rev Phys Chem.58: 35–55.

  34. 34

    Beri, V, JT Auletta, GM Maharvi, JF Wood, AH Fauq, and TL Rosenberry. 2013. Hydrolysis of low concentrations of the acetylthiocholine analogs acetyl(homo)thiocholine and acetyl(nor)thiocholine by acetylcholinesterase may be limited by selective gating at the enzyme peripheral site. Chem Biol Interactions.203(1): 38–43. doi:10.1016/j.cbi.2012.09.017.

  35. 35

    Zhou, HX, ST Wlodek, and JA McCammon. 1998. Conformation gating as a mechanism for enzyme specificity. Proc Natl Acad Sci U S A.95(16): 9280–3.

  36. 36

    Baker, NA, and JA McCammon. 1999. Non-Boltzmann rate distributions in stochastically gated reactions. J Phys Chem B.103(4): 615–7.

  37. 37

    Zhou, HX. 1990. On the calculation of diffusive reaction rates using Brownian dynamics simulations. J Chem Phys.92(5): 3092–95.

  38. 38

    Allen, M, and D Tildesley.1989. Computer Simulation of Liquids. Oxford: Clarendon Press.

  39. 39

    Andersson, CD, N Forsgren, C Akfur, A Allgardsson, L Berg, C Engdahl, and et al.2013. Divergent structure-activity relationships of structurally similar acetylcholinesterase inhibitors. J Med Chem.56(19): 7615–24. doi:10.1021/jm400990p.

  40. 40

    Li, H, AD Robertson, and JH Jensen. 2005. Very fast empirical prediction and rationalization of protein p K a values. Proteins Structure Funct Bioinformatics.61(4): 704–21. doi:10.1002/prot.20660.

  41. 41

    Dolinsky, TJ, P Czodrowski, H Li, JE Nielsen, JH Jensen, G Klebe, and et al.2007. PDB2PQR: expanding and upgrading automated preparation of biomolecular structures for molecular simulations. Nucleic Acids Res.35: 522–5. doi:10.1093/nar/gkm276.

  42. 42

    Dolinsky, TJ, JE Nielsen, JA McCammon, and NA Baker. 2004. PDB2PQR: an automated pipeline for the setup of poisson-boltzmann electrostatics calculations. Nucleic Acids Res.32: 665–7. doi:10.1093/nar/gkh381.

  43. 43

    Baker, NA, D Sept, S Joseph, MJ Holst, and JA McCammon. 2001. Electrostatics of nanosystems: Application to microtubules and the ribosome. Proc Natl Acad Sci U S A.98(18): 10037–41. doi:10.1073/pnas.181342398.

  44. 44

    Lange, R-J. 2012. Potential theory, path integrals and the Laplacian of the indicator. J High Energy Phys.2012(11): 1–46.

Download references

Acknowledgments

The authors gratefully acknowledge the funding support by the Applied Mathematics Program within the U.S. Department of Energy Office of Advanced Scientific Computing Research as part of the Collaboration on Mathematics for Mesoscopic Modeling of Materials (CM4). The research was performed using PNNL Institutional Computing at Pacific Northwest National Laboratory. The Pacific Northwest National Laboratory is operated for the U.S. Department of Energy by Battelle under Contract DE-AC06-76RL01830. This research was also supported by National Institutes of Health grant GM069702 to NAB.

Author information

Affiliations

  1. Advanced Computing, Mathematics and Data Division, Pacific Northwest National Laboratory, MSID K7-90, 99352, Richland, PO Box 999, WA, USA

    • Wenxiao Pan

  2. Chemical Physics and Analysis Division, Mathematics and Data Division, Pacific Northwest National Laboratory, MSID K1-83, 99352, Richland, PO Box 999, WA, USA

    • Michael Daily

  3. Computational and Statistical Analytics Division, Pacific Northwest National Laboratory, MSID K7-20, 99352, Richland, PO Box 999, WA, USA

    • Nathan A Baker

Authors

  1. Search for Wenxiao Pan in:

  2. Search for Michael Daily in:

  3. Search for Nathan A Baker in:

Corresponding author

Correspondence to Nathan A Baker.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

WP implemented the new SPH methods, performed simulations, and participated in the writing of the manuscript. MD performed simulations, and participated in the writing of the manuscript. NAB conceived of the study, participated in its design and coordination, and helped to draft the manuscript. All authors read and approved the final manuscript.

Rights and permissions

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Pan, W., Daily, M. & Baker, N.A. Numerical calculation of protein-ligand binding rates through solution of the Smoluchowski equation using smoothed particle hydrodynamics. BMC Biophys 8, 7 (2015) doi:10.1186/s13628-015-0021-y

Download citation

Keywords

  • Diffusion
  • Smoluchowski equation
  • Smoothed particle hydrodynamics
  • Protein-ligand interactions
  • Binding rates
  • Acetylcholinesterase

BMC Biophysics

ISSN: 2046-1682