An Algebraic Spline Model of Molecular Surfaces - ACM Digital Library
An Algebraic Spline Model of Molecular Surfaces Wenqi Zhao∗ ICES, U. Texas at Austin, Austin, TX 78712
(a)
Guoliang Xu† ICMSEC, Chinese Academy of Sciences, Beijing, 100080, China
(b)
Chandrajit Bajaj ‡ CS and ICES, U. Texas at Austin, Austin, TX 78712
(c)
(d)
Figure 1: Molecular models of a protein(1HIA). (a) The van der Waals surface model (693 atoms); (b) An initial triangulated solvent excluded surface (SES) model (27480 triangles); (c) The decimated triangulated SES model(7770 triangles); (d) Our algebraic spline molecular surface model (7770 patches) generated from (c).
Abstract
1
In this paper, we describe a new method to generate a smooth algebraic spline (AS) model approximation of the molecular surface (MS), based on an initial coarse triangulation derived from the atomic coordinate information of the biomolecule, resident in the PDB (Protein data bank). Our method first constructs a triangular prism scaffold Ps covering the PDB structure, and then generates piecewise polynomial Bernstein-Bezier (BB) spline function approximation F within Ps , which are nearly C1 everywhere. Approximation error and point sampling convergence bounds are also computed. An implicit AS model of the MS which is free of singularity, is extracted as the zero contours of F. Furthermore, we generate a polynomial parametrization of the implicit MS, which allows for an efficient point sampling on the MS, and thereby simplifies the accurate estimation of integrals needed for electrostatic solvation energy calculations.
The computation of electrostatic solvation energy (also known as polarization energy) for biomolecules plays an important role in the molecular dynamics simulation [Karplus and McCammon 2002], the analysis of stability in protein structure prediction [Srinivasan et al. 1998], and the protein-ligand binding energy calculation [Kuhn and Kollman 2000]. The explicit model of the solvent provides the most rigorous solvation energy calculation [Nina et al. 1997]. However, due to the large amount of solvent molecules, most of the computer time is spent on the trajectories of the solvent molecules, which greatly reduces the efficiency of this method [Roux and Simonson 1999]. An alternative method is to represent the solvent implicitly as a dielectric continuum [Schaefer and Karplus 1996] which can be modeled by the Poisson-Boltzmann (PB) equation [Baker et al. 2000; Madura et al. 1995]. A more efficient approximation to the PB electrostatic solvation energy is the generalized Born (GB) model [Still et al. 1990; Bashford and Case 2000; Lee et al. 2003], which approximates the electrostatic solvation energy ∆Gelec as qi q j τ Gpol = − ∑ (1) 2 i, j [r2 + R R exp(− ri2j )] 12
CR Categories: I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling—Curve, surface, solid, and object representations; J.3 [Life and Medical Sciences]: Biology and genetics
Introduction
ij
1 εp
i
j
FRi R j
1 εw ,
where τ = − ε p is the solute (low) dielectric constant, εw is the solvent (high) dielectric constant, qi is the atomic charge of atom i, ri j is the distance between atom i and j, F is an empirical factor (4 [Still et al. 1990], or 8 [Lee et al. 2003]), and Ri is
Keywords: Polynomial splines, molecular surfaces, prismatic scaffolds, Bernstein-Bezier basis, solvation, energetics, error bounds, rate of convergence ∗ Supported in part by NSF grant IIS-0325550, and NIH contract R01GM07308.Email: [email protected] † Support in part by NSFC grant 10371130, National Key Basic Research Project of China (2004CB318000), the J.T. Oden ICES visitor fellowship, and NSF grant IIS-0325550. Email:[email protected] ‡ Supported in part by NSF grants IIS-0325550, CNS0540033 and NIH contracts P20-RR020647, R01-9M074258, R01GM07308.Email:[email protected]
Copyright © 2007 by the Association for Computing Machinery, Inc. Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, to republish, to post on servers, or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from Permissions Dept, ACM Inc., fax +1 (212) 869-0481 or e-mail [email protected]. SPM 2007, Beijing, China, June 04 – 06, 2007. © 2007 ACM 978-1-59593-666-0/07/0006 $5.00
297
the effective Born radius of atom i. The effective Born radius reflects how deeply an atom is buried in the molecule: the deeper it is buried, the larger the radius is and consequently the less important to the polarization. The effective Born radii can be computed by the surface integral [Ghosh et al. 1998] R−1 i =
1 4π
Z Γ
(r − xi ) · n(r) dS |r − xi |4
However the error of the generated isosurface could be very large and result in inaccurate energy computation. Bajaj et al [Bajaj et al. 1997] presented a NURBS representation for the surface. Though it provides a parametric approximation to the SES, it still has singularity over the surface. Edelsbrunner [Edelsbrunner 1999] defines another paradigm of a smooth surface referred to as skin which is based on Voronoi, Delaunay, and Alpha complexes of a finite set of weighed points. The skin model has good geometric properties that it is free of singularity and it can be decomposed into a collection of quadratic patches. However when it is applied to the energetic computation, the skin triangulation has to be used [Cheng and Shi 2004; Cheng and Shi 2005], which only provides a linear approximation and therefor is dense. The dense triangulation causes oversampling on the surface and reduces the efficiency. Therefore it still remains a challenge to generate a model for the molecular surface which is accurate, smooth, and computable.
(2)
where Γ is the MS of the solute, xi is the center of atom i, and n(r) is the unit normal of the surface at r. The derivation of (2) and an algorithm for the fast evaluation of (2) based on the fast Fourier transform (FFT) is discussed in detail in [Bajaj and Zhao 2006]. The integration in (2) may be more accurately approximated by a Gaussian quadrature over the 3D MS. Since the Gaussian quadrature over a regular 2D planar domain can be easily achieved, this motivates us to define Γ analytically and parameterize Γ into a regular 2D domain.
The main contribution of this paper is a method to model the SES as piecewise algebraic spline patches with certain continuity at the boundary of the patches. Each patch has its dual implicit and parametric representation. Therefore the surface of higher degree is parameterized over a planer domain which makes it convenient to implement numerical quadrature involved in the energetic computation. The idea of constructing a simplical hull over the triangulation and generating the piecewise smooth patches is first introduced in [Dahmen 1989] where a series of tetrahedra are built and the quadric patches are generated. This idea is extended to cubic patches construction in [Guo 1991; Dahmen and Thamm-Schaar 1993], nonsingular and single sheeted cubic patches in [Bajaj et al. 1995]. Here the patches which are built based upon the prism scaffold surrounding the triangulation of the SES are defined implicitly and much more simple by the BB spline functions. Compared with the linear models, fewer number of triangles are needed to get the similar result. Besides, it is error bounded and in certain circumstances it is free of singularity.
Figure 2: The different molecular surfaces are shown for a two atom model in 2D. The VWS is shown as the red curve. The SAS is illustrated as the curve in purple. The SES is the envelope of the pink region. The solvent probe is illustrated as the blue ball. Three molecular surface well known are shown in Figure 2 in 2D. The van der Waals surface (VWS) is the union of a set of spheres with atomic van der Waals radii. The solvent accessible surface (SAS) is the union of a set of spheres with the van der Waals radii augmented by the solvent probe radius (normally taken as ˚ [Lee and Richards 1971]. The solvent excluded surface (SES, 1.4A) also called molecular surface or Connolly surface) is the boundary of the union of all possible solvent probes that do not overlap with the VWS [Connolly 1983; Richards 1977]. In practice SES can be got by starting with SAS and removing the volume occupied by the solvent probe. Following [Connolly 1983], the SES consists of convex spherical patches which are parts of the VWS, concave spherical patches which lie on the probes contacting three atoms, and toroidal patches that are defined by the rolling probe touching two atoms. The VWS and SAS lead to an overestimation and underestimation of the electrostatic solvation energy, respectively [Lee et al. 2003], however SES is advanced to the other two and is most often used as the MS in the energetic calculation. The SES is proper for the stable energy calculation except for one significant drawback: cusps caused by self-intersection of the rolling probe lead to singularity in the Born radii calculation and their differentiation computation as well.
The paper is organized as follows: Section 2 describes in detail the generation of the algebraic spline molecular surface (ASMS). Section 3 discusses the approximation of ASMS, its application to the energetic computation, and provides some examples.
2
Algebraic spline model
2.1
Algorithm Sketch
There are basically four steps in our ASMS construction algorithm: • construct a triangulation mesh of the MS (Section 2.2); • build a prism scaffold based on the triangulation (Section 2.3); • generate piecewise algebraic patches that match with certain continuity (Section 2.3); • parameterize the patches and implement the quadrature over the parametric domain (Section 2.5).
In the energetic computation, an analytical molecular surface model is needed and the singularity should be avoided. One way for generating a smooth model to approximate the SES is to define an analytical volumetric density function, such as the summation of Gaussian functions [Grant and Pickup 1995], Fermi-Dirac switching function [Lee et al. 2002], summation [Lee et al. 2003] or product [Im et al. 2003] of piecewise polynomials, and construct the surface as an iso-contour of the density function. Efficient methods have been developed to construct the iso-contour of the smooth kernel functions [Bajaj et al. 2005; Bajaj and Siddavanahalli 2006b].
2.2
Generation of a triangulation mesh of Γ
So far a lot of work has been done for the triangulation of Γ that approximates the MS [Akkiraju and Edelsbrunner 1996; Laug and Borouchaki 2002; Cheng and Shi 2005; Zhang et al. 2006; Bajaj and Siddavanahalli 2006a]. We could use any of these triangulation meshes as the initial of our algorithm. In the current research we choose to use the triangulation generated by a program embedded
298
in the software TexMol [Bajaj et al. 2004; Bajaj and Siddavanahalli 2006a]. Molecular features are well preserved in this triangulation. We further apply a procedure of triangulation decimation [Bajaj et al. 2002] to obtain a coarser mesh.
To obtain C1 continuity at the vertices, we let b210 − b300 = 1 3 ∇F(vi ) · (v j (λ ) − vi (λ )), where ∇F(vi ) = ni . Therefore
2.3
b120 , b201 , b102 , b021 , b012 are defined similarly.
1 b210 = λ + ni · (v j (λ ) − vi (λ )) 3
Generation of implicit/parametric patches
(5)
We also assume that S is C1 continuous at the midpoints of the edges of T . We define b111 using using the side-vertex scheme [Nielson 1979]:
Given the triangulation mesh T , let [vi v j vk ] be a triangle of T where vi , v j , vk are the vertices of the triangle. Suppose the unit normals of the surface at the vertices are also known, denoted as nl , (l = 1, j, k). Let vl (λ ) = vl + λ nl . First we define a prism Di jk := {p : p = b1 vi (λ ) + b2 v j (λ ) + b3 vk (λ ), λ ∈ Ii jk }, where (b1 , b2 , b3 ) are the barycentric coordinates of points in [vi v j vk ], and Ii jk is a maximal open interval such that 0 ∈ Ii jk and for any λ ∈ Ii jk , vi (λ ), v j (λ ) and vk (λ ) are not collinear and ni , n j and nk point to the same side of the plane Pi jk (λ ) := {p : p = b1 vi (λ ) + b2 v j (λ ) + b3 vk (λ )} (Figure 3).
(1)
(2)
(3)
b111 = w1 b111 + w2 b111 + w3 b111
(6)
where wi =
b2j b2k b22 b23 + b21 b23 + b21 b22
i = 1, 2, 3; i 6= j 6= k
,
(1)
(2)
(3)
Next we are going to explain how b111 , b111 and b111 are defined such that the C1 continuity is obtained at the midpoint of the edge v j vk , vi vk and vi v j . Consider the the edge vi v j in the prism Di jk . Recall that any point p := (x, y, z) in Di jk can be represented by (x, y, z)T = b1 vi (λ ) + b2 v j (λ ) + b3 vk (λ )
(7)
Therefore, after differentiate both sides of (7) with respect to x, y and z respectively, we get Figure 3: A prism Di jk constructed based on the triangle [vi v j vk ].
I3 =
Next we define a function over the prism Di jk in the BensteinBezier (BB) form: F(b1 , b2 , b3 , λ ) =
∑
bi jk (λ )Bnijk (b1 , b2 , b3 )
(3)
∂ b2 ∂x ∂ b2 ∂y ∂ b2 ∂z
∂λ ∂x ∂λ ∂y ∂λ ∂z
(vi (λ ) − vk (λ ))T T (v j (λ ) − vk (λ )) (b1 ni + b2 n j + b3 nk )T
where I3 is a 3 × 3 unit matrix. Let (vi (λ ) − vk (λ ))T T = (v j (λ ) − vk (λ ))T (b1 ni + b2 n j + b3 nk )T
i+ j+k=n
where Bnijk (b1 , b2 , b3 ) is the Bezier basis Bnijk (b1 , b2 , b3 ) =
∂ b1 ∂x ∂ b1 ∂y ∂ b1 ∂z
n! i j k b b b i! j!k! 1 2 3
(8)
(9)
Let A = vi (λ ) − vk (λ ), B = v j (λ ) − vk (λ ) and C = b1 ni + b2 n j + b3 nk , then T T = A, B, C By (8) we have
∂ b2 ∂x ∂ b2 ∂y ∂ b2 ∂z
∂λ ∂x ∂λ ∂y ∂λ ∂z
−1 =T =
1 (B ×C, C × A, A × B) det(T ) (10)
According to (3), at (b1 , b2 , b3 ) = ( 21 , 12 , 0) (the midpoint of vi v j ), we have ∂F (vi (λ ) − vk (λ ))T ∂b ni + n j ∂ F1 T ∂ b2 = (v j (λ ) − vk (λ )) 4 ∂F (ni + n j )T /2 ∂λ 3 2 (b210 − b111 ) 3 + (b120 − b111 )
Figure 4: The control coefficients of the cubic Bezier basis of function F We approximate the molecular surface by the zero contour of F, denoted as S. In order to make S smooth, n, the degree of the Bezier basis, should be no less than 3. To avoid redundancy, we consider the case of n = 3. The control coefficients should be defined properly such that the continuity of S is satisfied. Figure 4 shows the relationship of the control coefficients and the points of the triangle when n = 3. Next we are going to discuss how to define these coefficients.
2
1 2
By (6), at (b1 , b2 , b3 ) = ( 12 , 12 , 0) we have
Since S passes through the vertices vi , v j , vk , we know that b300 = b030 = b003 = λ
∂ b1 ∂x ∂ b1 ∂y ∂ b1 ∂z
(3)
b111 = b111
(4)
299
Therefore the gradient at ( 12 , 12 , 0) is
one can be uniquely map x to a point (b1 , b2 , b3 ) in the jth triangle ([vi v j vk ]), where b3 = 1 − b1 − b2 . Therefore, we have
∂F ∂F ∂F T , , ) ∂ b1 ∂ b2 ∂ λ ni + n j 1 (3) + [3(b210 − b111 )B ×C = 4 2 det(T )
∇F = T −1 (
(3)
+ 3(b120 − b111 )C × A + A × B]
x = x(b1 , b2 , ),
y = y(b1 , b2 ),
z = z(b1 , b2 )
Thus, Z
Z
p f (x(b1 , b2 ), y(b1 , b2 ), z(b1 , b2 )) EG − F 2 db1 db2
f (x) dS =
(11) Γi
σi
(16)
Let
where d1 (λ ) = v j (λ ) − vi (λ ) = B − A d2 (b1 , b2 , b3 ) = b1 ni + b2 n j + b3 nk = C d3 (b1 , b2 , b3 , λ ) = d1 × d2 = B ×C +C × A
∂x 2 ∂y 2 ∂z 2 ) +( ) +( ) ∂ b1 ∂ b1 ∂ b1 ∂x ∂x ∂y ∂y ∂z ∂z F= + + ∂ b1 ∂ b2 ∂ b1 ∂ b2 ∂ b1 ∂ b2 ∂x 2 ∂y 2 ∂z 2 G=( ) +( ) +( ) ∂ b2 ∂ b2 ∂ b2 E =(
(12)
Define 1 1 1 1 c = C( , , 0)d3 (λ ) = d3 ( , , 0, λ ) = B × c + c × A 2 2 2 2
(13)
Let ∇F = ∇F( 12 , 12 , 0). In order to make S be C1 at ( 12 , 12 , 0), we should have ∇F · d3 (λ ) = 0. Therefore, by (11) and (12), we have (3)
b111 =
d3 (λ )T (3b210 B × c + 3b120 c × A + A × B) 3||d3 (λ )||2 (1)
We apply the Gaussian quadrature to (16) Z
(14)
σi
2.4
(2)
3
(15)
Smoothness
Proof: We only need to show S at any point of the edges. Under the condition ni · (vi − v j ) = n j · (v j − vi ), we get ∇F on the edge vi v j is ∇F = ni (1 − t)2 + n j t 2
≤ |λ − λ 0 |(b1 ||ni || + b2 ||n j || + b3 ||nk ||) = |λ − λ 0 | Claim: Let h be the maximum side length of triangulation mesh T , P be the point on the true surface, P0 be the corresponding point on the approximation surface, then P0 converges to P at the rate of O(h3 ). i.e. There exists a constant C such that ||P − P0 || ≤ Ch3 . In order to show the error of S to the true surface, we test some exact surfaces S0 := z = f (x, y), (x, y) ∈ [0, 1]2 with f known. We first generate a triangulation mesh over the true surfaces. Based on the mesh, we construct the piecewise BB zero contour S. The error ||p−p0 || of S compared with the exact surface S0 is defined as max ||p|| , where p ∈ S0 , p0 ∈ S, and the pair of points (p, p0 ) are mapped to the same point (b1 , b2 , b3 ) on the triangle with different λ . With the
then Z j
Error bound and convergence
+ b3 ||vk (λ ) − vk (λ 0 )||
Patch parametrization and quadrature
Γ
Approximation of the ASMS model and an application
||p − p0 || ≤ b1 ||vi (λ ) − vi (λ 0 )|| + b2 ||v j (λ ) − v j (λ 0 )||
Proof of the Theorem 2.3 is similar to the proof of Theorem 2.2.
f (x) dS = ∑
2
Proof:
So S is also C1 on the edges. Theorem 2.3. S is C1 everywhere if ni = n j for every unit normal at the vertices of T .
Z
1
Lemma 3.1. The error of the approximation point p0 to the true point p is bounded by |λ − λ 0 |.
is C1
(r−xi )·n(r) , |r−xi |4
EG − F 2 |bk ,bk
Now we consider the error between S and the true surface. Suppose p := (x, y, z) is a point on the true surface within the prism Di jk . Because of (7), we know the volume coordinate (b1 , b2 , b3 , λ ). The approximation point p0 on S is (b1 , b2 , b3 , λ 0 ), where λ 0 is the solution to the equation F = 0 with (b1 , b2 , b3 ) known.
Proof: It is obvious that S is C1 at the vertices. At the midpoint of edge vi v j , we have ∇F = (ni + n j )/4. Theorem 2.2. S is C1 everywhere if ni · (vi − v j ) = n j · (v j − vi ) for every edge vi v j of T .
To compute (2), let f =
p
k=1
3.1
Theorem 2.1. The ASMS S is C1 at the vertices of T and the midpoints of the edges of T .
2.5
n
∑ Wk f (bk1 , bk2 )
(17) where (bk1 , bk2 ) and Wk are the Gaussian integration nodes and weights for the triangle.
Similarly, we may define b111 and b111 . For the surface evaluation, given the barycentric coordinates of a point (b1 , b2 , b3 ) in triangle [vi v j vk ], we solve equation F = 0 for λ by Newton’s method, where F is defined in (3). Then the corresponding point on S is (x, y, z)T = b1 vi (λ ) + b2 v j (λ ) + b3 vk (λ )
p f (b1 , b2 ) EG − F 2 db1 db2 ≈
f (x) dS Γj
where Γ j is the zero contour of the BB function over the jth triangle. For any point x = (x, y, z) on Γ j , by the inverse map of (15),
300
maximum edge length of the triangles h = 0.1, we sample (p, p0 ) on the surfaces and compute the maximum relative error (Table 1).
Protein ID
To show the convergence of S to S0 , we gradually refine the mess and compute the ratio of the maximum difference of λ and λ 0 to h3 . For each of these functions, as h decreases, we observe that the ratio converges to a constant C (Table 1) This result verifies our claim.
1CGI
Function (x, y) ∈ [0, 1]2 z=0 z = x2 + y2 z = x3 + y3 1 2 2 z = e− 4 [(x−0.5) +(y−0.5) ] cos(5.4y) z = 1.25 + 6+6(3x−1)2 z = tanh(9y − 9x) p z = 1 −p x2 − y2 z = [(2 − 1 − y2 )2 − x2 ]1/2
max{
||p−p0 || ||p|| }
0 2.450030e-05 1.063699e-04 5.286856e-07 2.555683e-04 1.196519e-02 8.614969e-05 1.418242e-05
1HIA
C 0 1.010636e-2 2.610113e-2 6.288604e-5 4.58608e-2 1.896754e-1 1.744051e-1 (h4 ) 1.748754e-02
1PPE
No. of Triangles 29108 8712 3674 27480 7770 3510 24244 6004 2748
G pol , kcal/mol ( timing, s) piecewise linear AS -1371.741894 (39.64) -1343.1496 (40.31) -1399.194841 (12.94) -1346.2230 (12.64) -1678.444735 (7.40) -1394.2270 (6.11) -1361.226603 (30.23) -1340.6384 (31.18) -1389.017538 (9.43) -1347.8067 (9.93) -1571.890827 (5.21) -1388.4665 (5.21) -835.563905 (17.27) -825.3252 (18.26) -852.713039 (5.09) -828.2158 (5.39) -933.956234 (2.74) -845.5085 (3.27)
Table 2: Comparison of the electrostatic solvation energy computed by the piecewise linear surface and the ASMS. With fewer number of triangles used, the ASMS gives a result better than the linear surface which needs more number of triangles. The running time contains the time needed for computing the integration nodes over the surfaces, computing the Born radii, and evaluating G pol .
Table 1: For some typical explicit surfaces (column 1), we compute the maximum relative error (the ratio of the Euclidian distance between p and p0 to the norm of p) in column 2 with h = 0.1, |λ −λ 0 | (x, y) ∈ [0, 1]2 . By letting h ↓ 0, we compute the limit of h3 and observe a constant C (column 3) which verifies the rate of convergence of S we claimed. For certain cases, the rate of convergence can be as good as O(h4 ).
3.2
(b)
(c)
(d)
Biomolecular energetic computation
We apply our ASMS to the GB model to compute the electrostatic solvation energies. As examples, we start from the initial triangulation derived from the PDB (Protein Data Bank) atomic structure for three protein models: 1CGI (852 atoms, Figure 5 (a)), 1HIA (693 atoms, Figure 1 (a)), and 1PPE (436 atoms, Figure 5 (c)). We generate the surface S and compute the electrostatic solvation energy based on S. For the numerical integration (17), 4-point Gaussian quadrature rule over a triangle given in [Dunavant 1985] is applied. We repeat the generation and computation for reduced number of triangles. Figure 1 (d), Figure 5 (b) and Figure 5 (d) show the piecewise algebraic surface generated from the decimated triangulation. Table 2 compares the computed polarization energy with different number of triangles for the proteins. As we see, while a coarser triangulation is used (roughly 1/3 of the original triangulation number), we still get similar results as the initial fine triangulation. We also do the computation on the piecewise linear surface. It turns out that much fewer number of triangles are needed if the piecewise algebraic surface is used.
4
(a)
Figure 5: (a) The atomic model of protein 1CGI (852 atoms); (b) The piecewise algebraic surface of 1CGI (8712 patches); (c) The atomic model of protein 1PPE (436 atoms); (d) The piecewise algebraic surface of 1PPE (6004 patches).
forces in molecular dynamics simulation. The main tasks of the force calculation also involve fast numerical integration. It is challenging because the integration domain is not a surface but a fat skin over the atoms.
Acknowledgment. We wish to thank Vinay Siddavanahalli for helping in the implementation of our ASMS within TexMol, a molecular modeling and visualization software of our center (http://cvcweb.ices.utexas.edu/software/#TexMol). A substantial part of this work in this paper was done when Guoliang Xu was visiting Chandrajit Bajaj at UT-CVC. His visit was supported by the J. T. Oden ICES visitor fellowship.
Conclusions
We have introduced a method to generate a model for the MS. The main steps of the method are the following: construct a triangular prism scaffold covering the PDB structure, and then define nearly C1 piecewise polynomials, BB functions, over the collection of prisms. The model to the molecular surface is the union of the zero contours of the piecewise BB functions. This model is extremely convenient for fast electrostatic solvation energy calculation. We should mention that, while not detailed in this paper, the algorithm of Section 2.3 can, by repeated evocation, yield a hierarchical multiresolution spline model of the MS. In the future research we could extend this algebraic patch model to the computation of electrostatic solvation forces which are the main driven
References A KKIRAJU , N., AND E DELSBRUNNER , H. 1996. Triangulating the surface of a molecule. Discrete Applied Mathematics 71, 5–22.
301
BAJAJ , C., AND S IDDAVANAHALLI , V., 2006. An adaptive grid based method for computing molecular surfaces and properties. ICES Technical Report TR-06-57.
G UO , B. 1991. Modeling arbitrary smooth objects with algebraic surfaces. PhD thesis, Cornell University. I M , W., L EE , M. S., AND B ROOKS , C. L. 2003. Generalized born model with a simple smoothing function. J. Comput Chem. 24, 1691–1702.
BAJAJ , C., AND S IDDAVANAHALLI , V., 2006. Fast error-bounded surfaces and derivatives computation for volumetric particle data. ICES Technical Report TR-06-06.
K ARPLUS , M., AND M C C AMMON , J. A. 2002. Molecular dynamics simulations of biomolecules. Nature Structural Biology 9, 646–652.
BAJAJ , C., AND Z HAO , W. 2006. Fast and accurate generalized born solvation energy computations. Manuscript. BAJAJ , C., C HEN , J., AND X U , G. 1995. Modeling with cubic A-patches. ACM Transactions on Graphics 14, 103–133.
K UHN , B., AND KOLLMAN , P. A. 2000. A ligand that is predicted to bind better to avidin than biotin: insights from computational fluorine scanning. J. Am. Chem. Soc 122, 3909–3916.
BAJAJ , C., L EE , H., M ERKERT, R., AND PASCUCCI , V. 1997. Nurbs based b-rep models from macromolecules and their properties. In Proceedings of Fourth Symposium on Solid Modeling and Applications, 217–228.
L AUG , P., AND B OROUCHAKI , H. 2002. Molecular surface modeling and meshing. Engineering with Computers 18, 199–210. L EE , B., AND R ICHARDS , F. M. 1971. The interpretation of protein structure: estimation of static accessiblilty. J. Mol. Biol. 55, 379–400.
BAJAJ , C., X U , G., H OLT, R., AND N ETRAVALI , A. 2002. Hierarchical multiresolution reconstruction of shell surfaces. Computer Aided Geometric Design 19, 89–112.
L EE , M. S., S ALSBURY, F. R., AND B ROOKS , C. L. 2002. Novel generalized born methods. J Chemical Physics 116, 10606– 10614.
BAJAJ , C., D JEU , P., S IDDAVANAHALLI , V., AND T HANE , A. 2004. Texmol: Interactive visual exploration of large flexible multi-component molecular complexes. Proc. of the Annual IEEE Visualization Conference, 243–250.
L EE , M. S., F EIG , M., S ALSBURY, F. R., AND B ROOKS , C. L. 2003. New analytic approximation to the standard molecular volume definition and its application to generalized born calculations. J. Comput Chem. 24, 1348–1356.
BAJAJ , C., C ASTRILLON -C ANDAS , J., S IDDAVANAHALLI , V., AND X U , Z. 2005. Compressed representations of macromolecular structures and properties. Structure 13, 463–471.
M ADURA , J. D., B RIGGS , J. M., WADE , R. C., DAVIS , M. E., L UTY, B. A., I LIN , A., A NTOSIEWICZ , J., G ILSON , M. K., BAGHERI , B., S COTT, L. R., AND M C C AMMON , J. A. 1995. Electrostatics and diffusion of molecules in solution: simulations with the university of houston brownian dynamics program. Computer Physics Communications 91, 57–95.
BAKER , N., H OLST, M., AND WANG , F. 2000. Adaptive multilevel finite element solution of the poisson-boltzmann equation ii. refinement at solvent-accessible surfaces in biomolecular systems. J. Comput Chem. 21, 1343–1352. BASHFORD , D., AND C ASE , D. A. 2000. Generalized born models of macromolecular solvation effects. Annu. Rev. Phys. Chem. 51, 129–152.
N IELSON , G. 1979. The side-vertex method for interpolation in triangles. J. Apprrox. Theory 25, 318–336.
C HENG , H., AND S HI , X. 2004. Guaranteed quality triangulation of molecular skin surfaces. IEEE Visualization, 481–488.
N INA , M., B EGLOV, D., AND ROUX , B. 1997. Atomic radii for continuum electrostatics calculations based on molecular dynamics free energy simulations. J. Phys. Chem. B 101, 5239– 5248.
C HENG , H., AND S HI , X. 2005. Quality mesh generation for molecular skin surfaces using restricted union of balls. IEEE Visualization, 51–58.
R ICHARDS , F. M. 1977. Areas, volumes, packing, and protein structure. Annu. Rev. Biophys. Bioeng. 6, 151–176.
C ONNOLLY, M. L. 1983. Analytical molecular surface calculation. J. Appl. Cryst 16, 548–558.
ROUX , B., AND S IMONSON , T. 1999. Implicit solvent models. Biophysical Chemistry 78, 1–20.
DAHMEN , W., AND T HAMM -S CHAAR , T.-M. 1993. Cubicoids: modeling and visualization. Computer Aided Geometric Design 10, 89–108.
S CHAEFER , M., AND K ARPLUS , M. 1996. A comprehensive analytical treatment of continuum electrostatics. J. Phys. Chem. 100, 1578–1599.
DAHMEN , W. 1989. Smooth piecewise quadratic surfaces. In Mathematical methods in computer aided geometric design, T. Lyche and L. Schumaker, Eds. Academic Press, Boston, 181– 193.
S RINIVASAN , J., C HEATHAM , T. E., C IEPLAK , P., KOLLMAN , P. A., AND C ASE , D. A. 1998. Continuum solvent studies of the stability of dna, rna, and phosphoramidate-dna helices. J. Am. Chem. Soc 120, 9401–9409.
D UNAVANT, D. 1985. High degree efficient symmetrical gaussian quadrature rules for the triangle. International Journal for Numerical Methods in Engineering 21, 1129–1148.
S TILL , W. C., T EMPCZYK , A., H AWLEY, R. C., AND H EN DRICKSON , T. 1990. Semianalytical treatment of solvation for molecular mechanics and dynamics. J. Am. Chem. Soc 112, 6127–6129.
E DELSBRUNNER , H. 1999. Deformable smooth surface design. Discrete Computational Geometry 21, 87–115.
Z HANG , Y., X U , G., AND BAJAJ , C. 2006. Quality meshing of implicit solvation models of biomolecular structures. Computer Aided Geometric Design 23, 510–530.
G HOSH , A., R APP, C. S., AND F RIESNER , R. A. 1998. Generalized born model based on a surface integral formulation. J. Phys. Chem. B 102, 10983–10990. G RANT, J. A., AND P ICKUP, B. T. 1995. A gaussian description of molecular shape. J. Phys. Chem. 99, 3503–3510.
302
(a)
Guoliang Xu† ICMSEC, Chinese Academy of Sciences, Beijing, 100080, China
(b)
Chandrajit Bajaj ‡ CS and ICES, U. Texas at Austin, Austin, TX 78712
(c)
(d)
Figure 1: Molecular models of a protein(1HIA). (a) The van der Waals surface model (693 atoms); (b) An initial triangulated solvent excluded surface (SES) model (27480 triangles); (c) The decimated triangulated SES model(7770 triangles); (d) Our algebraic spline molecular surface model (7770 patches) generated from (c).
Abstract
1
In this paper, we describe a new method to generate a smooth algebraic spline (AS) model approximation of the molecular surface (MS), based on an initial coarse triangulation derived from the atomic coordinate information of the biomolecule, resident in the PDB (Protein data bank). Our method first constructs a triangular prism scaffold Ps covering the PDB structure, and then generates piecewise polynomial Bernstein-Bezier (BB) spline function approximation F within Ps , which are nearly C1 everywhere. Approximation error and point sampling convergence bounds are also computed. An implicit AS model of the MS which is free of singularity, is extracted as the zero contours of F. Furthermore, we generate a polynomial parametrization of the implicit MS, which allows for an efficient point sampling on the MS, and thereby simplifies the accurate estimation of integrals needed for electrostatic solvation energy calculations.
The computation of electrostatic solvation energy (also known as polarization energy) for biomolecules plays an important role in the molecular dynamics simulation [Karplus and McCammon 2002], the analysis of stability in protein structure prediction [Srinivasan et al. 1998], and the protein-ligand binding energy calculation [Kuhn and Kollman 2000]. The explicit model of the solvent provides the most rigorous solvation energy calculation [Nina et al. 1997]. However, due to the large amount of solvent molecules, most of the computer time is spent on the trajectories of the solvent molecules, which greatly reduces the efficiency of this method [Roux and Simonson 1999]. An alternative method is to represent the solvent implicitly as a dielectric continuum [Schaefer and Karplus 1996] which can be modeled by the Poisson-Boltzmann (PB) equation [Baker et al. 2000; Madura et al. 1995]. A more efficient approximation to the PB electrostatic solvation energy is the generalized Born (GB) model [Still et al. 1990; Bashford and Case 2000; Lee et al. 2003], which approximates the electrostatic solvation energy ∆Gelec as qi q j τ Gpol = − ∑ (1) 2 i, j [r2 + R R exp(− ri2j )] 12
CR Categories: I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling—Curve, surface, solid, and object representations; J.3 [Life and Medical Sciences]: Biology and genetics
Introduction
ij
1 εp
i
j
FRi R j
1 εw ,
where τ = − ε p is the solute (low) dielectric constant, εw is the solvent (high) dielectric constant, qi is the atomic charge of atom i, ri j is the distance between atom i and j, F is an empirical factor (4 [Still et al. 1990], or 8 [Lee et al. 2003]), and Ri is
Keywords: Polynomial splines, molecular surfaces, prismatic scaffolds, Bernstein-Bezier basis, solvation, energetics, error bounds, rate of convergence ∗ Supported in part by NSF grant IIS-0325550, and NIH contract R01GM07308.Email: [email protected] † Support in part by NSFC grant 10371130, National Key Basic Research Project of China (2004CB318000), the J.T. Oden ICES visitor fellowship, and NSF grant IIS-0325550. Email:[email protected] ‡ Supported in part by NSF grants IIS-0325550, CNS0540033 and NIH contracts P20-RR020647, R01-9M074258, R01GM07308.Email:[email protected]
Copyright © 2007 by the Association for Computing Machinery, Inc. Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, to republish, to post on servers, or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from Permissions Dept, ACM Inc., fax +1 (212) 869-0481 or e-mail [email protected]. SPM 2007, Beijing, China, June 04 – 06, 2007. © 2007 ACM 978-1-59593-666-0/07/0006 $5.00
297
the effective Born radius of atom i. The effective Born radius reflects how deeply an atom is buried in the molecule: the deeper it is buried, the larger the radius is and consequently the less important to the polarization. The effective Born radii can be computed by the surface integral [Ghosh et al. 1998] R−1 i =
1 4π
Z Γ
(r − xi ) · n(r) dS |r − xi |4
However the error of the generated isosurface could be very large and result in inaccurate energy computation. Bajaj et al [Bajaj et al. 1997] presented a NURBS representation for the surface. Though it provides a parametric approximation to the SES, it still has singularity over the surface. Edelsbrunner [Edelsbrunner 1999] defines another paradigm of a smooth surface referred to as skin which is based on Voronoi, Delaunay, and Alpha complexes of a finite set of weighed points. The skin model has good geometric properties that it is free of singularity and it can be decomposed into a collection of quadratic patches. However when it is applied to the energetic computation, the skin triangulation has to be used [Cheng and Shi 2004; Cheng and Shi 2005], which only provides a linear approximation and therefor is dense. The dense triangulation causes oversampling on the surface and reduces the efficiency. Therefore it still remains a challenge to generate a model for the molecular surface which is accurate, smooth, and computable.
(2)
where Γ is the MS of the solute, xi is the center of atom i, and n(r) is the unit normal of the surface at r. The derivation of (2) and an algorithm for the fast evaluation of (2) based on the fast Fourier transform (FFT) is discussed in detail in [Bajaj and Zhao 2006]. The integration in (2) may be more accurately approximated by a Gaussian quadrature over the 3D MS. Since the Gaussian quadrature over a regular 2D planar domain can be easily achieved, this motivates us to define Γ analytically and parameterize Γ into a regular 2D domain.
The main contribution of this paper is a method to model the SES as piecewise algebraic spline patches with certain continuity at the boundary of the patches. Each patch has its dual implicit and parametric representation. Therefore the surface of higher degree is parameterized over a planer domain which makes it convenient to implement numerical quadrature involved in the energetic computation. The idea of constructing a simplical hull over the triangulation and generating the piecewise smooth patches is first introduced in [Dahmen 1989] where a series of tetrahedra are built and the quadric patches are generated. This idea is extended to cubic patches construction in [Guo 1991; Dahmen and Thamm-Schaar 1993], nonsingular and single sheeted cubic patches in [Bajaj et al. 1995]. Here the patches which are built based upon the prism scaffold surrounding the triangulation of the SES are defined implicitly and much more simple by the BB spline functions. Compared with the linear models, fewer number of triangles are needed to get the similar result. Besides, it is error bounded and in certain circumstances it is free of singularity.
Figure 2: The different molecular surfaces are shown for a two atom model in 2D. The VWS is shown as the red curve. The SAS is illustrated as the curve in purple. The SES is the envelope of the pink region. The solvent probe is illustrated as the blue ball. Three molecular surface well known are shown in Figure 2 in 2D. The van der Waals surface (VWS) is the union of a set of spheres with atomic van der Waals radii. The solvent accessible surface (SAS) is the union of a set of spheres with the van der Waals radii augmented by the solvent probe radius (normally taken as ˚ [Lee and Richards 1971]. The solvent excluded surface (SES, 1.4A) also called molecular surface or Connolly surface) is the boundary of the union of all possible solvent probes that do not overlap with the VWS [Connolly 1983; Richards 1977]. In practice SES can be got by starting with SAS and removing the volume occupied by the solvent probe. Following [Connolly 1983], the SES consists of convex spherical patches which are parts of the VWS, concave spherical patches which lie on the probes contacting three atoms, and toroidal patches that are defined by the rolling probe touching two atoms. The VWS and SAS lead to an overestimation and underestimation of the electrostatic solvation energy, respectively [Lee et al. 2003], however SES is advanced to the other two and is most often used as the MS in the energetic calculation. The SES is proper for the stable energy calculation except for one significant drawback: cusps caused by self-intersection of the rolling probe lead to singularity in the Born radii calculation and their differentiation computation as well.
The paper is organized as follows: Section 2 describes in detail the generation of the algebraic spline molecular surface (ASMS). Section 3 discusses the approximation of ASMS, its application to the energetic computation, and provides some examples.
2
Algebraic spline model
2.1
Algorithm Sketch
There are basically four steps in our ASMS construction algorithm: • construct a triangulation mesh of the MS (Section 2.2); • build a prism scaffold based on the triangulation (Section 2.3); • generate piecewise algebraic patches that match with certain continuity (Section 2.3); • parameterize the patches and implement the quadrature over the parametric domain (Section 2.5).
In the energetic computation, an analytical molecular surface model is needed and the singularity should be avoided. One way for generating a smooth model to approximate the SES is to define an analytical volumetric density function, such as the summation of Gaussian functions [Grant and Pickup 1995], Fermi-Dirac switching function [Lee et al. 2002], summation [Lee et al. 2003] or product [Im et al. 2003] of piecewise polynomials, and construct the surface as an iso-contour of the density function. Efficient methods have been developed to construct the iso-contour of the smooth kernel functions [Bajaj et al. 2005; Bajaj and Siddavanahalli 2006b].
2.2
Generation of a triangulation mesh of Γ
So far a lot of work has been done for the triangulation of Γ that approximates the MS [Akkiraju and Edelsbrunner 1996; Laug and Borouchaki 2002; Cheng and Shi 2005; Zhang et al. 2006; Bajaj and Siddavanahalli 2006a]. We could use any of these triangulation meshes as the initial of our algorithm. In the current research we choose to use the triangulation generated by a program embedded
298
in the software TexMol [Bajaj et al. 2004; Bajaj and Siddavanahalli 2006a]. Molecular features are well preserved in this triangulation. We further apply a procedure of triangulation decimation [Bajaj et al. 2002] to obtain a coarser mesh.
To obtain C1 continuity at the vertices, we let b210 − b300 = 1 3 ∇F(vi ) · (v j (λ ) − vi (λ )), where ∇F(vi ) = ni . Therefore
2.3
b120 , b201 , b102 , b021 , b012 are defined similarly.
1 b210 = λ + ni · (v j (λ ) − vi (λ )) 3
Generation of implicit/parametric patches
(5)
We also assume that S is C1 continuous at the midpoints of the edges of T . We define b111 using using the side-vertex scheme [Nielson 1979]:
Given the triangulation mesh T , let [vi v j vk ] be a triangle of T where vi , v j , vk are the vertices of the triangle. Suppose the unit normals of the surface at the vertices are also known, denoted as nl , (l = 1, j, k). Let vl (λ ) = vl + λ nl . First we define a prism Di jk := {p : p = b1 vi (λ ) + b2 v j (λ ) + b3 vk (λ ), λ ∈ Ii jk }, where (b1 , b2 , b3 ) are the barycentric coordinates of points in [vi v j vk ], and Ii jk is a maximal open interval such that 0 ∈ Ii jk and for any λ ∈ Ii jk , vi (λ ), v j (λ ) and vk (λ ) are not collinear and ni , n j and nk point to the same side of the plane Pi jk (λ ) := {p : p = b1 vi (λ ) + b2 v j (λ ) + b3 vk (λ )} (Figure 3).
(1)
(2)
(3)
b111 = w1 b111 + w2 b111 + w3 b111
(6)
where wi =
b2j b2k b22 b23 + b21 b23 + b21 b22
i = 1, 2, 3; i 6= j 6= k
,
(1)
(2)
(3)
Next we are going to explain how b111 , b111 and b111 are defined such that the C1 continuity is obtained at the midpoint of the edge v j vk , vi vk and vi v j . Consider the the edge vi v j in the prism Di jk . Recall that any point p := (x, y, z) in Di jk can be represented by (x, y, z)T = b1 vi (λ ) + b2 v j (λ ) + b3 vk (λ )
(7)
Therefore, after differentiate both sides of (7) with respect to x, y and z respectively, we get Figure 3: A prism Di jk constructed based on the triangle [vi v j vk ].
I3 =
Next we define a function over the prism Di jk in the BensteinBezier (BB) form: F(b1 , b2 , b3 , λ ) =
∑
bi jk (λ )Bnijk (b1 , b2 , b3 )
(3)
∂ b2 ∂x ∂ b2 ∂y ∂ b2 ∂z
∂λ ∂x ∂λ ∂y ∂λ ∂z
(vi (λ ) − vk (λ ))T T (v j (λ ) − vk (λ )) (b1 ni + b2 n j + b3 nk )T
where I3 is a 3 × 3 unit matrix. Let (vi (λ ) − vk (λ ))T T = (v j (λ ) − vk (λ ))T (b1 ni + b2 n j + b3 nk )T
i+ j+k=n
where Bnijk (b1 , b2 , b3 ) is the Bezier basis Bnijk (b1 , b2 , b3 ) =
∂ b1 ∂x ∂ b1 ∂y ∂ b1 ∂z
n! i j k b b b i! j!k! 1 2 3
(8)
(9)
Let A = vi (λ ) − vk (λ ), B = v j (λ ) − vk (λ ) and C = b1 ni + b2 n j + b3 nk , then T T = A, B, C By (8) we have
∂ b2 ∂x ∂ b2 ∂y ∂ b2 ∂z
∂λ ∂x ∂λ ∂y ∂λ ∂z
−1 =T =
1 (B ×C, C × A, A × B) det(T ) (10)
According to (3), at (b1 , b2 , b3 ) = ( 21 , 12 , 0) (the midpoint of vi v j ), we have ∂F (vi (λ ) − vk (λ ))T ∂b ni + n j ∂ F1 T ∂ b2 = (v j (λ ) − vk (λ )) 4 ∂F (ni + n j )T /2 ∂λ 3 2 (b210 − b111 ) 3 + (b120 − b111 )
Figure 4: The control coefficients of the cubic Bezier basis of function F We approximate the molecular surface by the zero contour of F, denoted as S. In order to make S smooth, n, the degree of the Bezier basis, should be no less than 3. To avoid redundancy, we consider the case of n = 3. The control coefficients should be defined properly such that the continuity of S is satisfied. Figure 4 shows the relationship of the control coefficients and the points of the triangle when n = 3. Next we are going to discuss how to define these coefficients.
2
1 2
By (6), at (b1 , b2 , b3 ) = ( 12 , 12 , 0) we have
Since S passes through the vertices vi , v j , vk , we know that b300 = b030 = b003 = λ
∂ b1 ∂x ∂ b1 ∂y ∂ b1 ∂z
(3)
b111 = b111
(4)
299
Therefore the gradient at ( 12 , 12 , 0) is
one can be uniquely map x to a point (b1 , b2 , b3 ) in the jth triangle ([vi v j vk ]), where b3 = 1 − b1 − b2 . Therefore, we have
∂F ∂F ∂F T , , ) ∂ b1 ∂ b2 ∂ λ ni + n j 1 (3) + [3(b210 − b111 )B ×C = 4 2 det(T )
∇F = T −1 (
(3)
+ 3(b120 − b111 )C × A + A × B]
x = x(b1 , b2 , ),
y = y(b1 , b2 ),
z = z(b1 , b2 )
Thus, Z
Z
p f (x(b1 , b2 ), y(b1 , b2 ), z(b1 , b2 )) EG − F 2 db1 db2
f (x) dS =
(11) Γi
σi
(16)
Let
where d1 (λ ) = v j (λ ) − vi (λ ) = B − A d2 (b1 , b2 , b3 ) = b1 ni + b2 n j + b3 nk = C d3 (b1 , b2 , b3 , λ ) = d1 × d2 = B ×C +C × A
∂x 2 ∂y 2 ∂z 2 ) +( ) +( ) ∂ b1 ∂ b1 ∂ b1 ∂x ∂x ∂y ∂y ∂z ∂z F= + + ∂ b1 ∂ b2 ∂ b1 ∂ b2 ∂ b1 ∂ b2 ∂x 2 ∂y 2 ∂z 2 G=( ) +( ) +( ) ∂ b2 ∂ b2 ∂ b2 E =(
(12)
Define 1 1 1 1 c = C( , , 0)d3 (λ ) = d3 ( , , 0, λ ) = B × c + c × A 2 2 2 2
(13)
Let ∇F = ∇F( 12 , 12 , 0). In order to make S be C1 at ( 12 , 12 , 0), we should have ∇F · d3 (λ ) = 0. Therefore, by (11) and (12), we have (3)
b111 =
d3 (λ )T (3b210 B × c + 3b120 c × A + A × B) 3||d3 (λ )||2 (1)
We apply the Gaussian quadrature to (16) Z
(14)
σi
2.4
(2)
3
(15)
Smoothness
Proof: We only need to show S at any point of the edges. Under the condition ni · (vi − v j ) = n j · (v j − vi ), we get ∇F on the edge vi v j is ∇F = ni (1 − t)2 + n j t 2
≤ |λ − λ 0 |(b1 ||ni || + b2 ||n j || + b3 ||nk ||) = |λ − λ 0 | Claim: Let h be the maximum side length of triangulation mesh T , P be the point on the true surface, P0 be the corresponding point on the approximation surface, then P0 converges to P at the rate of O(h3 ). i.e. There exists a constant C such that ||P − P0 || ≤ Ch3 . In order to show the error of S to the true surface, we test some exact surfaces S0 := z = f (x, y), (x, y) ∈ [0, 1]2 with f known. We first generate a triangulation mesh over the true surfaces. Based on the mesh, we construct the piecewise BB zero contour S. The error ||p−p0 || of S compared with the exact surface S0 is defined as max ||p|| , where p ∈ S0 , p0 ∈ S, and the pair of points (p, p0 ) are mapped to the same point (b1 , b2 , b3 ) on the triangle with different λ . With the
then Z j
Error bound and convergence
+ b3 ||vk (λ ) − vk (λ 0 )||
Patch parametrization and quadrature
Γ
Approximation of the ASMS model and an application
||p − p0 || ≤ b1 ||vi (λ ) − vi (λ 0 )|| + b2 ||v j (λ ) − v j (λ 0 )||
Proof of the Theorem 2.3 is similar to the proof of Theorem 2.2.
f (x) dS = ∑
2
Proof:
So S is also C1 on the edges. Theorem 2.3. S is C1 everywhere if ni = n j for every unit normal at the vertices of T .
Z
1
Lemma 3.1. The error of the approximation point p0 to the true point p is bounded by |λ − λ 0 |.
is C1
(r−xi )·n(r) , |r−xi |4
EG − F 2 |bk ,bk
Now we consider the error between S and the true surface. Suppose p := (x, y, z) is a point on the true surface within the prism Di jk . Because of (7), we know the volume coordinate (b1 , b2 , b3 , λ ). The approximation point p0 on S is (b1 , b2 , b3 , λ 0 ), where λ 0 is the solution to the equation F = 0 with (b1 , b2 , b3 ) known.
Proof: It is obvious that S is C1 at the vertices. At the midpoint of edge vi v j , we have ∇F = (ni + n j )/4. Theorem 2.2. S is C1 everywhere if ni · (vi − v j ) = n j · (v j − vi ) for every edge vi v j of T .
To compute (2), let f =
p
k=1
3.1
Theorem 2.1. The ASMS S is C1 at the vertices of T and the midpoints of the edges of T .
2.5
n
∑ Wk f (bk1 , bk2 )
(17) where (bk1 , bk2 ) and Wk are the Gaussian integration nodes and weights for the triangle.
Similarly, we may define b111 and b111 . For the surface evaluation, given the barycentric coordinates of a point (b1 , b2 , b3 ) in triangle [vi v j vk ], we solve equation F = 0 for λ by Newton’s method, where F is defined in (3). Then the corresponding point on S is (x, y, z)T = b1 vi (λ ) + b2 v j (λ ) + b3 vk (λ )
p f (b1 , b2 ) EG − F 2 db1 db2 ≈
f (x) dS Γj
where Γ j is the zero contour of the BB function over the jth triangle. For any point x = (x, y, z) on Γ j , by the inverse map of (15),
300
maximum edge length of the triangles h = 0.1, we sample (p, p0 ) on the surfaces and compute the maximum relative error (Table 1).
Protein ID
To show the convergence of S to S0 , we gradually refine the mess and compute the ratio of the maximum difference of λ and λ 0 to h3 . For each of these functions, as h decreases, we observe that the ratio converges to a constant C (Table 1) This result verifies our claim.
1CGI
Function (x, y) ∈ [0, 1]2 z=0 z = x2 + y2 z = x3 + y3 1 2 2 z = e− 4 [(x−0.5) +(y−0.5) ] cos(5.4y) z = 1.25 + 6+6(3x−1)2 z = tanh(9y − 9x) p z = 1 −p x2 − y2 z = [(2 − 1 − y2 )2 − x2 ]1/2
max{
||p−p0 || ||p|| }
0 2.450030e-05 1.063699e-04 5.286856e-07 2.555683e-04 1.196519e-02 8.614969e-05 1.418242e-05
1HIA
C 0 1.010636e-2 2.610113e-2 6.288604e-5 4.58608e-2 1.896754e-1 1.744051e-1 (h4 ) 1.748754e-02
1PPE
No. of Triangles 29108 8712 3674 27480 7770 3510 24244 6004 2748
G pol , kcal/mol ( timing, s) piecewise linear AS -1371.741894 (39.64) -1343.1496 (40.31) -1399.194841 (12.94) -1346.2230 (12.64) -1678.444735 (7.40) -1394.2270 (6.11) -1361.226603 (30.23) -1340.6384 (31.18) -1389.017538 (9.43) -1347.8067 (9.93) -1571.890827 (5.21) -1388.4665 (5.21) -835.563905 (17.27) -825.3252 (18.26) -852.713039 (5.09) -828.2158 (5.39) -933.956234 (2.74) -845.5085 (3.27)
Table 2: Comparison of the electrostatic solvation energy computed by the piecewise linear surface and the ASMS. With fewer number of triangles used, the ASMS gives a result better than the linear surface which needs more number of triangles. The running time contains the time needed for computing the integration nodes over the surfaces, computing the Born radii, and evaluating G pol .
Table 1: For some typical explicit surfaces (column 1), we compute the maximum relative error (the ratio of the Euclidian distance between p and p0 to the norm of p) in column 2 with h = 0.1, |λ −λ 0 | (x, y) ∈ [0, 1]2 . By letting h ↓ 0, we compute the limit of h3 and observe a constant C (column 3) which verifies the rate of convergence of S we claimed. For certain cases, the rate of convergence can be as good as O(h4 ).
3.2
(b)
(c)
(d)
Biomolecular energetic computation
We apply our ASMS to the GB model to compute the electrostatic solvation energies. As examples, we start from the initial triangulation derived from the PDB (Protein Data Bank) atomic structure for three protein models: 1CGI (852 atoms, Figure 5 (a)), 1HIA (693 atoms, Figure 1 (a)), and 1PPE (436 atoms, Figure 5 (c)). We generate the surface S and compute the electrostatic solvation energy based on S. For the numerical integration (17), 4-point Gaussian quadrature rule over a triangle given in [Dunavant 1985] is applied. We repeat the generation and computation for reduced number of triangles. Figure 1 (d), Figure 5 (b) and Figure 5 (d) show the piecewise algebraic surface generated from the decimated triangulation. Table 2 compares the computed polarization energy with different number of triangles for the proteins. As we see, while a coarser triangulation is used (roughly 1/3 of the original triangulation number), we still get similar results as the initial fine triangulation. We also do the computation on the piecewise linear surface. It turns out that much fewer number of triangles are needed if the piecewise algebraic surface is used.
4
(a)
Figure 5: (a) The atomic model of protein 1CGI (852 atoms); (b) The piecewise algebraic surface of 1CGI (8712 patches); (c) The atomic model of protein 1PPE (436 atoms); (d) The piecewise algebraic surface of 1PPE (6004 patches).
forces in molecular dynamics simulation. The main tasks of the force calculation also involve fast numerical integration. It is challenging because the integration domain is not a surface but a fat skin over the atoms.
Acknowledgment. We wish to thank Vinay Siddavanahalli for helping in the implementation of our ASMS within TexMol, a molecular modeling and visualization software of our center (http://cvcweb.ices.utexas.edu/software/#TexMol). A substantial part of this work in this paper was done when Guoliang Xu was visiting Chandrajit Bajaj at UT-CVC. His visit was supported by the J. T. Oden ICES visitor fellowship.
Conclusions
We have introduced a method to generate a model for the MS. The main steps of the method are the following: construct a triangular prism scaffold covering the PDB structure, and then define nearly C1 piecewise polynomials, BB functions, over the collection of prisms. The model to the molecular surface is the union of the zero contours of the piecewise BB functions. This model is extremely convenient for fast electrostatic solvation energy calculation. We should mention that, while not detailed in this paper, the algorithm of Section 2.3 can, by repeated evocation, yield a hierarchical multiresolution spline model of the MS. In the future research we could extend this algebraic patch model to the computation of electrostatic solvation forces which are the main driven
References A KKIRAJU , N., AND E DELSBRUNNER , H. 1996. Triangulating the surface of a molecule. Discrete Applied Mathematics 71, 5–22.
301
BAJAJ , C., AND S IDDAVANAHALLI , V., 2006. An adaptive grid based method for computing molecular surfaces and properties. ICES Technical Report TR-06-57.
G UO , B. 1991. Modeling arbitrary smooth objects with algebraic surfaces. PhD thesis, Cornell University. I M , W., L EE , M. S., AND B ROOKS , C. L. 2003. Generalized born model with a simple smoothing function. J. Comput Chem. 24, 1691–1702.
BAJAJ , C., AND S IDDAVANAHALLI , V., 2006. Fast error-bounded surfaces and derivatives computation for volumetric particle data. ICES Technical Report TR-06-06.
K ARPLUS , M., AND M C C AMMON , J. A. 2002. Molecular dynamics simulations of biomolecules. Nature Structural Biology 9, 646–652.
BAJAJ , C., AND Z HAO , W. 2006. Fast and accurate generalized born solvation energy computations. Manuscript. BAJAJ , C., C HEN , J., AND X U , G. 1995. Modeling with cubic A-patches. ACM Transactions on Graphics 14, 103–133.
K UHN , B., AND KOLLMAN , P. A. 2000. A ligand that is predicted to bind better to avidin than biotin: insights from computational fluorine scanning. J. Am. Chem. Soc 122, 3909–3916.
BAJAJ , C., L EE , H., M ERKERT, R., AND PASCUCCI , V. 1997. Nurbs based b-rep models from macromolecules and their properties. In Proceedings of Fourth Symposium on Solid Modeling and Applications, 217–228.
L AUG , P., AND B OROUCHAKI , H. 2002. Molecular surface modeling and meshing. Engineering with Computers 18, 199–210. L EE , B., AND R ICHARDS , F. M. 1971. The interpretation of protein structure: estimation of static accessiblilty. J. Mol. Biol. 55, 379–400.
BAJAJ , C., X U , G., H OLT, R., AND N ETRAVALI , A. 2002. Hierarchical multiresolution reconstruction of shell surfaces. Computer Aided Geometric Design 19, 89–112.
L EE , M. S., S ALSBURY, F. R., AND B ROOKS , C. L. 2002. Novel generalized born methods. J Chemical Physics 116, 10606– 10614.
BAJAJ , C., D JEU , P., S IDDAVANAHALLI , V., AND T HANE , A. 2004. Texmol: Interactive visual exploration of large flexible multi-component molecular complexes. Proc. of the Annual IEEE Visualization Conference, 243–250.
L EE , M. S., F EIG , M., S ALSBURY, F. R., AND B ROOKS , C. L. 2003. New analytic approximation to the standard molecular volume definition and its application to generalized born calculations. J. Comput Chem. 24, 1348–1356.
BAJAJ , C., C ASTRILLON -C ANDAS , J., S IDDAVANAHALLI , V., AND X U , Z. 2005. Compressed representations of macromolecular structures and properties. Structure 13, 463–471.
M ADURA , J. D., B RIGGS , J. M., WADE , R. C., DAVIS , M. E., L UTY, B. A., I LIN , A., A NTOSIEWICZ , J., G ILSON , M. K., BAGHERI , B., S COTT, L. R., AND M C C AMMON , J. A. 1995. Electrostatics and diffusion of molecules in solution: simulations with the university of houston brownian dynamics program. Computer Physics Communications 91, 57–95.
BAKER , N., H OLST, M., AND WANG , F. 2000. Adaptive multilevel finite element solution of the poisson-boltzmann equation ii. refinement at solvent-accessible surfaces in biomolecular systems. J. Comput Chem. 21, 1343–1352. BASHFORD , D., AND C ASE , D. A. 2000. Generalized born models of macromolecular solvation effects. Annu. Rev. Phys. Chem. 51, 129–152.
N IELSON , G. 1979. The side-vertex method for interpolation in triangles. J. Apprrox. Theory 25, 318–336.
C HENG , H., AND S HI , X. 2004. Guaranteed quality triangulation of molecular skin surfaces. IEEE Visualization, 481–488.
N INA , M., B EGLOV, D., AND ROUX , B. 1997. Atomic radii for continuum electrostatics calculations based on molecular dynamics free energy simulations. J. Phys. Chem. B 101, 5239– 5248.
C HENG , H., AND S HI , X. 2005. Quality mesh generation for molecular skin surfaces using restricted union of balls. IEEE Visualization, 51–58.
R ICHARDS , F. M. 1977. Areas, volumes, packing, and protein structure. Annu. Rev. Biophys. Bioeng. 6, 151–176.
C ONNOLLY, M. L. 1983. Analytical molecular surface calculation. J. Appl. Cryst 16, 548–558.
ROUX , B., AND S IMONSON , T. 1999. Implicit solvent models. Biophysical Chemistry 78, 1–20.
DAHMEN , W., AND T HAMM -S CHAAR , T.-M. 1993. Cubicoids: modeling and visualization. Computer Aided Geometric Design 10, 89–108.
S CHAEFER , M., AND K ARPLUS , M. 1996. A comprehensive analytical treatment of continuum electrostatics. J. Phys. Chem. 100, 1578–1599.
DAHMEN , W. 1989. Smooth piecewise quadratic surfaces. In Mathematical methods in computer aided geometric design, T. Lyche and L. Schumaker, Eds. Academic Press, Boston, 181– 193.
S RINIVASAN , J., C HEATHAM , T. E., C IEPLAK , P., KOLLMAN , P. A., AND C ASE , D. A. 1998. Continuum solvent studies of the stability of dna, rna, and phosphoramidate-dna helices. J. Am. Chem. Soc 120, 9401–9409.
D UNAVANT, D. 1985. High degree efficient symmetrical gaussian quadrature rules for the triangle. International Journal for Numerical Methods in Engineering 21, 1129–1148.
S TILL , W. C., T EMPCZYK , A., H AWLEY, R. C., AND H EN DRICKSON , T. 1990. Semianalytical treatment of solvation for molecular mechanics and dynamics. J. Am. Chem. Soc 112, 6127–6129.
E DELSBRUNNER , H. 1999. Deformable smooth surface design. Discrete Computational Geometry 21, 87–115.
Z HANG , Y., X U , G., AND BAJAJ , C. 2006. Quality meshing of implicit solvation models of biomolecular structures. Computer Aided Geometric Design 23, 510–530.
G HOSH , A., R APP, C. S., AND F RIESNER , R. A. 1998. Generalized born model based on a surface integral formulation. J. Phys. Chem. B 102, 10983–10990. G RANT, J. A., AND P ICKUP, B. T. 1995. A gaussian description of molecular shape. J. Phys. Chem. 99, 3503–3510.
302
