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

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 diffusionlimited behavior, complex geometry, and strong electrostatic influence has made AChE a useful target for both experimental and computational studies of biomolecular diffusion [4][5][6][7][8][9][10][11]. Two major classes of methods have been used to estimate diffusion rates in biomolecular systems. Mesoscopic coarse-grained methods like Monte Carlo [12][13][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][22][23][24][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][29][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][22][23][24][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.

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 Å.
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 × 10 3 . 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 12 M −1 min −1 for the Dirichlet BC, and is 8.24 × 10 11 M −1 min −1 for the Robin BC, both with L 2 errors less than 3% relative to the analytically evaluated ones.

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-  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.
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 nonzero 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.
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.
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 0 on is the effective reaction rate at zero ionic strength rate, k H on 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 3 which optimized agreement between computational and experimental results. 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. 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 0 on of 2.23

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(−α t 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].

Governing equation and boundary conditions
The time-dependent Smoluchowski equation can be written as: dp(x, t) dt where p(x, t) is the concentration distribution of the reactants, and the concentration flux J(x, t) is defined as: 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, 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 3 . For standardization, both computed and experimental data are fitted to the Debye-Hückel limiting law. 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: or 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. 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: subject to the reflective Neumann BC: The derivation of Eq. 9 is detailed in Appendix A, which demonstrates 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: and The normal unit vector n in Eq. 9 can be found in terms of a smoothed color functionφ as defined in Appendix A:

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: 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 In the present work, we use a cubic spline kernel as the weighting function 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 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: 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 V k = 1 d k is the volume of particle k and k∈ 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: 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: Otherwise, when the Dirichlet BC is enforced on the reactive boundary, the discretization of Eq. 41 in Appendix B is: where

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 subject to where F a and 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: and its smooth counterpart as The gradient ofφ can then be found from Eq. 29 as Using the definition of the surface delta function [44]: δ[ n(x s )·(x s −x)] = n(x)·∇φ(x), x ∈ ∪ a , for x s ∈ a , (31) and noting that we can rewrite the surface delta function in terms ofφ as: ∇φ(x), for x ∈ ∪ a , x s ∈ a .
The surface integral can then be rewritten as a volume integral through the surface delta function: To uniquely define P (x, t), we require it to vanish at a normal distance greater than h r from a and require that Comparing Eqs. 34 and 35 yields an expression for P (x, t) as: P (x, t) = αp r (x, t)n(x) · ∇φ(x), for x ∈ ∪ a . (36) Eq. 25 can then be rewritten by combining Eqs. 36 and 30 as: Since p r is not uniquely defined on a , we introduce a one-sided formulation by approximating Eq. 37 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.