This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Song, Y. Articles by Baker, N. A. Search for Related Content PubMed PubMed Citation Articles by Song, Y. Articles by Baker, N. A.

Finite Element Solution of the Steady-State Smoluchowski Equation for Rate Constant Calculations

Yuhua Song ^*, Yongjie Zhang , Tongye Shen , Chandrajit L. Bajaj , J. Andrew McCammon ^¶ and Nathan A. Baker ^*

^* Department of Biochemistry and Molecular Biophysics, Center for Computational Biology, Washington University in St. Louis, St. Louis, Missouri; Institute for Computational Engineering and Sciences, Center for Computational Visualization, The University of Texas at Austin, Austin, Texas; Department of Chemistry and Biochemistry, Center for Theoretical Biological Physics, University of California at San Diego, La Jolla, California; Department of Computer Sciences and Institute for Computational Engineering and Sciences, Center for Computational Visualization, The University of Texas at Austin, Austin, Texas; and ^¶ Department of Chemistry and Biochemistry, Center for Theoretical Biological Physics, Department of Pharmacology, Howard Hughes Medical Institute, University of California at San Diego, La Jolla, California

Correspondence: Address reprint requests to Nathan A. Baker, Dept. of Biochemistry and Molecular Biophysics, Center for Computational Biology, Washington University in St. Louis, 700 S. Euclid Ave., Campus Box 8036, St. Louis, MO 63110. Tel.: 314-362-2040; Fax: 314-362-0234; E-mail: baker{at}biochem.wustl.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORY AND ALGORITHMS VALIDATION OF THE SMOL... APPLICATION OF THE SMOL... CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
This article describes the development and implementation of algorithms to study diffusion in biomolecular systems using continuum mechanics equations. Specifically, finite element methods have been developed to solve the steady-state Smoluchowski equation to calculate ligand binding rate constants for large biomolecules. The resulting software has been validated and applied to mouse acetylcholinesterase. Rates for inhibitor binding to mAChE were calculated at various ionic strengths with several different reaction criteria. The calculated rates were compared with experimental data and show very good agreement when the correct reaction criterion is used. Additionally, these finite element methods require significantly less computational resources than existing particle-based Brownian dynamics methods.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY AND ALGORITHMS VALIDATION OF THE SMOL... APPLICATION OF THE SMOL... CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Diffusion plays a central role in numerous biological processes, governing the kinetic properties of events across a variety of length scales: from ligand binding (Antosiewicz et al., 1995, 1996b; Antosiewicz and McCammon, 1995; Lesyng and McCammon, 1993; McCammon and Karplus, 1977; Northrup et al., 1984; Tan et al., 1993; Tara et al., 1998; Wade et al., 1994) to protein-protein encounter (Elcock et al., 2001; Gabdoulline and Wade, 2001; Northrup and Erickson, 1992; Sheinerman et al., 2000; Zhou, 1997) to signal transmission at synaptic junctions (Franks et al., 2002; Kara and Friedlander, 1998; Roberts, 1994; Smart and McCammon, 1998; Tai et al., 2003; Zoli and Agnati, 1996). Biological simulations have been used to study such diffusion-controlled processes in a number of settings and have provided useful insight into the molecular determinants of the kinetic parameters. However, accurate modeling of diffusion within biomolecular systems while incorporating the effects of ionic strength, solvent, and protein charges, and applying to large biological systems with complex geometries, has proven to be the rate-limiting step for a variety of such simulations.

