^{1}, Michael D. Altman

^{2}, David J. Willis

^{3}, Shaun M. Lippow

^{4}, Bruce Tidor

^{5,a),b)}and Jacob K. White

^{6,a),c)}

### Abstract

Surface formulations of biophysical modeling problems offer attractive theoretical and computational properties. Numerical simulations based on these formulations usually begin with discretization of the surface under consideration; often, the surface is curved, possessing complicated structure and possibly singularities. Numerical simulations commonly are based on approximate, rather than exact, discretizations of these surfaces. To assess the strength of the dependence of simulation accuracy on the fidelity of surface representation, here methods were developed to model several important surface formulations using exact surface discretizations. Following and refining Zauhar’s work [J. Comput.-Aided Mol. Des.9, 149 (1995)], two classes of curved elements were defined that can exactly discretize the van der Waals, solvent-accessible, and solvent-excluded (molecular) surfaces. Numerical integration techniques are presented that can accurately evaluate nonsingular and singular integrals over these curved surfaces. After validating the exactness of the surface discretizations and demonstrating the correctness of the presented integration methods, a set of calculations are presented that compare the accuracy of approximate, planar-triangle-based discretizations and exact, curved-element-based simulations of surface-generalized-Born (sGB), surface-continuum van der Waals (scvdW), and boundary-element method(BEM)electrostatics problems. Results demonstrate that continuum electrostatic calculations with BEM using curved elements, piecewise-constant basis functions, and centroid collocation are nearly ten times more accurate than planar-triangle BEM for basis sets of comparable size. The sGB and scvdW calculations give exceptional accuracy even for the coarsest obtainable discretized surfaces. The extra accuracy is attributed to the exact representation of the solute-solvent interface; in contrast, commonly used planar-triangle discretizations can only offer improved approximations with increasing discretization and associated increases in computational resources. The results clearly demonstrate that the methods for approximate integration on an exact geometry are far more accurate than exact integration on an approximate geometry. A MATLAB implementation of the presented integration methods and sample data files containing curved-element discretizations of several small molecules are available online as supplemental material.

The authors are indebted to J.-H. Lee for useful discussions and to X. Wang for sharing his monomial integration software. This work was supported by the National Institutes of Health (GM065418 and CA096504), the Singapore-MIT Alliance, and the National Science Foundation. One of the authors (D.J.W.) gratefully acknowledges support from the Natural Sciences and Engineering Research Council of Canada. Two of the authors (J.P.B. and M.D.A.) contributed equally to this work.

I. INTRODUCTION

II. BACKGROUND

A. Surface formulations of biophysical problems

1. Molecular electrostatics

2. Surface-generalized Born

3. Continuum van der Waals

B. Defining molecule-solvent interfaces

III. SURFACE DISCRETIZATION

A. Toroidal element definition

B. Spherical element definition

IV. CURVED-ELEMENT INTEGRATION METHODS

A. Far-field quadrature

1. Generalized spherical triangle coordinate transformation

2. Toroidal element coordinate transformation

B. Near-field integration techniques

1. Single-layer potential

2. Double-layer potential

V. RESULTS

A. Problem geometries

1. Alanine tripeptide

2. Alanine dipeptide

3. Barnase-barstar complex

B. Validating the surface discretization

C. Validating curved boundary-element integration

D. Surface-generalized-Born calculations

E. Continuum van der Waals calculations

F. Poisson-Boltzmann electrostatics problems

1. Spherical geometry

2. Alanine dipeptide

G. Performance

VI. DISCUSSION

### Key Topics

- Boundary element methods
- 12.0
- Electrostatics
- 12.0
- Solvents
- 12.0
- Integral equations
- 9.0
- Integral transforms
- 9.0

## Figures

A mixed discrete-continuum model for biomolecule electrostatics. The surface represents the dielectric boundary between regions with dielectric constants and . Partial atomic charges are located in region I, with illustrative charges at and at . The Debye screening parameter is zero within region I and may be nonzero in region II. In work not described here, an ion-exclusion layer may also be treated (Refs. 29 and 30).

