New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
The weightedvolume derivative of a spacefilling diagram

Communicated by Michael Levitt, Stanford University School of Medicine, Stanford, CA (received for review December 12, 2002)
Abstract
Computing the volume occupied by individual atoms in macromolecular structures has been the subject of research for several decades. This interest has grown in the recent years, because weighted volumes are widely used in implicit solvent models. Applications of the latter in molecular mechanics simulations require that the derivatives of these weighted volumes be known. In this article, we give a formula for the volume derivative of a molecule modeled as a spacefilling diagram made up of balls in motion. The formula is given in terms of the weights, radii, and distances between the centers as well as the sizes of the facets of the power diagram restricted to the spacefilling diagram. Special attention is given to the detection and treatment of singularities as well as discontinuities of the derivative.
It is widely believed that characterizing the geometry of a protein molecule is essential for understanding its folding process as well as its interactions with other biomolecules and small ligands. Among all geometric measures, volume is probably the most fundamental property to study. The relevance of calculations of molecular volume is seen in the correlation of volume, or change in volume, with experimental factors such as energetics (1) as well as with physical chemical properties of the molecule such as charge, temperature, and pressure (2, 3). There is a large amount of current research that focuses on the calculations of standard atomic volumes (4–6). These atomic volumes have been considered as parameters that can be used to measure the quality of a threedimensional protein structure (7) as well as probes to quantify structural variations in proteins that occur during evolution (8). Atomic volumes have also been used in packing (9–11) and solvationenergy calculations (12, 13). The latter energy function is central to computational biology, because solvation is an important factor determining the structure and thermodynamics properties of molecules in aqueous solution. All solvent effects on a molecule can be included in an effective solvation potential, W = W_{elec} + W_{np}, in which the first term accounts for the molecule–solvent electrostatics polarization, and the second term accounts for the molecule–solvent van der Waals interactions and for the formation of a cavity in the solvent. W_{elec} is usually described by a continuum model such as the generalized Born model (13, 14). Interestingly, an essential element of such a model is the knowledge of the volume associated with the individual atoms of the solute. Many solvation models describe the nonpolar part of the solvation energy W_{np} as a weighted sum of the solventexposed or accessible surface area of each atom of the solute (15, 16). There is evidence, however, that for small solute the hydrophobic term W_{np} is not proportional to the surface area (17) but rather to the excluded volume of the molecule (18). A volumedependent solvation term was originally introduced by Gibson and Scheraga (19) as the hydration shell model. Biomolecular simulations with implicit solvent that rely on iterative energy calculations such as energy minimization and molecular dynamics require efficient, exact algorithms for the evaluation of the solvationenergy function and its derivatives. In this article, we provide a complete geometric description of the derivatives of the weighted volume of a molecule when its atoms are in motion.
Lee and Richards (20) define the accessible surface of a protein as the van der Waals surface of the molecule expanded by the radius of the solvent sphere about each atom center. The solventexcluded volume is the volume enclosed by the accessible surface. Computational methods that evaluate the solventexcluded volume can be divided into approximate and exact methods. Most of the approximate methods rely on numerical integration (21–23), whereas most analytical methods (9, 12, 24–27) point out that overlap volumes of four or more atoms do not occur in real molecules. An alternative approach for computing the volume and studying the packing of balls was originally developed by Voronoi (28) and later was applied to protein by Richards (9). It has been used successfully for the calculation of protein constituents, the description of protein motions, and the analysis of cavities in proteins (2, 8, 29, 30). A drawback of this approach is that it requires information on the first shell of hydration so that polyhedra may be defined for atoms at the surface of the molecule. It was also developed for equalsize balls and its applications to proteins with atoms that have different radii required approximations (30). Pavani and Ranghino (22) proposed a method for computing the volume of a molecule by using the inclusion–exclusion principle. In their implementation, only intersections of up to three balls were considered. This idea was generalized to any levels of intersections by Gibson and Scheraga (26), applying a theorem from mathematics that higherorder unions and intersections can always be reduced to lowerorder unions and intersections (31). Doing the reduction correctly, however, remains computationally difficult and expensive. The alphashape theory solves this problem exactly by using Delaunay triangulations and their filtrations (32). Alpha shapes have been used to compute the surface area and volume of proteins as well as for detecting and measuring cavities in proteins (33). Fig. 1 provides a simple overview of the application of the alphashape theory to compute the volume of a protein.
The distinction between approximate and exact computation also applies to existing methods for computing the derivatives of the solventexcluded volume with respect to atomic coordinates (34–36). All these methods have to take a large number of singularities into account, where approximations are usually required (36). This problem of singularities is even more acute for surfacearea calculation (37, 38). The alphashape theory proposes a robust solution to this problem by implementing arbitrary precision arithmetic to avoid numerical problems and systematically resolving all singularities without explicitly perturbing the positions of the ball centers. The latter method is referred to as Simulation of Simplicity (39). In this article we describe an extension of the alphashape method that implements the weighted volume derivative theorem described below to provide an efficient, robust, exact computation of the derivatives of volumes. There is an inherent difficulty in using a potential based on volume for energy minimization or molecular dynamics. Although the volume is continuous in the position of the atoms, its derivatives may not be. We examine this issue within the framework of the alphashape method and relate discontinuities to combinatorial changes in the subcomplex of the Delaunay triangulation that is dual to the spacefilling diagram of the molecule.
Background and Results introduces the background and states the main result of this paper. Derivation uses geometric arguments to prove the formula for the weightedvolume derivative. Discontinuities and Implementation discuss the continuity of the derivatives and the implementation of the theorem, and Discussion concludes the article.
Background and Results
As common in biology, we model an atom as a ball and a molecule as the union of a finite collection of such balls. This union is referred to as the spacefilling diagram of the molecule. In this section we explain how we approach the problem of computing the derivatives of the weighted volume of a threedimensional spacefilling diagram.
Geometric Structures.
Let B_{i} = (z_{i}, ϱ_{i}) be the ball with center z_{i} and radius ϱ_{i} for 0 ≤ i < n and call F = ∪_{i} B_{i} a spacefilling diagram. The corresponding state is the point z ∈ ℝ^{3n} that lists the 3n center coordinates in sequence. The weighted volume of F is a map W:ℝ^{3n} → ℝ, and its derivative at z is a linear map DW_{z}:ℝ^{3n} → ℝ defined by DW_{z}(t) = 〈w, t〉, where w ∈ ℝ^{3n} is the gradient of W at z. Specifying the derivative is equivalent to giving the gradient.
Denote the sphere bounding B_{i} by S_{i}. The power distance of a point x ∈ ℝ^{3} from S_{i} is π_{i}(x) = ∥x − z_{i}∥^{2} − ϱ, and the power cell of S_{i} is the set of points P_{i} for which S_{i} minimizes the power distance. We assume or simulate the generic case, in which every power cell is a threedimensional convex polyhedron, and every facet, edge, and vertex belongs to exactly two, three, and four power cells. The power diagram consists of all power cells and their facets, edges, and vertices. As illustrated in Fig. 1, the power diagram decomposes the spacefilling diagram into convex cells of the form P_{i} ∩ F, and these cells share facets, edges, and vertices, which are the intersections of F with the facets, edges, and vertices of the power diagram. The dual complex K of that decomposition consists of all vertices, edges, triangles, and tetrahedra that are convex hulls of sphere centers with cells that intersect in nonempty cells, facets, edges, and vertices.
Sizes and Distances.
We use fractions to express the sizes of geometric entities in the decomposition of the spacefilling diagram. For example, β_{i} = vol(P_{i} ∩ F)/vol(B_{i}) is the fraction of the ith ball that belongs to its power cell. The volume of F = ∪_{i} B_{i} is therefore (4π/3) ∑_{i} β_{i}ϱ, and given real weights α_{i}, the weighted volume is W = (4π/3) ∑_{i} α_{i}β_{i}ϱ. The formula for the derivative of W needs fractions of disks, which we now discuss. Two spheres intersect in a circle S_{ij} = S_{i} ∩ S_{j}, which bounds a disk B_{ij}. The fraction of the disk that belongs to the power diagram is 1 As explained in ref. 32, it is straightforward to compute these fractions from the dual complex K of F by using the principle of inclusion–exclusion. More details on how to perform this calculation are provided in Appendix. We will also need the area of the facet P_{i} ∩ P_{j} ∩ F, which we compute as β_{ij}A_{ij}. An expression for A_{ij} = area(B_{ij}) in terms of radii and distances between centers of spheres will be given in Derivation (see Eq. 7).
The derivative of W also requires the distance between the line L_{ij} passing through the centers z_{i} and z_{j} and the line L_{ijk}, the points of which have equal power distance from S_{i}, S_{j}, and S_{k}. This distance δ_{ij,k} is equal to the distance between the center z_{ij} of S_{ij} and the point y_{ijk} on L_{ijk} closest to z_{ij}. We note that y_{ijk} is the unique point in the plane of z_{i}, z_{j}, and z_{k} that has equal power distance to all three spheres.
Averages.
The derivative of W also refers to the weightedaverage vector from the center to the boundary of a facet P_{i} ∩ P_{j} ∩ F. The weight is the infinitesimal contribution to the area of the facet as we rotate the vector. That weight is constant over each edge and each arc in the boundary of the facet. We therefore compute the weightedaverage vector as the weighted sum of the average vectors of the edges and arcs, with the weights being the areas of the subtended sectors. The average vector from the center to an edge xy = P_{i} ∩ P_{j} ∩ P_{k} ∩ F is the vector to the midpoint, V_{xy} = (x + y)/2 − z_{ij}. Consider next an arc ab on the circle S_{ij}, and let ϕ = ϕ_{ab} be half the angle, as in Fig. 2. Assuming the center is the origin, the radius is 1, and the arc midpoint lies on the abscissa, we get the average by integration 2 To give a formula without the normalizing assumptions, we use the midpoint M_{ab} of the arc ab. The average vector is V_{ab} = v_{ab}(M_{ab} − z_{ij}).
In general, the facet P_{i} ∩ P_{j} ∩ F has a boundary that consists of several edges and arcs. The weightedaverage vector from the center to that boundary is the weighted sum of the individual averages, 3 where h_{xy} = δ_{ij,k} is the distance from z_{ij} to the line spanned by xy, and ϱ_{ij} is the radius of S_{ij}. The first sums range over all edges xy, and the second sums range over all arcs ab of the facet. Note that the denominator is twice the area of the facet. We provide a method to compute the weightedaverage vector using inclusion–exclusion in Appendix.
Main Result.
The main result of this article is a complete description of the weightedvolume derivative. We need some notation. The weight of the ith ball is α_{i} ∈ ℝ, which may be positive or negative. We write ζ_{ij} = ∥z_{i} − z_{j}∥ for the distance between two centers and U_{ij} = (z_{i} − z_{j})/ζ_{ij} for the unit vector between the same. Recall also the definitions of β_{ij}, A_{ij}, and V_{ij} given above.
WeightedVolume Derivative Theorem.
The derivative of the weighted volume of the spacefilling diagram with state z is DW_{z}(t) = 〈w, t〉, where 4for 0 ≤ i < n. The sum is over all edges z_{i}z_{j} in K.
We illustrate the formula with a small example consisting of three spheres that intersect in two points, the projection of which along the shared edge of the power diagram is shown in Fig. 3. In the unweighted case we have α_{i} = 1, w_{ij} = 1, and x_{ij} = 0 for all i and j. The derivative thus simplifies to a sum of vectors between the centers: ∑ β_{ij}A_{ij}U_{ij}. As shown by Csikós (40), this formula generalizes to dimensions different from three. We note that DW is not everywhere continuous. Specifically, it may not be continuous when two or more spheres coincide or three or more spheres intersect in a common circle (see Discontinuities).
Derivation
In this section we derive the formula claimed in the WeightedVolume Derivative Theorem. It is the sum of contributions from locally directionpreserving and locally distancepreserving components of the motion.
DirectionPreserving Motion.
A direction is determined by two spheres and is preserved if one sphere stays put and the other moves along that line. To study the volume derivative for such a motion, we restrict out attention to two spheres S_{i} and S_{j} with nonempty intersection S_{ij} and define F = B_{i} ∪ B_{j}. Let ζ_{i} and ζ_{j} be the signed distances of the centers from the plane of S_{ij}, as illustrated in Fig. 4. The sum of the two distances is ζ_{ij} = ζ_{i} + ζ_{j} = ∥z_{i} − z_{j}∥. We have ϱ − ϱ = ζ − ζ = ζ_{ij}(ζ_{i} − ζ_{j}), and therefore 5 and 6 It is geometrically obvious that the weightedvolume derivative is the area of the disk bounded by S_{ij} times the sum of the weighted derivatives of ζ_{i} and ζ_{j}. Because the square radius of the disk is ϱ = ϱ − ζ, the area of the disk is 7 where the product ranges over the four different ways of assigning signs. The derivatives of ζ_{i} and ζ_{j} with respect to ζ_{ij} can be computed from Eqs. 5 and 6, and the weightedvolume derivative is 8 Alternatively, we could work analytically, compute the weightedvolume contributions of the two balls by integration, and get dW/dζ_{ij} by differentiation. The result is the same.
DistancePreserving Motion.
The distance between the centers of two spheres S_{i} and S_{j} is preserved if one rotates about the other. We consider what happens to the facet P_{i} ∩ P_{j} ∩ F that separates the two corresponding cells. As illustrated in Fig. 5, that facet is the intersection of a disk and a convex polygon. The boundary consists of edges and circular arcs. Let T_{i} normal to U_{ij} be the motion vector of S_{i}, and let R_{i,j} be the line normal to both vectors that passes through the center z_{ij} of the disk. The infinitesimal motion of the facet consists of a rotation about R_{i,j} composed with a translation along T_{i}, but only the rotation has a nonzero contribution to the volume derivative. On one side of the line we trade volume of B_{i} for that of B_{j}, and on the other side we do it the other way round. The net gain or loss is the accumulation of the contributions of the sectors over the edges and arcs in the boundary of the facet. Here we count a sector over an edge negative if the line defined by the edge separates z_{ij} and the facet.
We first consider an edge xy defined by S_{i}, S_{j}, and a third sphere S_{k}. We compute the derivative of the signed volume swept out by the triangle z_{ij}xy with respect to the rotation angle θ, 9 where h_{xy} = δ_{ij,k} is the distance from z_{ij} to the line passing through x and y, and V_{xy} = (x + y)/2 − z_{ij} is the vector to the midpoint. Note that the result is two thirds the area of the triangle times the scalar product divided by the distance between the two ball centers.
Second, we consider a circular arc ab, which is part of the circle S_{ij}. We compute the derivative of the volume swept out by the disk sector z_{ij}ab with respect to θ, 10 where 2ϕ_{ab} is the angle of the arc, and V_{ab} is the average vector from the center to the arc. Again, the result is two thirds the area of the sector times the scalar product divided by the distance between the two ball centers. As explained in Background and Results, V_{ab} = v_{ab}(M_{ab} − z_{ij}), where M_{ab} is the midpoint of the arc and v_{ab} = (sin ϕ_{ab})/ϕ_{ab}.
Assembling the Relations.
Let t be the velocity vector of the state z and T_{i} = [t_{3i+1}, t_{3i+2}, t_{3i+3}]^{T} the velocity vector of z_{i}. Because the derivative is linear, we can decompose the motion into components and add the contributions. For every ordered pair i, j, we consider the directionpreserving component 〈U_{ij}, T_{i}〉U_{ij} and the distancepreserving component T_{i} − 〈U_{ij}, T_{i}〉U_{ij} of T_{i}.
The contribution of the directionpreserving component is the fraction of the disk B_{ij} that belongs to the power diagram times what is given in Eq. 8. This is β_{ij}A_{ij}w_{ij}〈U_{ij}, T_{i}〉, where w_{ij} is given in WeightedVolume Derivative Theorem. The contribution of the distancepreserving component is the weight difference times two thirds the area times the weighted average vector divided by the distance between the ball centers. The factor 2/3 is necessary to change the factor 1/2, which is used in the area computation of sectors, to 1/3, which arises in the volume computation of cones. The contribution is then x_{ij}β_{ij}A_{ij}〈V_{ij}, T_{i}〉, where x_{ij} is given in WeightedVolume Derivative Theorem. It is also the weight difference times the sum of contributions over all edges, given in Eq. 9, plus the sum of contributions over all arcs, given in Eq. 10. Adding the terms for all ordered pairs gives 11 for the weightedvolume derivative. We can write this more succinctly as DW_{z}(t) = 〈w, t〉, where w is as defined in WeightedVolume Derivative Theorem. We have β_{ij} = 0 unless z_{i}z_{j} is an edge of K, which implies that the summation can be limited to the edges in K.
Discontinuities
The weightedvolume derivative, DW, maps a state z to the gradient w = w_{z} of W at z. From the formula in WeightedVolume Derivative Theorem, we glean that DW is continuous except at states at which either two or more spheres coincide or three or more spheres intersect in a common circle. The set of such states is a (3n − 3)dimensional subset of ℝ^{3n}, which should be compared with the (3n − 1)dimensional subset at which the weightedarea derivative is discontinuous. To see that there are only these two types of discontinuities, we note that A_{ij} is everywhere continuous, x_{ij}, w_{ij}, and U_{ij} are continuous except when z_{i} = z_{j}, and β_{ij} and V_{ij} are continuous as long as the facet P_{i} ∩ P_{j} ∩ F varies continuously. States at which that facet changes abruptly either have z_{i} = z_{j} or contain at least one sphere, other than S_{i} and S_{j}, that contains the circle S_{ij}.
Every state z at which DW is discontinuous has the property that every open neighborhood contains states with combinatorially different dual complexes. In other words, when z passes through a discontinuity of DW, then the dual complex passes a moment at which it may undergo a combinatorial change. Note, however, that not every combinatorial change in the dual complex corresponds to a discontinuity in DW.
Implementation
To compute the volume and its derivatives with respect to the coordinates of the center of the balls, we first need the dual complex of the spacefilling diagram. For that purpose, we have written a new version of the alphashape software, alphavol, specific to molecular simulation applications. alphavol first performs a regular triangulation of the set of centers of the balls by using the algorithm described by Edelsbrunner and Shah (41). Singularities are resolved by using the Simulation of Simplicity (39), whereas numerical imprecisions are avoided by using arbitrary precision arithmetic. The software is kept efficient through the application of a floatingpoint filter. The dual complex is computed by using the same care for singularities and numerical uncertainties. alphavol includes implementations of the volume and weightedvolume derivative formulas. Computation of the volume and its derivatives of a 250residue protein by using alphavol requires ≈0.37 sec on a 1,000Mhz Pentium processor. The 0.37 sec roughly breaks down to 0.11 sec for computing the regular triangulation, 0.03 sec for generating the dual complex, 0.17 sec to compute the weighted volume, and 0.06 sec to compute the weightedvolume derivatives.
The weightedvolume derivative formula, Eq. 4, was tested against the numerically estimated derivatives, 12 where z_{3i+j} is the jth coordinate of the ith atom center. The average relative difference between analytical and numerical derivatives was computed as a relative root mean square error: 13 Using δ = 0.0001, we find μ(W) equal on the average to 5.1 × 10^{−8} and never larger than 9 × 10^{−8} for a collection of 100 proteins varying in size from 50 residues or ≈400 atoms to 500 residues or ≈4,000 atoms. The good agreement between the analytical and numerical derivatives of the weighted volume shows that our approach and implementation are correct. alphavol is available online at http://biogeometry.duke.edu/software/proshape.
Discussion
The software alphavol, which computes the weighted volume of a union of balls and its derivatives with respect to the coordinates of the centers of the balls, provides a fast, accurate, and robust method for computing the interaction of water with nonpolar atoms of a molecule in a hydration shell model (19). ALPHAVOL implements analytical formulas for computing the weighted volume and its derivatives and explicitly deals with the problem of discontinuities. We have made preliminary steps toward inserting alphavol into the molecular dynamics software encad (42) and gromacs (43), but it is too early to say anything about the corresponding results.
Measuring the volume occupied by individual atoms in macromolecular structures has been the subject of scientific investigation for several decades. In recent years, this interest has increased, because volumes of specific atom types are used in a wide range of implicit solvent models. One of the most widespread uses of atomic volumes in molecular mechanics calculations is in the generalized Born model of electrostatic solvation (14), where the atom volume appears in the calculation of the desolvation of one atom by another atom. In the present implementation of the generalized Born model, ball overlaps are ignored, corrected by scaling factors, or dealt with by using a Gaussian representation of the ball. We expect that the formalism introduced in this article for computing volumes and their derivatives that explicitly and accurately deal with ball overlaps should prove useful for the generalized Born model.
Voronoi polyhedra have been widely used to compute the volume occupied by the atoms of proteins (4, 6–9) and recently DNA (5) in solution. One difficulty in applying the Voronoipolyhedra approach to macromolecular structures is that the treatment of atoms at the surface requires special attention because their power cells are infinite. Possible solutions to this problem are the construction of a layer of solvent molecules around the macromolecule (44) or limiting the power diagram to the spacefilling diagram (29, 33). This article describes the latter approach and extends it by providing exact expression for the derivatives of the atomic volume. These derivatives are expected to be useful for including potential based on packing in molecular mechanics or molecular dynamics calculation.
The mathematical tools used in this paper to express the volume derivative of a spacefilling diagram have also been used to compute the weighted surfacearea derivatives (R. Bryant, H.E., P.K., and M. Levitt, unpublished data).
Acknowledgments
This work was supported in part by National Science Foundation Grant CCR0086013. P.K. acknowledges support from the Department of Energy (DEF60395ER62135).
Appendix
In this section we explain how to compute the fractions β_{ij} and the weightedaverage vectors V_{ij} needed to express the weightedvolume derivative. The descriptions are relevant only in the context of implementing the derivative once the dual complex K of the spacefilling diagram is known. We begin by introducing a notion of neighborhood within K: The star of a simplex τ ∈ K is the collection of simplices that contain τ, including τ itself, and the link of τ is the collection of faces of simplices in the star that are disjoint of τ. Writing the face relation between simplices with a smaller equal symbol, we can express these definitions more formally: 14 15 For example, the star of an interior edge contains the edge together with a ring of triangles and tetrahedra around the edge. The link is the cycle of vertices and edges that surrounds the ring away from the interior edge; for each triangle in the star it contains the opposite vertex, and for each tetrahedron it contains the opposite edge.
Computing β_{ij}.
Recall that B_{ij} is the disk spanned by the circle S_{ij}, which is the intersection of the spheres S_{i} and S_{j}. Formally, B_{ij} is the set of points x ∈ ℝ^{3} for which π_{i}(x) = π_{j}(x) ≤ 0. The number we are after is the fraction of that disk that belongs to the power diagram, 16 Given the dual complex K, we can use inclusion–exclusion to compute β_{ij} from simple pieces of the disk that are dual to simplices in the star of the edge z_{i}z_{j}. Let B be the segment of the disk on z_{k}'s side of the bisector. More formally, B is the set of points x ∈ B_{ij} for which π_{k}(x) ≤ π_{i}(x) = π_{j}(x). The general results in ref. 32 imply that β_{ij} can be computed by subtracting the segments B from the disk B_{ij} and adding back their pairwise intersections, 17 where ρ_{ij} is the radius of B_{ij}. The first sum ranges over all vertices z_{k}, and the second ranges over all edges z_{k}z_{l} in the link of z_{i}z_{j} in K.
Computing V_{ij}.
Recall that V_{ij} is the weightedaverage vector from the center z_{ij} of the disk B_{ij} to the boundary of the facet that is the portion of the disk belonging to the Voronoi diagram. The weight of each vector is its infinitesimal contribution to the area of the facet, which is constant along each edge and each arc. Hence, V_{ij} can be computed as the weighted sum of the average vectors of the edges xy and the arcs ab of the facet, 18 where h_{xy} is the distance from z_{ij} to the line spanned by the edge xy, and 2ϕ_{ab} is the angle of the arc ab. The first sums in the numerator and the denominator range over all edges xy, and the second sums range over all arcs ab of the facet. We note that the denominator is twice the area of the facet, which is 2β_{ij}A_{ij}. We thus focus on the numerator, which we will express as an alternating sum using again inclusion–exclusion over the link of the edge z_{i}z_{j} in K. We start with the average vector to the boundary of the entire disk, which is N_{ij} = 0. From this we subtract the contributions of disk segments and we add back contributions of pairs of disk segments.
Case 1 (disk segments).
Let x and y be the points at which the sphere S_{k} intersects the circle S_{ij}. The disk segment B is thus bounded by the edge xy and the arc ab connecting a = x with b = y along the circle as illustrated in Fig. 6 Left. We add the contribution of the edge and subtract the contribution of the arc. The former is twice the area of the triangle times the vector to the midpoint, which is 2ϱ cos ϕ_{ab}sin ϕ_{ab}V_{xy} = ϱ sin(2ϕ_{ab})V_{xy}. Assuming ϕ_{ab} ≠ ±90°, the average vector of the arc can be computed as V_{ab} = V_{xy}(tan ϕ_{ab})/ϕ_{ab}. This expression is undefined if ϕ_{ab} = ±90°, in which case we may go back to Eq. 2 and get V_{ab} = (2/π)(M_{ab} − z_{ij}), where M_{ab} is the midpoint of the arc. The weight of this vector is twice the area of the disk sector defined by ab and is given by 2ϱϕ_{ab}. The accumulated contribution of the disk segment therefore is 19 where we take the contribution of the arc positive and that of the edge negative, in preparation for subtracting N in the final sum.
Case 2 (pairs of disk segments).
Let x, y and u, v be the endpoints of the edges and arcs bounding the disk segments B and B, and define B = B ∩ B. We only need to consider the case in which the endpoints alternate along the circle. As shown in Fig. 6 Right, the intersection point of the two edges lies inside the disk, and the boundary of B consists of one arc and two edges that both end at that intersection point. Let N be the weighted sum of average vectors taking the contribution of the arc positive and those of the edges negative, as before.
We get the weighted average vector over the boundary of the facet P_{i} ∩ P_{j} ∩ F as 20 where N_{ij} = 0, the first sum ranges over all vertices z_{k}, and the second sum ranges over all edges z_{k}z_{l} in the link of z_{i}z_{j} in K.
 Received December 12, 2002.
 Accepted December 22, 2002.
 Copyright © 2003, The National Academy of Sciences
References
 ↵
 ↵
 ↵
 ↵
 ↵
 Nadassy K,
 TomásOliveira I,
 Alberts I,
 Janin J,
 Wodak S J
 ↵
 Tsai J,
 Gerstein M
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Ooi T,
 Oobatake M,
 Nemethy G,
 Scheraga H
 ↵
 ↵
 ↵
 Gibson K,
 Scheraga H
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Voronoi G F
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Gogonea V,
 Osawa E
 ↵
 ↵
 ↵
 Edelsbrunner H,
 Mücke E P
 ↵
 ↵
 ↵
 ↵
 Lindahl E,
 Hess B,
 van der Spoel D
 ↵