Currently, standard techniques for modeling diffusional processes can be loosely grouped into particle-based and continuum methods. Particle-based methods are typically stochastic in nature and include Monte Carlo (Berry, 2002; Genest, 1989; Saxton, 1992; Stiles and Bartol, 2000; Brownian dynamics (BD) (McCammon, 1987; Northrup et al., 1988a; Wade et al., 1993), and Langevin dynamics (Eastman and Doniach, 1998; Yeomans-Reyna and Medina-Noyola, 2001) simulations. The connection between BD simulations and of the calculation of association rate constants was established by Northrup, Allison, and McCammon (Northrup et al., 1984) and has been studied by numerous others (Antosiewicz et al., 1996a; Antosiewicz and McCammon, 1995; Chung et al., 2002; Northrup et al., 1988b; Tan et al., 1993; Tara et al., 1998; Wade et al., 1993; Zhou, 1993; Zhou et al., 1998a; Zhou and Szabo, 1996; Zou et al., 2000). In contrast to particle-based approaches, continuum methods describe diffusional processes in terms of probability or concentration profiles rather than simulating the stochastic motion of individual particles. Continuum methods are typically based on solutions of partial differential equations such as the diffusion or Smoluchowski equation (Chan and Halle, 1984; Gardiner, 1997; Lenzi et al., 2003; Smart and McCammon, 1998; Tai et al., 2003); these solutions can then be processed to determine ligand-protein binding (Agmon et al., 1991; Smart and McCammon, 1998; Tai et al., 2003; Zhou, 1990) or dissociation (Agmon, 1984). These methods have been particularly popular in the fields of ion channel (Coalson and Duncan, 1992; Gillespe et al., 2002; Im and Roux, 2002; Kurnikova et al., 1999a) and semiconductor (Selberherr, 1984) modeling.

Both particle-based and continuum diffusion methods have their relative strengths. Particle-based methods can deal with a wide range of diffusing molecular geometries, whereas continuum methods are restricted to spherical ligands. This spherical approximation is likely to be most appropriate for substrates with charge distributions with small multipole moments and reaction criteria that do not require a detailed fit between the substrate and macromolecule. Additionally, particle-based approaches permit the natural inclusion of stochastic reaction phenomena and other complicated boundary conditions. However, the stochastic nature of particle-based approaches can lead to convergence problems which are not present in the deterministic continuum method. Furthermore, continuum approaches facilitate the inclusion of other continuum phenomena such as elastic deformations and fluid flow. Finally, as illustrated in this article, the computational cost of the continuum simulations is significantly smaller than for particle-based methods.

Currently, there are a number of tools available for particle-based biomolecular diffusion simulations, including: SDA (Gabdoulline and Wade, 1997), UHBD (Briggs et al., 1995), MacroDox (Northrup et al., 1999), and MCell (Stiles and Bartol, 2000). However, there are no biomolecule-specific tools available for analyzing diffusion via continuum models and only a few general diffusion tools (Krissinel and Agmon, 1996). The objective of this study is to develop, validate, and apply algorithms to solve the steady-state Smoluchowski equation (SSSE) with finite element methods using realistic biomolecular geometries to determine the steady-state ligand binding rate constant. Specific aims in this study include: development of the adaptive meshing method to realistically describe biomolecular geometries; development of the finite element solver of the steady-state Smoluchowski equation to analyze the concentration of the diffusing particles and calculate the association rate constants; validation of the SSSE with a simple spherical biomolecular system through the comparison with the analytical results; and application of the validated SSSE solver to mouse acetylcholinesterase (mAChE) ligand binding.

	THEORY AND ALGORITHMS

TOP ABSTRACT INTRODUCTION THEORY AND ALGORITHMS VALIDATION OF THE SMOL... APPLICATION OF THE SMOL... CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
The steady-state Smoluchowski equation
The Smoluchowski equation describes the overdamped (i.e., instantaneous momentum relaxation) dynamics of multiple particles while neglecting interparticle interactions (Smoluchowski, 1917; Szabo et al., 1988; Zhou, 1990). For a stationary diffusion process, the Smoluchowski equation has the steady-state form of

(1)

where Lp(x) represents (t is the time), p(x) is the distribution function of the reactants, D(x) is the diffusion coefficient, ß = 1/kT is the inverse Boltzmann energy, k is the Boltzmann constant, T is the temperature, and W(x) is the potential mean force (PMF) for the diffusing particle. The above steady-state Smoluchowski equation (SSSE) can also be written in terms of the flux operator J, which generates vector-valued functions and is defined as

(2)

allowing Eq. 1 to be rewritten as

(3)

The SSSE can be solved to determine bimolecular diffusional encounter rates. Following the work of Zhou (1990), the application of the SSSE to this problem involves the solution of Eq. 3 in a three-dimensional domain , with the following boundary conditions: a bulk Dirichlet condition on the outer boundary ,

(4)

specifying the bulk concentration p_bulk; a reactive Robin or Dirichlet condition on the active site boundary ,

(5)

or

(6)

providing either an intrinsic reaction rate (x) or an absolute reactivity, respectively; and a reflective Neumann condition on the nonreactive boundary

(7)

The problem domain is depicted in Fig. 1; is a simply connected domain with boundary _b, which represents the volume containing the reactive object and the solvent. The domain is a simply connected region representing the reactive object with boundary _ar = _{a r} such that _{a r} = 0. The _r portion of this boundary represents the nonreactive surface of the biomolecular object whereas the _a boundary represents the reactive region of the object. The problem domain is the region in which substrate is allowed to diffuse, namely = - . To ensure the well-posedness of the SSSE Smoluchowski equation, we place several physically based restrictions on the problem data: the potential of mean force (PMF, the effect interaction potential) is bounded and dies off with increasing distance, the diffusion constant is positive definite and finite, and the temperature is positive definite and finite. Finally, the diffusion-influenced biomolecular reaction rate constant is obtained from the flux by integration over the active site boundary, as

(8)

View larger version (37K): [in this window] [in a new window]   FIGURE 1  Schematic of problem domain denoting the various surfaces and volumes described in the text.

(9)

where (y_i, r_i) are the coordinates and radii of the N atoms in the biomolecule and is the radius of the diffusing species. This accessibility function provides an abstraction of the biomolecular structure which can then be used as input for advanced volume isocontouring methods. Specifically, (x) provides a (grid-based) dataset which is then isocontoured (at a value of 0.5) via dual contouring methods into a three-dimensional tetrahedral mesh (Zhang et al., 2003b).

The goal of the dual contouring methods used in this work is to tetrahedralize the interval volume between the biomolecular surface and an outer boundary sphere S₁ which is usually 40 times that of the biomolecule (compare to Figs. 2 and 3). The resulting tetrahedral mesh is spatially adaptive and preserves molecular surface features while minimizing the number of simplices. The four main steps of our adaptive tetrahedral meshing of the problem domain are described in the following sections.

View larger version (68K): [in this window] [in a new window]   FIGURE 2  Adaptive tetrahedral meshes for mouse acetylcholinesterase. (Bottom) The molecular surface and outer sphere S0; the active site gorge is shown in greater detail inside the red box. (Middle) Magnification of the red box in the bottom picture. (Top) The tetrahedral mesh of the interval volume between the molecular surface and the outer sphere S0 (cross section).

View larger version (34K): [in this window] [in a new window]   FIGURE 3  Data scaling for the adaptive mesh; surfaces are described in the text.

Tetrahedral mesh extraction
Dual contouring (Ju et al., 2003) uses an octree data structure to analyze mesh edges, called sign-change edges, which have endpoints lying on different sides of the isosurface of interest. The mesh adaptivity is determined during a top-down octree construction. Each sign change edge is shared by either three (adaptive case) or four (uniform case) octree cells. A minimizer point is calculated for each cell by minimizing a predefined quadratic error metric (Garland and Heckbert, 1998). Finally, for each sign change edge, a quadrilateral or triangle is constructed by connecting the minimizer points; these resulting quadrilaterals and triangles provide a dual piecewise linear approximation of the isosurface. This dual contouring method has already been extended to extract tetrahedral meshes from volumetric scalar fields (Zhang et al., 2003b) following a similar procedure to the two-dimensional case.

Three-dimensional dual contouring is used to construct the tetrahedral mesh between the biomolecular surface and the sphere S₀. When used as described above, the resulting mesh is finest around the molecular surface, and gradually gets coarser away from the molecule, toward the boundary S₀. However, the adaptive nature of this method can also be exploited to provide additional detail in biologically important regions of the molecule such as active sites. As illustrated in Figs. 2 and 6 for mouse acetylcholinesterase, we control the mesh resolution to generate a very fine mesh in the active site gorge region and a relatively coarse mesh elsewhere near the biomolecular surface.

View larger version (38K): [in this window] [in a new window]   FIGURE 6  Reactive surfaces 1–6 of mAChE (bottom to top) from the reactive boundary definition described in the text. The views of these surfaces start from the outside of the protein (bottom) and move to the interior (top) as the surfaces move inside the active site gorge.

Quality detection and improvement
The tetrahedral mesh generated by the dual contouring method may have isolated vertices, nonsimple components, or overlapping tetrahedra—problems which must be removed before the mesh can be used to solve the SSSE. Isolated vertices and nonsimple components are identified and removed using the connectivity information (i.e., edges and simplices) intrinsic to the finite element mesh. Overlapping simplices are corrected using methods previously described by Zhang et al. (2003b).

In addition to correcting these topological problems, the geometric quality of the mesh must be checked and, if necessary, improved. Poor quality tetrahedra called "slivers" are common in most tetrahedral mesh generation methods. The presence of such slivers can often confound finite element solvers and therefore these poor-quality tetrahedra must be removed (Cheng and Dey, 2002; Cheng et al., 2000) before using the generated meshes to solve the SSSE. We combine edge contraction and smoothing methods to improve the quality of the meshes based on tetrahedral-edge ratios, Joe-Liu shape metrics (Liu and Joe, 1994), and minimum volume bounds. Edge contraction (Cheng and Dey, 2002; Cheng et al., 2000) improves tetrahedral shape by merging vertices to combine smaller simplices into larger, better quality tetrahedra. Additionally, we use a smoothing method based on multilinear averaging (Zhang et al., 2003b) to improve the Joe-Liu shape metric and the minimum volume bounds.

Finite element discretization of the Smoluchowski equation
To solve the SSSE numerically as a finite system, it is necessary to truncate and discretize the infinitely large problem domain implicit in Eqs. 1–7. We solve Eq. 1 using finite element methods (Axelsson and Barker, 2001; Braess, 1997) inside a domain that is 40 times the scale of the size of the biomolecule. Since the effects of the PMF are not included beyond the outer boundary, this large size is typically necessary for electrostatic forces to decay to 0 and/or avoid more complicated outer boundary conditions. However, due to the adaptive nature of the finite element meshes used to discretize the system and our multiresolution approach to the calculation of the potential mean force (see below and Baker et al., 2001), this large mesh does not impose a significant computational burden. The outer boundary of the diffusion domain is subject to the Dirichlet boundary condition per Eq. 4. The inner (molecular) boundary is assigned reactive and reflective conditions as shown in Fig. 1 and described in Eqs. 5–7. The whole domain is discretized on an adaptively generated tetrahedral mesh using the methods outlined above.

The resulting tetrahedral mesh forms the structure over which we define a function space V_h = span{v_i}, where {v_i} is the set of piecewise-linear finite element basis functions defined over each tetrahedral vertex. The solution to the SSSE is approximated by a function, constructed from the linear combination of basis functions,

(10)

The trace function is not explicitly constructed, but is assumed to satisfy the Dirichlet boundary conditions in Eqs. 4 and 6. For this construction of p_h from piecewise-linear functions to be successful, we must restate the SSSE equations in their weak form. Clearly, piecewise-linear functions do not provide a well-defined second derivative, as required by Eq. 1. This difficulty can be overcome by integrating the SSSE with a test function v(x), as

(11)

Integrating the above equation by parts via Green's theorem gives

(12)

and incorporation of the various boundaries conditions Eqs. 4–7, allows the boundary integral to be simplified as

(13a)

or

(13b)

depending on the choice of reactive boundary condition (compare to Eq. 5 and Eq. 6, respectively). Finally, the above integrals can be rewritten in terms of the bilinear form F(p),v,

(14a)

(14b)

to give the weak form of the SSSE,

(15)

This form of the SSSE requires only one order of differentiation under an integral and is therefore a weaker formulation of the SMOL equation than the original second-order differential Eq. 1. Equation 15 is the foundation for the finite element solvers used in this work.

In the Smoluchowski equation, the diffusion coefficient D(x) and electrostatic potential mean force W(x) all are spatially dependent. Ignoring hydrodynamic interactions (which could be included as part of the spatial dependence of D), the diffusion coefficient is treated as a constant. In most BD models, the electrostatic potential is treated independently from the diffusing species—eliminating the possibility of screening by ligand, substrate inhibition, etc. This is also the model we will employ for these initial studies, implying that Eq. 15 is a linear equation. Discretization of Eq. 15 with a finite element basis (per Eq. 10) leads to a linear system of (sparse) equations which can be solved using standard linear algebra methods. However, it is important to note that this model could easily be extended to the related Poisson-Nernst-Planck (PNP) system where the PMF is coupled to the diffusing species. For PNP equations, well-established nonlinear methods will need to be employed (Holst and Saied, 1995; Koumanov et al., 2003; Kurnikova et al., 1999b; Schuss et al., 2001). Additionally, like many PNP methods, the present model does not consider correlations between diffusing species. Such effects have been shown to be important in confined spaces such as ion channels (Boda et al., 2002; Gillespe et al., 2002) or at higher concentrations of multivalent ions (Holm et al., 2001). For the low concentrations of ligand used in many protein-ligand diffusive encounter simulations, these correlation effects are likely to be relatively small; however, it is important that such effects are considered when performing diffusion simulations. Finally, the present model does not include dielectric boundary or apolar forces which have been shown to be important factors for diffusion behavior in some biological systems (Nadler et al., 2003). Such effects are very important for diffusion of charged species in confined spaces (e.g., ion channels) but have a much less significant impact for ligand-protein binding at relatively low concentrations (where dielectric boundary and apolar forces often cancel).

Implementation
We have developed a software package called SMOL which uses the finite element method to solve the SSSE for biomolecular systems. The SMOL program is based on the LBIE-Mesh software (http://www.ices.utexas.edu/CCV/software), FEtk finite element software library (M. Holst, Adaptive Multilevel Finite Element Methods on Manifolds and their Implementation in FEtk; in preparation; currently available as a technical report and User's guide to the FEtk software; see http://www.fetk.org/) and the APBS software package (http://agave.wustl.edu/apbs/, Baker et al., 2001). Specifically, SMOL uses LBIE-mesh for construction of the initial mesh, FEtk's finite element infrastructure for the discretization and solution of the linear system of equations. APBS is used to calculate the electrostatic component of the potential mean force by solution of the Poisson-Boltzmann equation. The SMOL software is still actively under development and will be released under the open-source GNU public license (see http://www.gnu.org/). Adventurous readers are welcome to contact the authors for access to early versions of the software; others should visit http://agave.wustl.edu/ for more information about user-friendly public releases.

The sequence of operations followed by the SMOL program during a typical solution of the SSSE for biomolecular rate calculation is outlined in the following sections.

Parameterization
The biomolecular coordinates from the Protein Data Bank must be parameterized with appropriate atomic radii and charges. The current work uses the CHARMM22 force field (Brooks et al., 1983); conversion tools for several other force fields are provided with the APBS software or available through the PDB2PQR web service (http://nbcr.sdsc.edu/pdb2pqr/).

Mesh construction
First, a grid-based form characteristic accessibility function (x) is generated with APBS using the diffusing species' radius (compare to Eq. 9) and the parameterized biomolecular structure as input data. This discretized characteristic function is then used with the adaptive meshing tools described above to generate quality tetrahedral meshes of the problem domain. Finally, boundary conditions as in Eqs. 4–7 are assigned to the mesh based on user-defined spatial criteria using knowledge of the active site location and other structural features.

Potential of mean force
Only electrostatic contributions to the PMF are included in this version of the software. The electrostatic potential is calculated over the entire diffusion domain of the biomolecule with APBS. This data is stored in a parallel, multiresolution format where (sometimes multiple) high-resolution meshes provide PMF values near the molecular surface while coarser meshes provide values away from the biomolecule. This decomposition of the potential closely follows the parallel-focusing methods described elsewhere (Baker et al., 2001).

Solution of the SSSE and analysis
Given the above PMF data, finite element mesh, and user-defined values for the diffusing particle's charge, concentration, and diffusion constant data, the SSSE is solved with the finite element methods described above. After the successful solution of the differential equation, output is provided for the steady-state rate constant and (if requested) the concentration profile of ligand around the biomolecule.

	VALIDATION OF THE SMOL PROGRAM WITH A SPHERICAL TEST CASE

TOP ABSTRACT INTRODUCTION THEORY AND ALGORITHMS VALIDATION OF THE SMOL... APPLICATION OF THE SMOL... CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Before applying the SMOL program to a biomolecular system with complex geometry, we first tested it with the classic spherical system (Krissinel and Agmon, 1996) and compared the calculated result with the known analytical solution. For this test case, we chose a fixed sphere with an 8 Å radius and +1 e charge and a diffusing sphere with a 2 Å radius and variable charge. The partial domain for this spherical molecule is shown in Fig. 4 A. The entire problem domain was discretized with 1,024,752 tetrahedral elements. A detailed view of the surface mesh for the stationary sphere is also shown in Fig. 4 B. While solving the Smoluchowski equation, the whole surface of biomolecule sphere _a was treated with a perfectly absorbing zero Dirichlet reactive boundary condition as Eq. 6. The diffusing particle's dimensionless bulk concentration was set to 1. Ignoring hydrodynamic interactions, the diffusion constant D is calculated as 7.8 x 10⁴ Ų/µs using the Stokes-Einstein equation with a hydrodynamic radius of 3.5 Å, solvent viscosity of 0.891 x 10^-3 kg/(m s), and 298 K temperature.

View larger version (44K): [in this window] [in a new window]   FIGURE 4  Illustration of the discretized problem domain for the spherical test case. (a) Partial domains of the fixed sphere; the outer boundary of the domain is 40 times the radius of the sphere. (b) Subset of the mesh near the fixed sphere surface.

(16)

with boundary conditions

(17)

where r₁ and r₂ are the radii for the reactive and outer boundaries. The analytical expression for the association rate constant obtained from the solution of Eqs. 16 and 17 is

(18)

This analytical form of the solution was evaluated using the parameters listed above with varying ligand charges in the infinite dilution limit: r₂ . Table 1 presents the numerical values of these analytical reaction rates as a function of effective ligand charge.

	APPLICATION OF THE SMOL PROGRAM TO MOUSE ACETYLCHOLINESTERASE

TOP ABSTRACT INTRODUCTION THEORY AND ALGORITHMS VALIDATION OF THE SMOL... APPLICATION OF THE SMOL... CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
One of the major advantages of continuum methods such as the SSSE is the ability to simulate diffusion to large biological systems with complex geometries with significantly lower computational cost than Brownian dynamics techniques. This section demonstrates the use of SSSE to study the ligand binding kinetics of acetylcholinesterase (AChE, E.C. 3.1.1.7) (Quinn et al., 1995). AChE is a serine esterase that terminates the activity of acetylcholine (ACh) within the cholinergic synapse by hydrolysis of the ACh ester bond to produce acetate and choline (Berg et al., 1995). Hydrolysis of ACh occurs in the active site of AChE, which lies at the base of a 20 Å-deep gorge within the enzyme. The rate-limiting step of ACh hydrolysis by AChE is the diffusional encounter (Bazelyansky et al., 1986; Berman et al., 1991; Nolte et al., 1980), making the system a popular target for both experimental (Bourne et al., 1995; Radic et al., 1997; Velsor et al., 2003) and computational diffusion studies (Tan et al., 1993; Tara et al., 1998).

To provide additional data for assessment of SMOL and the SSSE in rate constant calculations, we based our analysis of AChE kinetics on previous BD work by Tara et al. (1998). Specifically, we used a mouse AChE (mAChE) structure adopted from the crystal structure of the mAChE-fasciculin 2 complex (1MAH) (Bourne et al., 1995) and perturbed by Tara and co-workers via molecular dynamics simulations with an ACh-like ligand in the active site gorge (Tara et al., 1998) to produce gorge conformations with wider widths than the original x-ray structure. The diffusing ligand was modeled as a sphere with an exclusion radius of 2.0 Å and a diffusion constant of 7.8 x 10⁴ Ų/µs. This perturbation was necessary for computational diffusion simulations with a fixed biomolecular structure. Related work (Baker and McCammon, 1999; Zhou et al., 1998b) supports the use of this open gorge state due to its relevance in conformational gating dynamics of the enzyme. The geometry of mAChE and its active site gorge are shown in Figs. 5 and 6.

View larger version (77K): [in this window] [in a new window]   FIGURE 5  Geometry of mAChE finite element mesh (left) and its active site gorge (right).

Variables and parameters for the biomolecule and diffusion domain
For the purpose of the PMF (electrostatics) calculation for mAChE, partial charges and radii were assigned from the CHARMM22 force field. Poisson-Boltzmann electrostatics calculations were performed using dielectric values of 4 and 78 for the protein and solvent, respectively; a solvent probe radius of 1.4 Å; and an ion exclusion layer of 2.0 Å. APBS was used to solve the nonlinear Poisson-Boltzmann equation using the parallel focusing method (Baker et al., 2001). Sets of nested potential grids were obtained with the finest grid having dimensions 76 x 66 x 91 Å with 161 grid points in each direction. For the study of the effect of ionic strength on reaction rate, separate calculations were performed at ionic strengths of 0.00, 0.050, 0.100, 0.150, 0.200, 0.250, 0.300, 0.450, 0.600, and 0.670 M.

The problem domain definition for the mAChE system is shown in Fig. 1. The outer boundary _b has a radius of 40 times the size of mAChE molecule, and is given a bulk concentration of unity. Reactive surface _a is assigned zero (perfectly absorbing) Dirichlet boundary conditions (compare to Eq. 6) and the remainder of the biomolecular surface is treated as the reflective boundary _r with Neumann condition (Eq. 7). The adaptive meshing method described above was used to discretize the problem domain. To accurately reproduce the topology of the 20 Å-deep active site gorge, a finer mesh was used for discretization near catalytic site and gorge; other areas of the molecule are discretized with relatively coarser mesh. The resulting mesh contained 656,823 simplices and 121,670 vertices. A subset of the domain is shown in Figs. 3, 5, and 6 to illustrate the general features of the mesh. Detailed views of the mesh near the active site gorge are shown in Figs. 5 and 6.

The diffusing particle was treated as a sphere with a +1 e charge, a 2.0 Å exclusion radius, and a diffusion constant of 7.8 x 10⁴ Ų/µs (parameters were taken from Tara et al., 1998, which are similar to those for the TFK⁺ ligand).

Calculating the reaction rate for mAChE-TFK⁺ encounter with SMOL
Reaction rates for mAChE-TFK⁺ diffusional encounter were calculated from the SSSE via the SMOL program for ionic strengths ranging from 0.000 to 0.670 M NaCl and for the reactive boundary definitions described above. All calculations were performed on Intel Pentium 4 CPU 2.20-GHz computer systems with 1.5 GB RAM running RedHat Linux 9. SMOL executables were compiled with the Intel FORTRAN and C compilers using the -O2 optimization. The SMOL calculations took an average of 55 s per reactive surface. The resulting reaction rates are presented in Fig. 7 and Table 2; analysis of this data and comparison with both experimental and BD results are presented below.

View larger version (12K): [in this window] [in a new window]   FIGURE 7  Reaction rates of mAChE calculated with the SMOL FE program with reactive surface 1 (•), calculated from Brownian dynamics with reactive surface 1 (), and from experimental data (Radic et al. 1997) fit to the Debye-Hückel limiting law (solid line).

To compare the results from SSSE and BD calculations for the mAChE system (Tara et al., 1998), the binding rate constants for surfaces 1–6 and ionic strengths (0–670 mM) were plotted in Fig. 8. This comparison shows that most of the SSSE calculations give smaller reaction rates than BD simulations and therefore lead to the disagreement between the methods in the choice of the best reactive surface. There are multiple reasons for this disagreement. Although both the particle-based and continuum methods have the same underlying physics, they employ very different implementations. First and foremost, the particle-based approach uses discrete particle displacements and therefore should only give perfect agreement with the continuum simulations in the limit of infinitesimal step sizes. Second, BD simulations enforce reflective boundary conditions (compare to Eq. 7) only approximately by simply rejecting any steps which take particles into nonreactive surfaces. Given the obvious role of the reflective biomolecular surface in the gorge, it is not surprising that there are significant differences between BD and SSSE results. Finally, we have not implemented any of the error-based adaptive refinement methods available for finite element methods (Axelsson and Barker, 1984; Braess, 1997).

View larger version (16K): [in this window] [in a new window]   FIGURE 8  Correlation plots for calculations of reaction rate using the SMOL program and using BD with spherical surface reactive boundary definitions. For each method, calculations are made for six reactive surfaces and ionic strength ranges 0.00, 0.05, 0.10, 0.15, 0.3, 0.45, 0.6, and 0.67 M. For • data points, BD results are plotted on the x axis, SMOL results are plotted on the y axis. The solid line is the linear fit of the data, and its slope is 0.84 with a Pearson correlation coefficient of 0.91.

(19)

where k_on is the observed binding rate constant, is the effective infinite ionic strength limiting rate, k_on⁰ is the effective 0 ionic strength rate, z_E is the effective enzyme charge, z_I is the effective inhibitor charge, and I is the ionic strength. Following the analysis of Radic et al. (1997) of their experimental data, nonlinear regression analysis was used to fit the calculated results at concentrations between 0.050 and 0.450 M and obtain values for k_on⁰ and z_E. The value of z_I was held fixed at +1 e and k_on^H was set to the value of k_on calculated at 0.670 M ionic strength. The results of the rate calculated and the fit parameters are given in Table 2. The data show that the calculated rates fit Eq. 19 very well over the ionic strength range for which the Debye-Hückel limiting law is valid (0.05 M to 0.5 M).

Table 2 shows that the rate constants for reactive surfaces 3–6 do not decrease monotonically with the decreasing reactive patch size. This interesting behavior appears to be influenced by the biomolecular electrostatics; the expected monotonic decrease is observed when the ligand charge is set to 0. A number of residues influence the movement of the ligand within the gorge. Previous work (Tara et al., 1998) showed that residue Asp-74, immediately below reactive surface 3 and above surface 4, plays an important role for the ligand kinetics within the gorge and is a likely candidate for the source of this behavior.

When performing computational diffusion simulations, both particle-based and continuum methods need to choose the appropriate reaction criterion for calculation binding rate constants. This choice is usually made by comparing the simulation results with available experimental data. Fig. 7 indicates that SSSE-derived reaction rates show the best overall agreement with experiment for surface 1. Table 2 also shows that, whereas the calculated results consistently underestimate k_on⁰, the best agreement is observed for surface 1. The results obtained at surfaces 3–6 are somewhat less than the experimental values but do show good agreement with one another, indicating that the value of k_on⁰ is insensitive to the choice of surface definition in this region. Likewise, surfaces 3–6 all give similar values for the effective enzyme charge z_E; however, unlike k_on⁰, the values of z_E from surfaces 3–6 all agree well with experiment. As for k_on⁰, reactive surface 1 gives a calculated value, z_E, which matches experimental observations.

The above results show that the SSSE gives very good agreement with experiment when the correct reaction criterion is used. In the case of mAChE, this criterion corresponds to the reactive boundary located at surface 1. The BD results by repeating the simulation of Tara et al. (1998) also used spherical reaction criteria similar to those used in the present study. However, the reaction rates results suggest that surface 3 provides the best agreement with experiment. The differences between BD and SSSE results were discussed above (see also Fig. 8). It is not surprising that the two methods give different rates at the same surfaces; it is encouraging that both methods offer reaction conditions which give good agreement with the experimental results.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION THEORY AND ALGORITHMS VALIDATION OF THE SMOL... APPLICATION OF THE SMOL... CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
In this study, we describe new continuum-based methods for studying diffusion in biomolecular systems. Specifically, we present the SMOL software package, a finite element-based software package for solving the SSSE to calculate ligand binding rate constants for large biomolecules. Additionally, we describe new adaptive meshing methods developed to discretize biomolecular systems into finite element meshes which respect the geometry of the biomolecule. Although not presented in this study, it is important to note that the new meshing methods could be useful in a variety of biological simulations including computational studies of biomolecular electrophoresis (Allison, 2001), elasticity (Zhang et al., 2003a), and electrostatics (Baker et al., 2000; Cortis and Friesner, 1997a,b; Holst et al., 2000). The SMOL software was validated using a spherical analytical test case. Rates for inhibitor binding to mAChE were calculated at various ionic strengths with several different reaction criteria. The calculated rates were compared with experimental data (Radic et al., 1997) and show very good agreement with experiment while requiring substantially less computational effort than existing particle-based Brownian dynamics methods.

Since the inception of this work, new Brownian dynamics methods have been developed which dramatically improve the speed and convergence of these particle-based simulations (Zou and Skeel, 2003). However, one of the key ingredients to these new methods is the construction of an a priori approximate probability distribution to bias and accelerate the BD calculations. One interesting future direction of the present research would be the use of these continuum methods to develop approximate biasing functions for the new biased BD methodology of Zou and Skeel.

This initial research lays the groundwork for several new directions. Of particular interest is the integration of molecular-scale information into simulations of cellular-scale systems such as the neuromuscular junction (Smart and McCammon, 1998; Tai et al., 2003). Additionally, this new finite element framework should facilitate the incorporation of other continuum mechanics phenomena into biomolecular simulations. The ultimate goal of this work is to develop scalable methods and theories that will allow researchers to begin to study biological macromolecules in a cellular context.

	ACKNOWLEDGEMENTS

TOP ABSTRACT INTRODUCTION THEORY AND ALGORITHMS VALIDATION OF THE SMOL... APPLICATION OF THE SMOL... CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
N.A.B. and Y.H.S. thank Steve Bond for helpful discussions, Mike Holst for advice and access to the FEtk software, and the National Partnership for Advanced Computational Infrastructure and the National Biomedical Computation Resource for financial and computational support. N.A.B. is an Alfred P. Sloan Research Fellow. Y.H.S. also thanks Todd Dolinsky for help during the development of the SMOL software.

Work at University of California at San Diego is supported in part by the National Institutes of Health, the National Science Foundation, the Howard Hughes Medical Institute, the National Science Foundation Center for Theoretical Biological Physics, the National Biomedical Computation Resource, the W.M. Keck Foundation, and Accelrys, Inc. The work at the University of Texas was supported in part by National Science Foundation grants ACI-9982297, 0220037, and CCR-9988357, and a subcontract from the University of California at San Diego as part of the National Science Foundation-National Partnership for Advanced Computational Infrastructure.

Submitted on August 25, 2003; accepted for publication October 14, 2003.

	REFERENCES

TOP ABSTRACT INTRODUCTION THEORY AND ALGORITHMS VALIDATION OF THE SMOL... APPLICATION OF THE SMOL... CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Agmon, N. 1984. Unimolecular dissociation as diffusion with a radiation boundary condition. J. Chem. Phys. 80:5049–5054.

Agmon, N., H. Schnorer, and A. Blumen. 1991. Competitive reversible binding: a bimolecular boundary condition for the diffusion equation. J. Phys. Chem. 95:7326–7330.

Allison, S. A. 2001. Boundary element modeling of biomolecular transport. Biophys. Chem. 93:197–213.[Medline]

Antosiewicz, J., J. M. Briggs, and J. A. McCammon. 1996a. Orientational steering in enzyme-substrate association: ionic strength dependence of hydrodynamic torque effects. Eur. Biophys. J. 24:137–141.[Medline]

Antosiewicz, J., S. T. Wlodek, and J. A. McCammon. 1996b. Acetylcholinesterase: role of the enzyme's charge distribution in steering charged ligands toward the active site. Biopolymers. 39:85–94.[Medline]

Antosiewicz, J., M. K. Gilson, I. H. Lee, and J. A. McCammon. 1995. Acetylcholinesterase: diffusional encounter rate constants for dumbbell models of ligand. Biophys. J. 68:62–68.[Abstract/Free Full Text]

Antosiewicz, J., and J. A. McCammon. 1995. Electrostatic and hydrodynamic orientational steering effects in enzyme-substrate association. Biophys. J. 69:57–65.[Abstract/Free Full Text]

Axelsson, O., and V. A. Barker. 2001. Finite Element Solution of Boundary Value Problems. Theory and Computation. Society for Industrial and Applied Mathematics, Philadelphia, PA.

Axelsson, O., and V. A. Barker. 1984. Finite Element Solution of Boundary Value Problems. Theory and Computation. Academic Press, San Diego, CA.

Baker, N. A., and J. A. McCammon. 1999. Non-Boltzmann rate distributions in stochastically gated reactions. J. Phys. Chem. B. 103:615–617.

Baker, N., M. Holst, and F. Wang. 2000. Adaptive multilevel finite element solution of the Poisson-Boltzmann equation II: refinement at solvent accessible surfaces in biomolecular systems. J. Comput. Chem. 21:1342–1352.

Baker, N. A., D. Sept, M. J. Holst, and J. A. McCammon. 2001. The adaptive multilevel finite element solution of the Poisson-Boltzmann equation on massively parallel computers. IBM J. Res. Dev. 45:427–438.

Bazelyansky, M., E. Robey, and J. F. Kirsch. 1986. Fractional diffusion-limited component of reactions catalyzed by acetylcholinesterase. Biochemistry. 25:125–130.[Medline]

Berg, J. M., J. L. Tymoczko, and L. Stryer. 1995. Biochemistry. W. H. Freeman & Co., New York, NY.

Berman, H. A., K. Leonard, and M. W. Nowak. 1991. Cholinesterases: Structure, Function, Mechanism, Genetics and Cell Biology. J. Massoulie, editor. American Chemical Society, Washington, DC.

Berry, H. 2002. Monte Carlo simulations of enzyme reactions in two dimensions: fractal kinetics and spatial segregation. Biophys. J. 83:1891–1901.[Abstract/Free Full Text]

Boda, D., D. Busath, B. Eisenberg, D. Henderson, and W. Nonner. 2002. Monte Carlo simulations of ion selectivity in a biological Na⁺ channel: charge-space competition. Phys. Chem. Chem. Phys. 4:5154–5160.

Bourne, Y., P. Taylor, and P. Marchot. 1995. Acetylcholinesterase inhibition by fasciculin: crystal structure of the complex. Cell. 83:503–512.[Medline]

Braess, D. 1997. Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics. Cambridge University Press, New York.

Briggs, J., J. Madura, M. Davis, M. Gilson, J. Antosiewicz, B. Luty, R. Wade, B. Bagheri, A. Ilin, R. Tan, and J. A. McCammon. 1995. University of Houston Brownian Dynamics Program User's Guide and Programmer's Manual Release 5.1.

Brooks, B. R., R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and M. Karplus. 1983. A program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 4:187–217.

Chan, D. Y., and B. Halle. 1984. The Smoluchowski-Poisson-Boltzmann description of ion diffusion at charged interfaces. Biophys. J. 46:387–407.[Abstract/Free Full Text]

Cheng, S.-W., and T. K. Dey. 2002. Quality meshing with weighted Delaunay refinement. Proc. 13th Annual Symposium on Discrete Algorithms. ACM Press, San Francisco, CA. 137–146.

Cheng, S.-W., T. K. Dey, H. Edelsbrunner, M. A. Facello, and S. Teng. 2000. Sliver exudation. J. ACM. 47:883–904.

Chung, S. H., T. W. Allen, and S. Kuyucak. 2002. Modeling diverse range of potassium channels with Brownian dynamics. Biophys. J. 83:263–277.[Abstract/Free Full Text]

Coalson, R. D., and A. Duncan. 1992. Systematic ionic screening theory of macroions. J. Chem. Phys. 97:5653–5661.

Cortis, C. M., and R. A. Friesner. 1997a. An automatic three-dimensional finite element mesh generation system for the Poisson-Boltzmann equation. J. Comput. Chem. 18:1570–1590.

Cortis, C. M., and R. A. Friesner. 1997b. Numerical solution of the Poisson-Boltzmann equation using tetrahedral finite-element meshes. J. Comput. Chem. 18:1591–1608.

Eastman, P., and S. Doniach. 1998. Multiple time step diffusive Langevin dynamics for proteins. Proteins. 30:215–227.[Medline]

Elcock, A. H., D. Sept, and J. A. McCammon. 2001. Computer simulation of protein-protein interactions. J. Phys. Chem. B. 105:1504–1518.

Franks, K. M., T. M. Bartol, Jr., and T. J. Sejnowski. 2002. A Monte Carlo model reveals independent signaling at central glutamatergic synapses. Biophys. J. 83:2333–2348.[Abstract/Free Full Text]

Gabdoulline, R. R., and R. C. Wade. 1997. Simulation of the diffusional association of barnase and barstar. Biophys. J. 72:1917–1929.[Abstract/Free Full Text]

Gabdoulline, R. R., and R. C. Wade. 2001. Protein-protein association: investigation of factors influencing association rates by Brownian dynamics simulations. J. Mol. Biol. 306:1139–1155.[Medline]

Gardiner, C. W. 1997. Handbook of Stochastic Methods For Physics, Chemistry and the Natural Sciences. H. Haken, editor. Springer-Verlag, Berlin, Heidelberg, New York.

Garland, M., and P. S. Heckbert. 1998. Simplifying surfaces with color and texture using Quadric error metrics. In Proceedings of the Conference on Visualization '98. D. Ebert, H. Hagen, and H. Rushmeir editors. IEEE Computer Society Press, Research Triangle Park, NC. 263–270.

Genest, D. 1989. A Monte Carlo simulation study of the influence of internal motions on the molecular conformation deduced from two-dimensional NMR experiments. Biopolymers. 28:1903–1911.[Medline]

Gillespe, D., W. Nonner, and R. S. Eisenberg. 2002. Coupling Poisson-Nernst-Planck and density functional theory to calculate ion flux. J. Phys. 14:12129–12145.

Holm, C., P. Kekicheff, and R. Podgornik. editors. 2001. Electrostatic Effects in Soft Matter and Biophysics. Kluwer Academic Publishers, Boston, MA.

Holst, M., N. Baker, and F. Wang. 2000. Adaptive multilevel finite element solution of the Poisson-Boltzmann equation I: algorithms and examples. J. Comput. Chem. 21:1319–1342.

Holst, M., and F. Saied. 1995. Numerical solution of the nonlinear Poisson-Boltzmann equation: developing more robust and efficient methods. J. Comput. Chem. 16:337–364.

Im, W., and B. Roux. 2002. Ion permeation and selectivity of OmpF porin: a theoretical study based on molecular dynamics, Brownian dynamics, and continuum electrodiffusion theory. J. Mol. Biol. 322:851–869.[Medline]

Ju, T., F. Losass, S. Schaefer, and J. Warren. 2003. Dual contouring of hermite data. SIGGRAPH. 2003:339–346.

Kara, P., and M. J. Friedlander. 1998. Dynamic modulation of cerebral cortex synaptic function by nitric oxide. Prog. Brain Res. 118:183–198.[Medline]

Koumanov, A., U. Zachariae, H. Engelhardt, and A. Karshikoff. 2003. Improved 3D continuum calculations of ion flux through membrane channels. Eur. Biophys. J. 32:689–702.[Medline]

Krissinel, E. B., and N. Agmon. 1996. Spherical symmetric diffusion problem. J. Comput. Chem. 17:1085–1098.

Kurnikova, M. G., R. D. Coalson, P. Graf, and A. Nitzan. 1999a. A lattice relaxation algorithm for 3D Poisson-Nernst-Planck theory with application to ion transport through the gramicidin A channel. Biophys. J. 76:646–656.

Kurnikova, M. G., R. D. Coalson, P. Graf, and A. Nitzan. 1999b. A lattice relaxation algorithm for three-dimensional Poisson-Nernst-Planck theory with application to ion transport through the gramicidin A channel. Biophys. J. 76:642–656.[Abstract/Free Full Text]

Lenzi, E. K., R. S. Mendes, and C. Tsallis. 2003. Crossover in diffusion equation: anomalous and normal behaviors. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 67:031104.[Medline]

Lesyng, B., and J. A. McCammon. 1993. Molecular modeling methods. Basic techniques and challenging problems. Pharmacol. Ther. 60:149–167.[Medline]

Liu, A., and B. Joe. 1994. Relationship between tetrahedron shape measures. BIT. 34:268–287.

McCammon, J. A. 1987. Computer-aided molecular design. Science. 238:486–491.[Abstract/Free Full Text]

McCammon, J. A., and M. Karplus. 1977. Internal motions of antibody molecules. Nature. 268:765–766.[Medline]

Nadler, B., U. Hollerbach, and R. S. Eisenberg. 2003. Dielectric boundary force and its crucial role in gramicidin. Phys. Rev. E. 68:021905.

Nolte, H. J., T. L. Rosenberry, and E. Neumann. 1980. Effective charge on acetylcholinesterase active sites determined from the ionic strength dependence of association rate constants with cationic ligands. Biochemistry. 19:3705–3711.[Medline]

Northrup, S. H., and H. P. Erickson. 1992. Kinetics of protein-protein association explained by Brownian dynamics computer simulation. Proc. Natl. Acad. Sci. USA. 89:3338–3342.[Abstract/Free Full Text]

Northrup, S. H., S. A. Allison, and J. A. McCammon. 1984. Brownian dynamics of diffusion-influenced biomolecular reactions. J. Chem. Phys. 80:1517–1524.

Northrup, S. H., J. O. Boles, and J. C. Reynolds. 1988a. Brownian dynamics of cytochrome c and cytochrome c peroxidase association. Science. 241:67–70.[Abstract/Free Full Text]

Northrup, S. H., J. A. Luton, J. O. Boles, and J. C. Reynolds. 1988b. Brownian dynamics simulation of protein association. J. Comput. Aided Mol. Des. 1:291–311.