| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |


* Department of Biochemistry and Molecular Biophysics, Center for Computational Biology, Washington University in St. Louis, St. Louis, Missouri; and
Institute for Computational Engineering and Sciences, Center for Computational Visualization, and
Department of Computer Sciences and Institute for Computational Engineering and Sciences, Center for Computational Visualization, The University of Texas at Austin, Austin, Texas
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 |
|---|
|
|
|---|
| INTRODUCTION |
|---|
|
|
|---|
In this work, we apply adaptive finite element methods using a posteriori error estimation to describe binding of substrate to wild-type and mutant mouse acetylcholinesterases (mAChEs). The AChE system has been a popular research target for both computational model and experimental studies because its hydrolysis of acetylcholine is diffusion-controlled and strongly influenced by electrostatics (Anglister et al., 1995
; Radic et al., 1997
). Previous computational studies of AChE ligand binding used the traditional spherical reactive surface (Radic et al., 1997
; Song et al., 2003
; Tara et al., 1998
). In this work, we introduce a reactive boundary based on the molecular surface, thereby permitting the mapping of "active site" residues directly to the reactive boundary conditions.
Adaptive finite element solution of SE
Original and discretized steady-state SE
A detailed description of the steady-state SE, its application to bimolecular rate constant calculations, and its solution by finite element discretization of the SE were provided in the previous work (Song et al., 2003
). Here, we present a brief review of the SE and the calculation of rate constants from its solutions.
For a stationary diffusion process, the SE has the following (steady-state) form:
![]() | (1) |
is the probability flux, D(x) is the scalar diffusion coefficient,
is the inverse thermal energy, and W(x) is the potential of mean force (PMF). Calculation of the reaction rate involves the solution of the above equation in a three-dimensional domain
with the following boundary conditions. First, we specify the bulk concentration pbulk via a Dirichlet condition on the outer boundary
![]() | (2) |
for either a finite reactivity
(x) via the Robin condition,
![]() | (3) |
![]() | (4) |
via the Neumann condition
![]() | (5) |
Our observable is the diffusion-influenced biomolecular reaction rate constant k, which can be calculated by integration of the flux over the active site boundary:
![]() | (6) |
![]() | (7) |
![]() | (8) |
is a trace function satisfying the Dirichlet boundary conditions, and Vh is the function space spanned by the discrete basis set.
Error estimation and mesh refinement
As demonstrated in the previous work (Song et al., 2003
), Eq. 8 can be used to solve the SE on a given finite element mesh. However, the quality of the resulting approximate solution depends strongly on the underlying finite element discretization of the problem. As their name implies, error estimation methods assess the accuracy of finite element solutions by providing estimates of the difference between the approximate and "true" solutions (Axelsson and Barker, 1984
; Braess, 1997
). These methods are often used in finite element solutions to guide selective refinement of the finite element mesh (Baker et al., 2000
, 2001a
; Holst et al., 2000
; Holst, 2001
; M. Holst and D. Bernstein, unpublished) and thereby adaptively improve the quality of the numerical solution.
In this article, adaptive mesh refinement methods are used in a "solve-estimate-refine" algorithm as described by Holst and co-workers (Baker et al., 2000
, 2001a
; Holst et al., 2000
; Holst, 2001
; M. Holst and D. Bernstein, unpublished) and implemented in the FEtk software (http://www.fetk.org/). The first step of this procedure (solve) is the calculation of an approximate solution to the SE ph (x) on the current finite element mesh (Song et al., 2003
). In the second step (estimate), this solution is used to provide a per-simplex residual-based a posteriori error estimate
s of the form (Holst, 2001
):
![]() | (9) |
is the current numerical estimate of the flux, f
s denotes a face of simplex, hf is the size of the face f,
denotes the jump across the face of some function v,
is the component of the flux normal to element face f, and the Lesbegue norms are defined as
![]() | (10) |
Finally, in the third step (refine), this per-simplex error estimate
s is used to identify simplices of the finite element mesh where the error is above a particular tolerance. Simplices with a high error estimate value are refined by longest-edge bisection. This entire "solve-estimate-refine" cycle is repeated until the global error
is reduced to an acceptable user-defined level.
As described previously, methods for solving the SE have been implemented in a software package called "SMOL." This software uses the Holst group FEtk toolkit (http://www.fetk.org/) for finite element geometric routines, multilevel solvers, and the residual-based error estimation protocol outlined above (Holst, 2001
).
Validation of the adaptive SMOL finite element solver with a spherical system
To validate the new adaptive finite element features of the SMOL software, we examined a classic spherical test case (Krissinel and Agmon, 1996
) and compared calculated rate constants with the analytical results. For this test case, we chose a fixed sphere with an 8 Å radius and a diffusing sphere with a 2 Å radius; both spheres had variable charge. The smaller sphere's diffusion constant D(x) was chosen as a constant 7.8 x 104 Å2 µs, a value obtained from the Stokes-Einstein relationship for a 23 Å substrate (tetramethyl ammonium) in water (Tara et al., 1998
). The PMF (W(x) in Eq. 1) for these calculations was obtained from Coulomb's law for a homogeneous dielectric of 78.54.
The outer boundary of the diffusion domain was chosen to be 40 times the combined size of the fixed and diffusing spheres (i.e., a 400 Å radius) and the inner boundary (at 10 Å) was uniformly reactive with the perfectly-absorbing boundary condition described by Eq. 4. The entire domain was initially discretized into 445,488 tetrahedral elements using the contouring methods of Zhang and co-workers (Song et al., 2003
; Zhang et al., 2003
). Note that this is a much coarser mesh than used in the previous work (Song et al., 2003
). However, the surface area of the spherical inner boundary of this coarser mesh (1257.98 Å2) differed by only 0.11% from the actual area of a 10 Å sphere (1256.64 Å2), indicating that the finite element mesh realistically represents the boundary geometry.
Fig. 1 presents the binding rate constant calculations as a function of ligand charge for fixed spheres with +1 e and +10 e charges. Results from +1 e fixed-sphere charge show that the adaptive methods generate results which are in much better agreement (<1% relative error) with the analytical solutions than nonadapted results (8% relative error). Furthermore, for the +10 e fixed sphere charge, the adaptive methods can provide up to
25% improvement in the rate constants compared to the nonadapted calculations.
|
30 increase in computation time represents an extreme case of adaptive refinementthe current case was specifically chosen to demonstrate the ability of the method to refine from a very coarse initial mesh to the correct answer. As described in the Conclusions of this work, the "real world" application of the adaptive method would likely start from a much finer mesh and therefore require substantially fewer rounds of adaptive refinement and lower overall computation time.
Rate constant calculations for mAChE ligand binding
AChE (E. C. 3.1.1.1.7) is a serine esterase which hydrolyzes the neurotransmitter acetylcholine at diffusion-limited rates (Anglister et al., 1995
). Previously, we investigated the ligand-binding kinetics of mAChE by calculating steady-state rate constants using a nonadaptive version of the SMOL software and the spherical reactive boundary condition commonly used in Brownian dynamics reaction rate calculations (Song et al., 2003
). Here, we use our new adaptive finite element scheme to investigate the binding kinetics of mutant and wild-type mAChE using a new reactive boundary based on the biomolecular surface.
mAChE domain geometry
Like the previous calculations (Song et al., 2003
), the starting geometry for these calculations is the "open" mAChE structure used by Tara et al. (1998)
to study mAChE binding kinetics. Using an outer boundary 40 times the size of the biomolecule (an ellipsoid with dimensions 3130 Å x 2770 Å x 3680 Å), this domain was discretized using the dual contouring methods described previously (Song et al., 2003
; Zhang et al., 2003
) into an initial mesh containing 656,823 tetrahedral elements. Fig. 2 shows a cross section of the mesh between the mAChE surface and the outer sphere generated from the LBIE-Mesh software (http://www.ices.utexas.edu/CCV/software). Note that the mesh near the biomolecule is extremely fine and captures the details of the biomolecular surface; mesh elements increase in size with increasing distance from the biomolecule.
|
|
The electrostatic potential used for our PMF was obtained from the Poisson-Boltzmann equation using the APBS software (http://agave.wustl.edu/apbs/) (Baker et al., 2001b
). The CHARMM22 force field was used to assign the partial changes and radii of the atoms for mAChE, the dielectric values of 4 and 78 were assigned for the protein and solvent, a solvent probe radius is 1.4 Å, and an ion exclusion layer is 2.0 Å. Various ionic strengths (between 0 and 0.670 M) were used in the PMF calculations.
Additionally, the adaptive finite element methods described above were used to calculate the reaction rates for both the molecular-surface-based (Fig. 3 a) and spherical (Fig. 3 b) reactive boundary No. 1 as a function of ionic strength. Iterative error-based refinement of the initial 656,823-simplex mesh was performed until the global error was <107, a problem-value chosen to provide reaction rates which did not change appreciably upon further refinement (see Fig. 4). The reaction rate results from these calculations are shown in Fig. 5. As this figure illustrates, the spherical and molecular-surface-based reaction criteria both give results that are in good overall agreement with each other and experiment. The comparison (at 150 mM ionic strength) between the molecular and spherical boundary definitions for the 6 reactive surfaces is shown in Fig. 6. Again, the two methods are in good overall agreement but do show some differences at surfaces No. 1 and No. 2 where the differences between the reactive boundariesin particular, their surface areasare most extreme (see Fig. 3).
|
|
|
10100% from the nonadapted calculation for the molecular reactive surface and by
110% for the spherical reactive surface. Therefore, although the mesh used in the previous work (Song et al., 2003
|
Reaction rate calculation for mutant mAChE with a molecular-surface reactive boundary
One of the distinct advantages of the biomolecular reactive surface definition is the ability to directly map molecular information about the protein onto the reaction criteria. In this work, we demonstrate the ability of the biomolecular surface reactive boundary definition to correctly describe the kinetics of mAChE active site and gorge mutants. Specifically, we examined the effects of E202Q, D74N, and D74N/E202Q mutations on the mAChE ligand binding rate. The location of D74 and E202 are indicated in Fig. 3 a. D74 is located between reactive surface 3 and 4, whereas E202 is located within reactive surface 6, close to the active site. The chosen mutations constitute changes adjacent to the catalytic triad (E202Q) and in the active site gorge (D74N). These mutants have been characterized both experimentally (Radic et al., 1997
) and computationally (Tara et al., 1999
), using BD methods. The purpose of this work is to determine the ability of the molecular reactive surface boundary to quantitatively capture the effects of these mutations on the reaction rate.
The simulation protocol for these mutations is the same as for the ionic strength calculations described above, with one important difference. New electrostatic PMFs are recalculated for each of the three mutants at 150 mM ionic strength. Due to the isosteric nature of the mutations, all other parameters (particularly the mesh) remain unchanged.
To compare the continuum diffusion rate constants with BD results, simulations similar to those of Tara et al. (1998)
were repeated using the UHBD software (Madura et al., 1995
) to calculate the 150 mM ionic strength reaction rates at each of the 6 reactive surfaces for D74N, E202Q, and D74N/E202Q mutations. In particular, for each set of conditions, 5 BD runs of 200 trajectories each were simulated. The ligand was modeled as a sphere with +1 charge and a 2.0 Å radius. Each trajectory was started at a random location on a spherical surface of 55 Å (centered on the protein) and terminated when the ligand either passed the reactive boundary (see above) or a second spherical surface of radius 300 Å. The BD equation of motion was integrated using the standard Ermak-McCammon algorithm (Ermak and McCammon, 1978
) and variable time steps: 5 fs at 100 Å from the mAChE center, 1 ps at 100175 Å from the mAChE center, and 5 ps at 175300 Å from the mAChE center.
Fig. 7 presents the log-ratios of reaction rates for the wild-type, E202Q, D74N, and D74N/E202Q mAChE calculated from the adaptive SE calculations and BD simulations at 150 mM ionic strength. This figure also includes experimental values (Radic et al., 1997
) for reference. These results indicate that the continuum SE calculations generate rates with similar trends as BD; however, the actual values can differ by as much as two orders of magnitude between the two methods. Such discrepancies in predicted rates are not too surprising and simply reflect differences in methods, particularly their definitions of reactions and reactive boundaries.
|
| CONCLUSIONS AND DISCUSSIONS |
|---|
|
|
|---|
The timing information provided for the adaptive method illustrates that it is definitely more computationally-demanding than the nonadaptive calculations and, in some cases, even more expensive than the traditional BD simulations. However, judicious use of these expensive calculations can save substantial time and still provide an overall gain over BD simulations. As mentioned earlier, the initial mesh is generated based on the biomolecular geometry and does not necessarily provide the best possible basis set for solution of the SE. Therefore, adaptive refinement should be a standard part of finite element SE calculations to ensure the most accurate results. In particular, a "benchmark" adaptive refinement calculation is needed for each system studied to determine the error in rates calculated on the initial mesh and to obtain the global error tolerance at which the calculated rates are converged (cf. Fig. 4). However, much of the refinement performed in the each of the various adaptive calculations described above is redundant; i.e., the same regions of the initial mesh are refined for each system. Therefore, substantial time could be saved by reusing the refined mesh from the benchmark calculation as a starting point for subsequent simulations and thereby avoiding some of the expensive refinement steps.
Furthermore, as mentioned in the previous work, these continuum diffusion methods are designed to address a different scale of simulation from traditional BD methods; specifically, laying the groundwork for integration of molecular-scale information into cellular-scale systems (Smart and McCammon, 1998
; Tai et al., 2003
). In particular, this ability of FE-based continuum methods to integrate scales is demonstrated by performing the diffusion calculations on a finite element mesh of 0.31 µm x 0.28 µm x 0.37 µm, or
5 times the length of the BD domain (600 Å long). The adaptive methods developed in this work further facilitate our ultimate goal of multiscale modeling by enabling the efficient solution of the SE through adaptively allocating the unknowns based on error estimates.
Finally, one particularly useful aspect of the molecular surface boundary definition introduced in this work is the ability to directly connect the reactive surfaces of the simulation with the underlying biomolecule. For example, surface 6 maps directly onto the catalytic triad and are the most intuitive reactive boundary definitions available for this system. Additionally, these surfaces performed well in quantifying the impact of mutation on binding-rate constants. Such molecular-based reactive definitions suggest future possibilities of connecting coarse-grained simulations of diffusion with more detailed descriptions of enzyme function; in particular molecular dynamics and quantum mechanics worked describing the details of biomolecular binding and catalysis (Luty et al., 1993
; Zhang et al., 2002
).
| ACKNOWLEDGEMENTS |
|---|
|
|
|---|
The work at Washington University was supported by the National Science Foundation-National Partnership for Advanced Computational Infrastructure, National Biomedical Computation Resource (National Institutes of Health P41 RR08605), and an Alfred P. Sloan Research Fellowship to N.A.B. The work at the University of Texas was supported in part by National Science Foundation grants ACI-0220037, CCR-9988357, and EIA-0325550; a University of Texas-M. D. Anderson Cancer Center Whitaker grant; and a subcontract from the University of California, San Diego 1018140 as part of the National Science Foundation-National Partnership for Advanced Computational Infrastructure project.
Submitted on February 11, 2004; accepted for publication June 1, 2004.
| REFERENCES |
|---|
|
|
|---|
Allison, S. A., and J. A. McCammon. 1985. Dynamics of substrate binding to copper zinc superoxide dismutase. J. Phys. Chem. 89:10721074.[CrossRef]
Anglister, L., J. R. Stiles, B. Haesaert, J. Eichler, and M. M. Salpeter. 1995. Acetylcholinesterase at neuromuscular junctions. In Enzymes of the Cholinesterase Family. D. M. Quinn, A. S. Balasubramanian, and B. P. Doctor, editors. Plenum Press, New York.
Antosiewicz, J., and J. A. McCammon. 1995. Electrostatic and hydrodynamic orientational steering effects in enzyme-substrate association. Biophys. J. 69:5765.
Antosiewicz, J., J. A. McCammon, S. T. Wlodek, and M. K. Gilson. 1995. Simulation of charge-mutant acetylcholinesterases. Biochemistry. 34:42114219.[CrossRef][Medline]
Antosiewicz, J., S. T. Wlodek, and J. A. McCammon. 1996. Acetylcholinesteraserole of the enzyme's charge distribution in steering charged ligands toward the active site. Biopolymers. 39:8594.[CrossRef][Medline]
Axelsson, O., and V. A. Barker. 1984. Finite Element Solution of Boundary Value Problems. Theory and Computation. Academic Press, San Diego.
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:13431352.[CrossRef]
Baker, N. A., D. Sept, M. J. Holst, and J. A. McCammon. 2001a. The adaptive multilevel finite element solution of the Poisson-Boltzmann equation on massively parallel computers. IBM J. Res. Dev. 45:427438.
Baker, N. A., D. Sept, S. Joseph, M. J. Holst, and J. A. McCammon. 2001b. Electrostatics of nanosystems: Application to microtubules and the ribosome. Proc. Natl. Acad. Sci. USA. 98:1003710041.
Braess, D. 1997. Finite Elements. Theory, Fast Solvers, and Applications in Solid Mechanics. Cambridge University Press, Cambridge.
Elcock, A. H., M. J. Potter, D. A. Matthews, D. R. Knighton, and J. A. McCammon. 1996. Electrostatic channeling in the bifunctional enzyme dihydrofolate reductase-thymidylate synthase. J. Mol. Biol. 262:370374.[CrossRef][Medline]
Elcock, A. H., D. Sept, and J. A. McCammon. 2001. Computer simulation of protein-protein interactions. J. Phys. Chem. B. 105:15041518.
Ermak, D. L., and J. A. McCammon. 1978. Brownian dynamics with hydrodynamic interactions. J. Chem. Phys. 69:13521360.[CrossRef]
Gabdoulline, R. R., and R. C. Wade. 1998. Brownian dynamics simulation of protein-protein diffusional encounter. Methods 14:329341.[CrossRef][Medline]
Holst, M. 2001. Adaptive numerical treatment of elliptic systems on manifolds. Adv. Comp. Math. 15:139191.[CrossRef]
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:13191342.[CrossRef]
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:851869.[CrossRef][Medline]
Krissinel, E. B., and N. Agmon. 1996. Spherical symmetric diffusion problem. J. Comput. Chem. 17:10851098.[CrossRef]
Kurnikova, M. G., R. D. 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:642656.
Luty, B. A., S. Elamrani, and J. A. McCammon. 1993. Simulation of the bimolecular reaction between superoxide and superoxide dismutase: synthesis of the encounter and reaction steps. J. Am. Chem. Soc. 115:1187411877.[CrossRef]
Madura, J. D., J. M. Briggs, R. C. Wade, M. E. Davis, B. A. Luty, A. Ilin, J. Antosiewicz, M. K. Gilson, B. Bagheri, L. R. Scott, et al. 1995. Electrostatics and diffusion of molecules in solutionsimulations with the University of Houston Brownian Dynamics Program. Comput. Phys. Commun. 91:5795.[CrossRef]
Northrup, S. H., S. A. Allison, and J. A. McCammon. 1984. Brownian dynamics simulation of diffusion-influenced biomolecular reactions. J. Chem. Phys. 80:15171524.[CrossRef]
Radic, Z., P. D. Kirchhoff, D. M. Quinn, J. A. McCammon, and P. Taylor. 1997. Electrostatic influence on the kinetics of ligand binding to acetylcholinesterasedistinctions between active center ligands and fasciculin. J. Biol. Chem. 272:2326523277.
Roux, B., and T. Simonson. 1999. Implicit solvent models. Biophys. Chem. 78:120.[CrossRef][Medline]
Schuss, Z., B. Nadler, and R. S. Eisenberg. 2001. Derivation of Poisson and Nernst-Planck equations in a bath and channel from a molecular model. Phys. Rev. E. 64:036116.[CrossRef]
Smart, J. L., and J. A. McCammon. 1998. Analysis of synaptic transmission in the neuromuscular junction using a continuum finite element model. Biophys. J. 75:16791688.
Song, Y., Y. Zhang, T. Shen, C. L. Bajaj, J. A. McCammon, and N. A. Baker. 2003. Finite element solution of the steady-state Smoluchowksi equation for rate constant calculations. Biophys. J. 86:20172029.
Stiles, J. R., and T. M. Bartol. 2000. Monte Carlo methods for simulating realistic synaptic micro-physiology using MCell. In Computational Neuroscience: Realistic Modeling for Experimentalists. E. De Schutter, editor. CRC Press, New York.
Tai, K., S. D. Bond, H. R. MacMillan, N. A. Baker, M. J. Holst, and J. A. McCammon. 2003. Finite element simulations of acetylcholine diffusion in neuromuscular junctions. Biophys. J. 84:22342241.
Tara, S., A. H. Elcock, P. D. Kirchhoff, J. M. Briggs, Z. Radic, P. Taylor, and J. A. McCammon. 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:465474.[CrossRef][Medline]
Tara, S., T. P. Straatsma, and J. A. McCammon. 1999. Mouse acetylcholinesterase unliganded and in complex with huperzine A: a comparison of molecular dynamics simulations. Biopolymers. 50:3543.[CrossRef][Medline]
Zhang, Y., C. Bajaj, and B. S. Sohn. 2003. Adaptive and Quality 3-D Meshing from Imaging Data. Proc. ACM SM 2003. ACM Press, Seattle. 286291. (Full version to appear in CMAME).
Zhang, Y., J. Kua, and J. A. McCammon. 2002. Role of the catalytic triad and oxyanion hole in acetylcholinesterase catalysis: An ab initio QM/MM study. J. Am. Chem. Soc. 124:1057210577.[CrossRef][Medline]
This article has been cited by other articles:
![]() |
Y. Cheng, J. K. Suen, D. Zhang, S. D. Bond, Y. Zhang, Y. Song, N. A. Baker, C. L. Bajaj, M. J. Holst, and J. A. McCammon Finite Element Analysis of the Time-Dependent Smoluchowski Equation for Acetylcholinesterase Reaction Rate Calculations Biophys. J., May 15, 2007; 92(10): 3397 - 3406. [Abstract] [Full Text] [PDF] |
||||
![]() |
D. Zhang, J. Suen, Y. Zhang, Y. Song, Z. Radic, P. Taylor, M. J. Holst, C. Bajaj, N. A. Baker, and J. A. McCammon Tetrameric Mouse Acetylcholinesterase: Continuum Diffusion Rate Calculations by Solving the Steady-State Smoluchowski Equation Using Finite Element Methods Biophys. J., March 1, 2005; 88(3): 1659 - 1665. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |