Computer Aided Drug Design
Electrostatic potentials of proteins in water: a structured continuum approach
1 Center for Bioinformatics, Saarland University PO 15 11 50, 66041 Saarbrücken, Germany
2 Interdisciplinary Research Institute c/o IEMN, Cité Scientifique BP 69 F-59652 Villeneuve d'Ascq, France
3 Department of Mathematics, Saarland University PO 15 11 50, 66041 Saarbrücken, Germany
4 Department for Simulation of Biological Systems, WSI/ZBIT, University of Tübingen, Sand 14 72070 Tübingen, Germany
*To whom correspondence should be addressed.
 |
ABSTRACT |
|---|
 TOP ABSTRACT 1 INTRODUCTION2 INTEGRATING WATER STRUCTURE...3 THE BOUNDARY ELEMENT...4 RESULTS5 DISCUSSIONREFERENCES |
|---|
Electrostatic interactions play a crucial role in many biomolecular processes, including molecular recognition and binding. Biomolecular electrostatics is modulated to a large extent by the water surrounding the molecules. Here, we present a novel approach to the computation of electrostatic potentials which allows the inclusion of water structure into the classical theory of continuum electrostatics. Based on our recent purely differential formulation of nonlocal electrostatics [Hildebrandt, et al. (2004) Phys. Rev. Lett., 93, 108104] we have developed a new algorithm for its efficient numerical solution. The key component of this algorithm is a boundary element solver, having the same computational complexity as established boundary element methods for local continuum electrostatics. This allows, for the first time, the computation of electrostatic potentials and interactions of large biomolecular systems immersed in water including effects of the solvent's structure in a continuum description. We illustrate the applicability of our approach with two examples, the enzymes trypsin and acetylcholinesterase. The approach is applicable to all problems requiring precise prediction of electrostatic interactions in water, such as protein–ligand and protein–protein docking, folding and chromatin regulation. Initial results indicate that this approach may shed new light on biomolecular electrostatics and on aspects of molecular recognition that classical local electrostatics cannot reveal.
Contact: anhi{at}bioinf.uni-sb.de
|
1 INTRODUCTION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 INTEGRATING WATER STRUCTURE...3 THE BOUNDARY ELEMENT...4 RESULTS5 DISCUSSIONREFERENCES |
|---|
Electrostatic interactions constitute a fundamental driving force underlying many biological processes on molecular scales and play an important role in molecular recognition, e.g. in enzymatic processes, or more generally in protein–ligand or protein–protein interactions. A physically precise and, at the same time, computationally efficient method is therefore highly desirable, also for applications in structure-based drug design.
Despite many technical and computational advances in recent years, the precise determination of electrostatic fields of proteins or other biological macromolecules is still a challenge, both from a theoretical and a computational point of view. The main problem lies in the fact that biomolecules reside in a watery environment, the dielectric properties on a molecular scale of which are strongly determined by the dipolar character of the water molecules and the complex spatial network they form.
Due to their high computational demand, approaches to include water structure around proteins based on Molecular Dynamic simulations of proteins with explicit water typically suffer from long run times. On the other hand, the application of macroscopic continuum electrostatics (Nicholls and Honig, 1990; Baker et al., 2001; Chen and Chipman, 2003; Totrov and Abagyan, 2001) completely neglects the effects of the solvent's structure, which is known to be of fundamental importance.
We have developed a novel approach to protein electrostatics based on an extension of classical continuum electrostatics accounting for the structure of the water surrounding a biomolecule (Hildebrandt et al., 2004). The main feature of this approach is an extension of the concept of the dielectric constant which can be consistently generalized to a dielectric function of two spatial arguments, capturing the dielectric effect of the water–water correlation. The resulting continuum theory of electrostatics, known in the physics literature as ‘nonlocal electrostatics’, can be formulated as a system of partial differential equations which generalizes the classical local theory in a transparent way. Furthermore, this formulation readily allows for a discretization in terms of a boundary element method (BEM), in contrast to the original nonlocal formulation which is based on difficult-to-handle integro-differential Equations. Our current implementation has the same computational complexity as previously published BEM approaches to local protein electrostatics.
To investigate the properties of our novel nonlocal electrostatic model, we have derived analytical solutions for geometrically feasible situations, i.e. for the electrostatic potential and the free energy of solvation for monoatomic ions. Comparison of the results to experimental data shows excellent agreement (Hildebrandt et al., 2004). In (Hildebrandt, 2005) we compared the analytical solutions for this model with the numerical solutions computed with our BEM solver and obtained an almost perfect correlation. Application to the free energy of solvation for small molecules (neutral amino acid side chain analogs) also resulted in highly accurate prediction of the experimental values (Hildebrandt, 2005). In the present work, we demonstrate that the efficient algorithm presented here allows the application of this method to large macromolecular systems (e.g. whole proteins or protein complexes) as well. We have studied the nonlocal electrostatic potentials of trypsin and mouse acetylcholinesterase (AChE) and compared them with their local counterparts, yielding interesting insights into the nature of the nonlocal effect. Our results indicate that the water structure has two main effects on the electrostatic potential. First, as expected, the nonlocal electrostatic potentials show a greatly extended visibility—the potentials extend much further into space than in the local case. Second, even for larger distances from the molecular surface, the electrostatic potential features a much more pronounced structure in the nonlocal case than in the local case, due to water–water correlation. Both effects have significant implications for molecular recognition.
The structure of this paper is as follows. In Section 2, we briefly introduce the theoretical background of our structured continuum approach. The BEM method developed for the new theory is described in Section 3. In Section 4, we present results on two test systems. A discussion in Section 5 concludes the paper.
|
2 INTEGRATING WATER STRUCTURE INTO CLASSICAL ELECTROSTATICS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 INTEGRATING WATER STRUCTURE...3 THE BOUNDARY ELEMENT...4 RESULTS5 DISCUSSIONREFERENCES |
|---|
Classical electrostatics is based on the Poisson equation
 | (1) |