A mixed discrete-continuum model for biomolecule electrostatics. The surface represents the dielectric boundary between regions with dielectric constants and . Partial atomic charges are located in region I, with illustrative charges at and at . The Debye screening parameter is zero within region I and may be nonzero in region II. In work not described here, an ion-exclusion layer may also be treated (Refs. 29 and 30).

Three definitions of solute-solvent boundaries, shown for a two-atom case: (a) van der Waals surface, (b) the Lee and Richards solvent-accessible surface, and (c) the Richards solvent-excluded (molecular) surface. The dotted lines in (b) and (c) denote the van der Waals surface.

Three definitions of solute-solvent boundaries, shown for a two-atom case: (a) van der Waals surface, (b) the Lee and Richards solvent-accessible surface, and (c) the Richards solvent-excluded (molecular) surface. The dotted lines in (b) and (c) denote the van der Waals surface.

Specification of a torus and a torus element with and .

Specification of a torus and a torus element with and .

A generalized spherical triangle (GST) with one bounding edge belonging to the circle centered at the blue dot. The remaining edges belong to great circles on the sphere.

A generalized spherical triangle (GST) with one bounding edge belonging to the circle centered at the blue dot. The remaining edges belong to great circles on the sphere.

(a) The standard unit triangle in parametric coordinate space. (b) A GST viewed from the negative axis. The angle is measured relative to the positive axis. Each is mapped to one plane with normal along the axis; the plane intersects the sphere and defines a circle. (c) A GST viewed from the positive axis. The dashed lines indicate the circle of intersection between the sphere surface and the plane specified by . The angle specifies the rotation about the axis. The image of the standard-triangle vertices under the coordinate transformation are labeled.

(a) The standard unit triangle in parametric coordinate space. (b) A GST viewed from the negative axis. The angle is measured relative to the positive axis. Each is mapped to one plane with normal along the axis; the plane intersects the sphere and defines a circle. (c) A GST viewed from the positive axis. The dashed lines indicate the circle of intersection between the sphere surface and the plane specified by . The angle specifies the rotation about the axis. The image of the standard-triangle vertices under the coordinate transformation are labeled.

Schematic of the approach for evaluating the potential induced by a distribution of monopole charge on a generalized spherical triangle. A planar reference element viewed edge on (blue) is defined to be tangent to the original GST (black solid arc) at the GST centroid. The reference element vertices are defined to be the projection of the GST vertices to the plane tangent to the GST centroid.

Schematic of the approach for evaluating the potential induced by a distribution of monopole charge on a generalized spherical triangle. A planar reference element viewed edge on (blue) is defined to be tangent to the original GST (black solid arc) at the GST centroid. The reference element vertices are defined to be the projection of the GST vertices to the plane tangent to the GST centroid.

The Newman approach to calculating the potential induced by a uniform distribution of a normally oriented dipole charge layer (Ref. 13). The cross at the center of the sphere denotes the point at which the potential is to be determined; the thin arcs form the edges of a GST; the thick lines represent the projection of the GST bounding arcs to the sphere. The double-layer potential, which equals the solid angle bounded by the thick lines, is directly proportional to the bounded area.

The Newman approach to calculating the potential induced by a uniform distribution of a normally oriented dipole charge layer (Ref. 13). The cross at the center of the sphere denotes the point at which the potential is to be determined; the thin arcs form the edges of a GST; the thick lines represent the projection of the GST bounding arcs to the sphere. The double-layer potential, which equals the solid angle bounded by the thick lines, is directly proportional to the bounded area.

Edges may be projected to the unit sphere regardless of their position relative to it. (left) The edges are outside the unit sphere. (right) The edges are interior to the unit sphere.

Edges may be projected to the unit sphere regardless of their position relative to it. (left) The edges are outside the unit sphere. (right) The edges are interior to the unit sphere.

Generalized-Born radii calculated by volume integration and by evaluating surface integrals based on the GB model proposed by Grycuk (Ref. 39). The volume radii are plotted as large squares and the surface GB radii are plotted with circles, triangles, crosses, and diamonds. (a) Alpha-helix blocked alanine tripeptide. (b) Beta-sheet blocked alanine tripeptide.

