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.
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 R1. The entire domain was confined by the outer boundary Γb determined as a spherical surface with the radius of R2=125 Å.
For the first test case, we let R1=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 L2 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=R2 at different times with the relative L2 errors for different resolutions. Specifically, L2=0.0326 for Δx=8Å(green square), L2=0.0180 for Δx=4Å(blue circle), and L2=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 R1=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=R1). 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×1012M−1min−1 for the Dirichlet BC, and is 8.24×1011M−1min−1 for the Robin BC, both with L2 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=R1 and r=R2 at steady-state with the relative L2 errors for different resolutions. Specifically, L2=0.0666 for Δx=8Å(green square), L2=0.0321 for Δx=4Å(blue circle), and L2=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=R1 and Dirichlet BC on r=R2 at steady-state with the relative L2 errors for different resolutions. Specifically, L2=0.00914 for Δx=8Å(green square), L2=0.00598 for Δx=4Å(blue circle), and L2=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 R2=128 Å. Figure 4B and 4C depict the constructed reactive boundaries 1 and 4.
Figure 4
Panel A shows the discretized domain with R2=128 Å and the mAChE molecule in the center with the reactive boundary shown in purple. Light blue indicates the outer boundary (R2), 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×103) 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×103) 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 kon(t) on reactive boundary 1 at different ionic strengths. For this boundary, kon(t) converges within 150 ns for all ionic strengths. The right panel shows kon(t) on reactive boundaries 1-4, respectively, at 0.15 M ionic strength.
Figure 7
(Left) kon 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) kon 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
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 kon calculated at 0.67 M ionic strength, zE is the effective enzyme charge, and zI 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
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 kon, while the FEM produces good agreement with experimental kon at RMSD = 0.37 M−1min−1. With the Dirichlet BC, SPH predicts kon with RMSD of 0.57 M−1min−1, intermediate between FEM and BD results. However, with the Robin BC, SPH predicts kon with RMSD of 0.33 M−1min−1, 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×103. 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
RMSD exp of simulation results to experimental kon is calculated over the range of ionic strengths between 0 and 0.67 M. The unit of \(\textit {k}_{\text {on}}^{\text {0}}\), \(\textit {k}_{\text {on}}^{\text {H}}\) and RMSD is 1011M-1min-1. And the error is the standard deviation of parameter fits using nonlinear least squares.
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 kon 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:
where D(x) is the diffusion coefficient; for simplicity, it is assumed to be constant. β=1/kBT is the inverse Boltzmann energy with the Boltzmann constant kB 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,
specifying a Dirichlet BC on the outer boundary Γb where the concentration is equal to a bulk concentration pbulk. 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 R2≈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 pbulk=1.
The active site boundary Γa was modeled using either reactive Robin or absolute absorbing Dirichlet BC:
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
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]:
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:
In Eq. 9, the normalized kernel function, w(x), is a positive bell-shaped function with at least first continuous derivative and compact support κhr such that w(|r|>κhr)=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:
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 ri (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 di as the prescribed number density at ri. The discretization is based on a meshless interpolation scheme:
where, Ai=A(ri) is a function defined at particle i, Aj=A(rj) is the function defined at particle i’s neighboring particles j with distances rij=ri−rj, and w(rij,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 ri have a non-zero contribution to the summation. Spatial derivatives of A can be calculated as
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:
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 rij is the magnitude of the vector rij. If the reactive Robin BC on reactive boundary is imposed, Equation 9 is then solved instead and its corresponding SPH discretization form is:
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:
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 2hr distance from particle i. The SPH expression for calculating the normal unit vector is obtained as:
In the simulations presented below, we set hr=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 kon(t) can be calculated by Equation 42 in Appendix B and its corresponding SPH discretization form is:
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
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
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:
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 pbulk=1, the expression for the reaction rate can be rewritten as
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.
Authors’ Affiliations
(1)
Advanced Computing, Mathematics and Data Division, Pacific Northwest National Laboratory
(2)
Chemical Physics and Analysis Division, Mathematics and Data Division, Pacific Northwest National Laboratory
(3)
Computational and Statistical Analytics Division, Pacific Northwest National Laboratory
References
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.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleADSGoogle Scholar
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.View ArticleADSGoogle Scholar
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.View ArticleADSGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleADSGoogle Scholar
Davis, ME, JD Madura, J Sines, BA Luty, SA Allison, and JA McCammon. 1991. Diffusion-controlled enzymatic-reactions. Methods Enzymology.202: 473–97.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
Saxton, MJ. 1992. Lateral diffusion and aggregation - a Monte-Carlo study. Biophys J.61(1): 119–28.View ArticleADSGoogle Scholar
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.View ArticleMATHMathSciNetGoogle Scholar
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.View ArticleADSGoogle Scholar
Eastman, P, and S Doniach. 1998. Multiple time step diffusive Langevin dynamics for proteins. Proteins-Structure Funct Genet.30(3): 215–27.View ArticleGoogle Scholar
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.Google Scholar
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.Google Scholar
Senapati, S, CF Wong, and JA McCammon. 2004. Finite concentration effects on diffusion-controlled reactions. J Chem Phys.121(16): 7896–900.View ArticleADSGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleADSGoogle Scholar
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.View ArticleADSMATHMathSciNetGoogle Scholar
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.View ArticleGoogle Scholar
Holst, MJ. 2001. Adaptive numerical treatment of elliptic systems on manifolds. Adv Comput Mathematics.15(1): 139–91. doi:10.1023/a:1014246117321.View ArticleMATHMathSciNetGoogle Scholar
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.View ArticleADSMATHMathSciNetGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleADSMATHMathSciNetGoogle Scholar
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.View ArticleADSMathSciNetGoogle Scholar
Gillespie, DT. 2007. Stochastic simulation of chemical kinetics. Annu Rev Phys Chem.58: 35–55.View ArticleADSGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleADSGoogle Scholar
Baker, NA, and JA McCammon. 1999. Non-Boltzmann rate distributions in stochastically gated reactions. J Phys Chem B.103(4): 615–7.View ArticleGoogle Scholar
Zhou, HX. 1990. On the calculation of diffusive reaction rates using Brownian dynamics simulations. J Chem Phys.92(5): 3092–95.View ArticleADSGoogle Scholar
Allen, M, and D Tildesley.1989. Computer Simulation of Liquids. Oxford: Clarendon Press.Google Scholar
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.View ArticleGoogle Scholar
Li, H, AD Robertson, and JH Jensen. 2005. Very fast empirical prediction and rationalization of protein p Ka values. Proteins Structure Funct Bioinformatics.61(4): 704–21. doi:10.1002/prot.20660.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleGoogle Scholar
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.View ArticleADSGoogle Scholar
Lange, R-J. 2012. Potential theory, path integrals and the Laplacian of the indicator. J High Energy Phys.2012(11): 1–46.View ArticleADSGoogle Scholar
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