is the Laplacian, 
is the electrostatic potential, 
the density of charges, and the coefficients 
0 and are the dielectric constants of the vacuum and of the water surrounding the protein, respectively. In this form the equation is justifiable only for macroscopic bodies in which a spatial average over the exterior volume is assumed. This is clearly not the case when we want to use Equation (1) to determine electrostatic potentials on smaller, let alone atomic, scales.
One way to attack this problem is to include water structure in the dielectric properties by introducing an effective dielectric function (r, r') which in general depends on two spatial arguments. It describes the spatial dispersion of dielectric properties on length scales on which the hydrogen-bonded network of the water molecules displays significant correlations in its orientational polarization. As a result, the function (r, r') will depend on the length scale 
of these correlations, which is typically in the range of 10–20 Å. While the theory does not allow to distinguish the contributions of individual molecules, but rather describes a mean-field effect over a range set by the parameter , it allows to obtain a more realistic insight into the overall effect of the surrounding water on the solvated protein.
The introduction of (r, r') into Equation (1) comes at a price, since the equation now reads as
 | (2) |
The Poisson equation has turned into an integro-differential equation which is technically very demanding. It can be tackled by analytical means only for cases of high symmetry (sphere, cylinder). For more general geometries, the complexity of the above equation renders even numerical methods infeasible.
We have recently developed a novel approach to Equation (2) which makes it tractable to numerical methods, as shown below. The basic setup is illustrated in Figure 1. By treating the protein in water as a body with a dielectric constant 
2–5 embedded in a nonlocal medium obeying Equation (2) with 
0 and taking into account the correct interface conditions—which also become nonlocal due to the presence of (r, r')—we have shown that one can recast the standard formulation of nonlocal electrostatics into a system of partial differential equations, if
- the dielectric function (r, r') can be related to the Green's function of a known differential operator L (this condition is fulfilled under very general conditions);
- the dielectric displacement field D can be decomposed into an irrotational and a solenoidal part, i.e. D = –
+ 
;
- the electrostatic potential is accompanied by the electrostatic potential which fulfills the Poisson equation inside the protein and is coupled to outside the protein by a third independent partial differential equation involving the explicit form of L.
|
The full form of the novel formulation of nonlocal electrostatics for the special case of
was derived in (Hildebrandt et al., 2004) and was illustrated there for the test case of singly and doubly charged ions. Simplifying one of the two interface conditions, the resulting system of differential equations looks as follows:
|
| (3) |
 |
. |
3 THE BOUNDARY ELEMENT SOLVER |
|---|
|
TOPABSTRACT1 INTRODUCTION2 INTEGRATING WATER STRUCTURE...3 THE BOUNDARY ELEMENT...4 RESULTS5 DISCUSSIONREFERENCES |
|---|
While the purely differential nature of our novel formulation of nonlocal electrostatics greatly simplifies analytical computations in the nonlocal framework, it is much more important from a practical point of view that it allows, for the first time, the development of efficient and highly accurate numerical solvers. It is this property which makes the application of the nonlocal theory to large and structurally complex systems possible, arising in real-world applications in computational biology—something that seemed completely infeasible in the classical formulation.
In principle, several numerical techniques can be applied to solve the above system of differential equations, such as finite difference methods, finite element methods and BEMs. For the conventional local electrostatics of biomolecular systems, it has been shown that boundary element schemes can provide very accurate numerical solvers that can be made very efficient (Chen and Chipman, 2003; Totrov and Abagyan, 2001). Hence, we decided to develop such a scheme for the novel formulation of nonlocal electrostatics.
Boundary element methods replace differential equations in the volume by integral equations on the boundary surface of the system, i.e. in our case on the surface 
separating the biomolecule from the solvent. But this reduction in the dimensionality of the problem comes at a cost: the method is not generally applicable. The differential equation of interest must possess a fundamental solution which must be known analytically—fortunately, this holds in our case. Furthermore, the explicit formulation requires lengthy and complicated analytical preparations. The differential equation has to be multiplied by its fundamental solution and the differential operator has to be shifted from the unknown function we seek to the known fundamental solution using analogues of Green's theorem. Due to space constraints, the complete derivation of these results cannot be presented here. The interested reader is referred to (Hildebrandt, 2005).
The boundary integral operators appearing in the resulting equations are subsequently discretized by approximating the potentials of interest as low-rank polynomials over the triangles of the surface discretization. This step reduces the involved systems of integral equations to simple linear systems of equations which can be solved efficiently using standard techniques from computational linear algebra. The resulting system has the following form:
|
| (4) |
, 
, 
, 
, 
and 
are matrices of size n x n, where n denotes the number of triangles of the surface discretization and u, q, w are vectors of size n. The solution of the resulting linear system then directly yields the nonlocal electrostatic potential on the boundary, i.e. on the molecular surface, from which all relevant potentials and fields can be computed easily and efficiently everywhere in space, using the so-called representation formulae.
Comparing system (4) with the well-known local case,
 | (5) |
We have implemented a solver for system (4) in C++ that uses the ATLAS library (Whaley et al., 2001) for efficient solution of the linear system of equations.
|
4 RESULTS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 INTEGRATING WATER STRUCTURE...3 THE BOUNDARY ELEMENT...4 RESULTS5 DISCUSSIONREFERENCES |
|---|
We validated the implementation described above by comparing numerical results for spherical ions with the analytical solution taken from (Hildebrandt et al., 2004). Numerical and analytical results were in perfect agreement with a correlation coefficient of r
0.99992. We have also tested the BEM solver by predicting small molecule solvation free energies (Hildebrandt, 2005). We chose amino acid side chain analogs as a test case. Again, we obtained a very good agreement with the experimental values (data not shown).
With the current implementation, the application to systems with 
12 000–15 000 triangles is already feasible. While a larger number of triangles would be desirable (see Section 5), this is sufficient for qualitative results on small- to medium-sized proteins. In the present work, we discuss the application to the two enzymes trypsin and AChE.
During these tests, we found that generating a triangulation of the molecular surface suited to our numerical purposes is not as straightforward as it might seem. While, in principle, a variety of programs for the computation of molecular surfaces exist, most of those are optimized for visualization and not for numerical purposes. In our numerical experiments, we achieved the best results with the program SMART (Zauhar, 1995). While the triangulations generated by SMART are generally of a very high quality from a numerical point of view, the number of triangles they contain turned out to be too large for our currently implemented solver. Thus, we coarsened the resulting meshs with the help of the GTS library (http://gts.sourceforge.net/) until they contained about 12 000 triangles. The triangulations generated in this way were then checked for correctness using the program Netgen (http://www.hpfem.jku.at/netgen/).
Finally, the resulting triangulations were used as input to our numerical boundary element solver, and the nonlocal electrostatic potential was computed on a three-dimensional grid surrounding the respective proteins. Setting up and solving the boundary element systems took about 2 h on an AMD Opteron processor. Subsequent computation of the potential grids with a spatial resolution of 1 Å took about 2 h on a cluster of 36 Intel Pentium 4 machines.
4.1 Trypsin
The structure of bovine trypsin was taken from the PDB (Berman et al., 2000) [PDB-ID 2PTC
[PDB]
(Marquart et al., 1983)]. Charges and radii were taken from (Sitkoff et al., 1994), and hydrogens were added with the BALL library (Kohlbacher and Lenhof, 2000). The qualitative effects of the inclusion of water structure on the electrostatic potential of a protein can be seen in Figure 2, where the contour surfaces of the electrostatic potential in the conventional local electrostatic framework are compared with the contour surfaces obtained from the nonlocal model. The contour surfaces correspond to values of = –0.1V (red) and = 0.1V (blue). Comparing the local with the nonlocal results, drastic changes can be observed.
|
As expected, the contour surfaces of the electrostatic potential obtained from the nonlocal approach reach out much further into space than their counterparts in the local description. This is due to water–water correlations, leading to a reduced orientational polarizability of the solvent and thus to a reduction in shielding of the electrostatic potential, an effect that has often been proposed in the literature (Mehler and Solmajer, 1991). We can thus state that the first ‘nonlocal’ effect is a greatly extended electrostatic visibility for possible binding partners.
In addition to this increase in electrostatic visibility, we can also recognize an interesting electrostatic structure in the vicinity of trypsin's binding pocket: two negative ‘arms’ stick out of the ‘wall’ of positive potential, which builds a bridge-like structure in this area of the active site. This seems to indicate that a second (and in our view more profound) effect of the inclusion of water–water correlations is the appearance of more pronounced structures of the electrostatic potential.
4.2 Acetylcholinesterase
Structure, residue protonation state, atom charges and radii of murine AChE were taken from the set of test cases of APBS (Baker et al., 2001) and are based on the PDB ID 1MAH
[PDB]
(Bourne et al., 1995).
AChE catalyzes the reaction of the neurotransmitter acetylcholine to choline and acetate at an extremely high reaction rate. It has been argued that one reason for this enormous efficiency is the particularly strong electrostatic field over AChE's active site (Antosiewicz et al., 1996; Tan et al., 1993).
To demonstrate the effect of the solvent structure on the electrostatics of AChE, the local and nonlocal approximations of the potential (given in Volt) are plotted color coded on a cutting plane through the active site of the enzyme in Figure 3, where the left-hand side represents the local potential and the right-hand side the nonlocal one. The interior of AChE is colored in dark blue, and the location of its active site gorge is indicated by the red arrow.
|
To simplify an intuitive understanding of the magnitude of the effect, the contour lines shown in the plot correspond to regions of multiples of the thermal energy of a singly charged probe at room temperature (kBT/e). The number of contour lines visible in the local case in Figure 3 is lower than in the nonlocal case, as in the local case most of the contour lines (those for larger potentials) coincide with the protein's surface.
As in the case of trypsin, we see that the solvent–solvent correlations, i.e. the water structure around the protein, have a drastic influence on the range of ‘electrostatic visibility’: the contour lines extend much further into space in the nonlocal case than in the local approximation.
But even more intruiging is the electrostatic structure in the vicinity of the binding pocket: not only is the relative strength of the potential greatly increased as compared with a structureless description, but also the density of contour lines. From this result, we can directly infer that the gradient of the potential at these positions, and thus the electric field, will be much greater than in the local description, leading to much stronger forces toward the interior of the binding pocket.
|
5 DISCUSSION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 INTEGRATING WATER STRUCTURE...3 THE BOUNDARY ELEMENT...4 RESULTS5 DISCUSSIONREFERENCES |
|---|
In this work, we have presented a method for the efficient computation of electrostatic potentials in a nonlocal structured continuum approach. This approach accounts for the structure of water arising from water–water correlation and thus yields a more accurate picture of electrostatics close to biomolecular surfaces than classical (local) continuum approaches. To our knowledge, this is the first such approach and the first implementation of a boundary element solver for the nonlocal Poisson equation for arbitrary geometries. Its efficiency allows the application to large-scale systems of biological interest.
The solver yields quantitatively accurate results for small molecules and although its application to larger systems, e.g. proteins, is computationally very demanding, the results are consistent with what one would expect from theoretical considerations. The approach is thus promising and opens up new avenues toward the accurate modeling and simulation of water structure on biomolecular electrostatics. Our experiments hint at a largely increased electrostatic visibility. Furthermore, the potentials display interesting structures not seen in the results obtained from local theory. We suspect that these structures might help explain enzyme efficiency and reaction rates due to larger visibility and the large forces often exerted on docking partners.
Two clear disadvantages of our current approach are the large memory and CPU requirements. We are currently working on a more efficient implementation based on the adaptive cross-approximation technique (Benendorf and Rjasanow, 2003a,b). Furthermore, we are working on the application of our approach to protein-docking problems, in the expectation that it will lead to improved energetic filters.
|
REFERENCES |
|---|
|
TOPABSTRACT1 INTRODUCTION2 INTEGRATING WATER STRUCTURE...3 THE BOUNDARY ELEMENT...4 RESULTS5 DISCUSSIONREFERENCES |
|---|
Antosiewicz, J., et al. (1996) Acetylcholinesterase: role of the enzyme's charge distribution in steering charged ligands toward the active site. Biopolymers, 39, 85–94[CrossRef][Web of Science][Medline].
Baker, N., et al. (2001) Electrostatics of nanosystems: application to microtubules and the ribosome. Proc. Natl Acad. Sci. USA, 98, 10037–20041
Bebendorf, M. and Rjasanow, S. (2003a) Adaptive low–rank approximation of collocation matrices. Computing, 70, 1–24[CrossRef].
Bebendorf, M. and Rjasanow, S. (2003b) Numerical simulation of exhaust systems in car industry—efficient calculation of radiation heat transfer. In Jäger, W. and Krebs, H.-J. (Eds.). Mathematics Key Technology for the Future, , Berlin Springer, pp. 55–62.
Berman, H., et al. (2000) The Protein Data Bank. Nucleic Acids Res, . 28, 235–242
Bourne, Y., et al. (1995) Acetylcholinesterase inhibition by fasciculin: crystal structure of the complex. Cell, 83, 503[CrossRef][Web of Science][Medline].
Chen, F. and Chipman, D. (2003) Boundary element methods for dielectric cavity construction and integration. J. Chem. Phys, . 119, 10289–10297[CrossRef].
Hildebrandt, A., et al. (2004) Novel formulation of nonlocal electrostatics. Phys. Rev. Lett, . 93, 108104[CrossRef][Medline].
Hildebrandt, A. (2005) Biomolecules in a structured solvent—a novel formulation of nonlocal electrostatics and its numerical solution. , Germany PhD Thesis, Saarland University, Saarbrücken.
Kohlbacher, O. and Lenhof, H.-P. (2000) BALL—rapid software prototyping in computational molecular biology. Bioinformatics, 16, 815–824
Marquart, M., et al. (1983) The geometry of the reactive site and of the peptide groups in trypsin, trypsinogen and its complexes with inhibitors. Acta Crystallogr. B, 39, 480[CrossRef].
Mehler, E. and Solmajer, T. (1991) Electrostatic effects in proteins: comparison of dielectric and charge models. Protein Eng, . 8, 903–910.
Nicholls, A. and Honig, B. (1990) A rapid finite difference algorithm, utilizing successive over-relaxation to solve the Poisson-Boltzmann equation. J. Comput. Chem, . 12, 435–445[CrossRef].
Sitkoff, D., et al. (1994) Accurate calculation of hydration free energies using macroscopic solvent models. J. Phys. Chem, . 98, 1978–1988[CrossRef].
Tan, R., et al. (1993) Acetylcholinesterase: electrostatic steering increases the rate of ligand binding. Biochemistry, 32, 401–403[CrossRef][Medline].
Totrov, M. and Abagyan, R. (2001) Rapid boundary element solvation electrostatics calculations in folding simulations: successful folding of a 23–residue peptide. Biopolymers, 60, 124–133[CrossRef][Web of Science][Medline].
Whaley, R.C., et al. (2001) Automated empirical optimization of software and the ATLAS project. Parallel Comput, . 27, 3–35[CrossRef].
Zauhar, R.J. (1995) Smart: a solvent-accessible triangulated surface generator for molecular graphics and boundary element applications. J. Comput. Aided Mol. Des, . 9, 149–159[CrossRef][Web of Science][Medline].
CiteULike 
Connotea 
Del.icio.us What's this?
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||