Generalized-Born radii calculated by volume integration and by evaluating surface integrals based on the GB model proposed by Grycuk (Ref. 39). The volume radii are plotted as large squares and the surface GB radii are plotted with circles, triangles, crosses, and diamonds. (a) Alpha-helix blocked alanine tripeptide. (b) Beta-sheet blocked alanine tripeptide.

Convergence of solvation free energies for a centrally located charge in a sphere, calculated by BEM numerical solution of the Yoon and Lenhoff integral equations. For both cases and . (a) . (b) .

Convergence of solvation free energies for a centrally located charge in a sphere, calculated by BEM numerical solution of the Yoon and Lenhoff integral equations. For both cases and . (a) . (b) .

Solvation free energies for four conformers of the alanine dipeptide; atom centers are those presented in Ref. 61 and PARSE atomic radii and partial charges have been used (Ref. 62). (a) c5 geometry. (b) geometry. (c) c7ax geometry. (d) c7eq geometry.

Solvation free energies for four conformers of the alanine dipeptide; atom centers are those presented in Ref. 61 and PARSE atomic radii and partial charges have been used (Ref. 62). (a) c5 geometry. (b) geometry. (c) c7ax geometry. (d) c7eq geometry.

## Tables

Comparison of discretized surface areas with analytical molecular (solvent-excluded) surface area. Probe radius is taken to be . All area quantities are in and have been rounded to the nearest .

Comparison of discretized surface areas with analytical molecular (solvent-excluded) surface area. Probe radius is taken to be . All area quantities are in and have been rounded to the nearest .

Comparison of discretized surface areas with analytical solvent-accessible surface area. Probe radius is taken to be . All area quantities are in and have been rounded to the nearest .

Comparison of discretized surface areas with analytical solvent-accessible surface area. Probe radius is taken to be . All area quantities are in and have been rounded to the nearest .

Comparison of pit, belt, and cap areas computed by analytical, direct quadrature, and polynomial-fitting methods, using the molecular surface discretizations of Sec. V B. All area quantities are in and have been rounded to the nearest .

Comparison of pit, belt, and cap areas computed by analytical, direct quadrature, and polynomial-fitting methods, using the molecular surface discretizations of Sec. V B. All area quantities are in and have been rounded to the nearest .

Solute-solvent van der Waals interaction energies estimated using a volume integration scheme and using a surface formulation of the continuum van der Waals model of Levy *et al.* (Ref. 8) and curved surface elements. All energies are in kcal/mol and have been rounded to the nearest .

Solute-solvent van der Waals interaction energies estimated using a volume integration scheme and using a surface formulation of the continuum van der Waals model of Levy *et al.* (Ref. 8) and curved surface elements. All energies are in kcal/mol and have been rounded to the nearest .

Approximate number of panel integrals per second computable using a C implementation of the presented techniques. All planar-triangle integrals have been computed by the methods of Hess and Smith (Ref. 12) or Newman (Ref. 13). All far-field curved-element integrals have been computed using a 16-point quadrature rule. Self-term and near-field GST integrals have been computed using the polynomial-fitting scheme. Self-term and near-field toroidal element integrals have been computed by recursive subdivision. A single-layer integral has been defined to be in the far field if the separation between the evaluation (field) point and the curved element exceeds four times the length of the element’s longest edge; a double-layer integral is defined to be in the far field if the separation exceeds twice the length of the longest edge.

Approximate number of panel integrals per second computable using a C implementation of the presented techniques. All planar-triangle integrals have been computed by the methods of Hess and Smith (Ref. 12) or Newman (Ref. 13). All far-field curved-element integrals have been computed using a 16-point quadrature rule. Self-term and near-field GST integrals have been computed using the polynomial-fitting scheme. Self-term and near-field toroidal element integrals have been computed by recursive subdivision. A single-layer integral has been defined to be in the far field if the separation between the evaluation (field) point and the curved element exceeds four times the length of the element’s longest edge; a double-layer integral is defined to be in the far field if the separation exceeds twice the length of the longest edge.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content