- Email: [email protected]
Triangulation of molecular surfaces based on extracting surface atoms
Computers & Graphics 38 (2014) 291–299
Contents lists available at ScienceDirect
Computers & Graphics journal homepage: www.elsevier.com/locate/cag
Special Section on CAD/Graphics 2013
Triangulation of molecular surfaces based on extracting surface atoms Jingqiao Zhang, Zhe Shi n School of Computer Engineering and Science of Shanghai University, Shanghai 200444, China
art ic l e i nf o
a b s t r a c t
Article history: Received 5 August 2013 Received in revised form 7 November 2013 Accepted 7 November 2013 Available online 22 November 2013
In this paper, we propose a novel surface atom extraction algorithm. By calculating the relative positions of the extended spheres, this algorithm detects all surface atoms reachable by the probe. Our algorithm overcomes the drawbacks of current methods and can be implemented in real time. This algorithm significantly improved the efficiency of calculating the solvent-excluded surface because only the values for the surface atoms are calculated instead of for all atoms. We also construct a triangular mesh model of the molecular solvent-excluded surface in terms of the surface atoms. By subdividing the saddle patches into more uniform triangular meshes, we improved the existing methods for triangulating a molecular surface composed of triangular concave spherical patches, saddle patches and convex patches. In addition, we introduce a new method that seamlessly joins saddle patches to convex patches. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Solvent-excluded surface Extended sphere Surface atom Triangulation
1. Introduction The three-dimensional structure of a molecule is one of the most important factors determining its function. The molecular surface is a primary part of the protein structure and is crucial to protein docking and visualization. Currently, most research has focused on three molecular surface models [1–6]: the van der Waals surfaces, solvent-accessible surfaces (SASs) and solventexcluded surfaces (SESs). Solvent-excluded surfaces are also commonly known as solvent contact surfaces (SCSs), Connolly surfaces or smooth molecular surfaces. We consider the SES in this work. The SES was first proposed by Richards [2] in 1977, and it considers a probe sphere representing a solvent molecule rolling over the atoms in a molecule. SESs, which have become a wellknown definition for the molecular surface, consist of three types of patches according to Connolly [7], namely, convex spherical patches, saddle patches and triangular concave spherical patches, as shown in Fig. 1. Throughout the rest of this article the term molecular surface refers to SES unless otherwise indicated. In [7,8], Connolly showed how to calculate the three types of patches as well as the molecular area and volume and briefly discussed triangulating the molecular surface. Varshney et al. [9,10] presented a parallelizable algorithm for computing the molecular surface, which needs O(k log k) over n processors, where n is the number of atoms in the molecule and k is the average number of neighbors per atom. Note that the extended
n
Corresponding author. Tel.: þ 86 18817350768. E-mail addresses: [email protected] (J. Zhang), [email protected] (Z. Shi).
0097-8493/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cag.2013.11.007
sphere defined by Varshney et al. [9,10] is important to our algorithm for extracting molecular surface atoms. In 1996, Sanner et al. [11] proposed a reduced surface consisting of a polyhedron with faces, edges and vertices that indicate whether the probe is in contact with more than two atoms, two atoms or only one atom. This reduced surface can be computed for a molecule from a binary spatial division tree, which provides an efficient algorithm for computing the molecular surface. Sanner et al. [11] claimed that the binary spatial division tree can be constructed using O(n log n) computations in the worst case and the reduced surface can be computed using O(n log n) computations on average, where n is the number of atoms in the molecule. In 2006, Ryu et al. [12,13] introduced a β-shape, which is a generalization of the well-known α-shape [14,15] that can be used to efficiently calculate molecular surfaces. Later, Ryu et al. [16,17] discussed the uniform triangulation of the three patch types. However, constructing the β-shape is difficult. First, for a molecule with n atoms, the Voronoi diagram of atoms requires O(n3) calculations in the worst case (the average complexity is O(n)). Second, for a given Voronoi diagram, the β-shape can be computed in O(log n þk) calculations in the worst case, where k is the number of edges and faces in the resulting β-shape [13,18,19]. Nevertheless, the Voronoi diagram of the atoms and β-shape are powerful tools for other important applications such as protein docking problems. The reduced surface [11] is equivalent to one β-shape. In this paper, we improve Ryu's methods for triangulating saddle patches [17]. Our improved method generates a uniform mesh even though the saddle patch is narrow. Several other researchers have studied the calculation and visualization of various molecular surface models in other ways [20–29].
292
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
Theorem 1. The extended spheres of two atoms intersect (except when tangent) if and only if the two atoms are neighbors. Proof. The extended spheres for two given atoms, Ai ðc i ; r i Þ and Aj ðc j ; r j Þ, are Si ðci ; r i þ r p Þ and Sj ðc j ; r j þ r p Þ, respectively, and the distance between their center is D. If these two extended spheres intersect (except when tangent), we obtain D o L ¼ r i þ r p þ r j þ r p . According to Definition 2, Ai and Aj are neighboring atoms. In contrast, two neighbor atoms have D o L ¼ r i þ r p þ r j þr p ; namely, their extended spheres intersect (except when tangent).
Fig. 1. The solvent-excluded surface.
Most previous works researched various molecular surface models based on all of the atoms in the molecule. However, if the surface atoms can be effectively extracted, the scale of the problem could decrease. Because the SES is constructed independently of the internal atoms, the calculation efficiency for SES can be improved by only accounting for the surface atoms. In this paper, we define surface atoms as those that can be reached by a probe rolling over the molecule and internal atoms as those that cannot be reached. For this surface atom extraction, Wang et al. [30] proposed a geometric algorithm named the Protein Surface Atom Geometrical Algorithm (PSAGA) in 2007 to determine the molecular surface atoms according to a molecule's 3D geometric characteristics. However, the results extracted by this algorithm are not exact. Moreover, it requires the radii of all the atoms to be the same. In 2011, Zhang et al. [31] introduced a new extraction algorithm called the Molecular Surface Atom Rolling Algorithm (MSARA), which is more accurate than the PSAGA while being effective for atoms with differing radius. However, this algorithm cannot find certain surface atoms, and its efficiency is not ideal. Here, we present a new, real-time extraction algorithm that overcomes the drawbacks of the existing methods. Our algorithm consists of three steps: (1) extract the surface atoms; (2) calculate and triangulate the three molecule surface patch types; and (3) connect these patches seamlessly. This paper is organized as follows. Section 2 introduces the new extraction algorithm. Section 3 provides the calculations for the SES based on the extracted surface atoms and discusses the triangulation and connection of the three SES patches. The time complexity is analyzed in Section 4. We test our algorithms in Section 5. Conclusions and future work are presented in Section 6.
Theorem 2. If the surface of an extended sphere for one atom in a molecule is completely covered by the interior of other atoms' extended spheres, it is an internal atom; otherwise, it is a surface atom. Proof. If the surface of an atom's extended sphere is not completely covered, the probe can reach this atom without colliding with other atoms whenever the probe's center lies on the uncovered portion of the atom's extended sphere surface. Thus, the atom is a surface atom. In contrast, if the probe cannot reach this atom without colliding with other atoms, it is an internal atom. Theorem 2 provides a method for determining whether an atom is a surface atom. Fig. 2(a) and (b) shows a molecule and the extended spheres of its atoms, respectively. Because the extended sphere surface illustrated by the red dashed circle in Fig. 2(c) is completely covered, the corresponding atom in Fig. 2(d) is internal. On the other hand, because the extended sphere surfaces, which are illustrated by the black solid circles in Fig. 2(c), are not completely covered, their corresponding atoms in Fig. 2(d) are surface atoms. Moreover, according to Theorem 1, whether an atom is a surface atom only depends on it and its neighbors instead of all atoms. We propose the following method to simplify determining whether the surface of an extended sphere is completely covered.
2. Surface Atom Extraction based on Extended Spheres (SAEES) In the rest of this article, we denote r p as the radius of the probe. Definition 1. Let Aðc ; rÞ be an atom, where c is the center of the atom and r is the radius of the atom. The extended sphere of an atom is a sphere Sðc ; r þ r p Þ with the center at c and of the radius ðr þ r p Þ. The definition of an extended sphere indicates that the center of the probe cannot be inside of an atom's extended sphere; otherwise a collision would occur between the probe and the atom. Definition 2. Given two atoms, Ai ðci ; r i Þ and Aj ðc j ; r j Þ, let D be the Euclidean distance between c i and c j and L ¼ r i þ r j þ 2r p . Ai and Aj are neighbors to each other if D o L.
Fig. 2. Determine whether an atom is a surface atom based on the extended spheres. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.). (a) A molecule, (b) The molecule and the extended spheres of all the atoms, (c) All the extended spheres and (d) Surface atoms and an internal atom.
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
Let Await be an atom with an extended sphere of Swait . A1 with the extended sphere S1 neighbors Await . According to Theorem 1, Swait and S1 intersect. As shown in Fig. 3(a), the intersection between Swait 's surface and S1 's surface is a circle named the intersection circle. The planar point set surrounded by the intersection circle is called the cutting face and is denoted with C. P wait and P 1 are the centers of Swait and S1 , respectively. The vector v is the unit vector of P wait P 1 and can be used as the unit normal vector of C. Imagine Swait is a watermelon with a thin peel; the cutting face is used to cut off part of Swait and v points towards the removed parts, as shown in Fig. 3(b). Suppose Await has m neighboring atoms. There are m cutting faces, denoted by C1, C2… Cm, formed by Swait and m extended spheres from Await 's neighboring atoms. These m cutting faces are stored in a cutting-face list, with each node storing the cutting face's center and a unit normal vector pointing to the removed portion of Swait . Thus, Swait is a watermelon cut by Ci(i¼ 1,…,m). After this “cutting”, if some peel of the Swait remains, that is, the surface of Swait is not completely covered, Await is a surface atom. This method works even when the atomic radii are different and can identify the surface atoms at the boundary of internal voids. Before “cutting” Swait , a preprocess is implemented to reduce the number of cutting faces or in some cases directly determine whether Await is a surface atom. During preprocessing, the intersection line L for the plane of any two cutting faces, C1 and C2, (if they are not parallel) is calculated. If the radius of Swait is longer than the distance from the center of Swait to L, C1 and C2 intersect. Fig. 4(a) shows that two cutting faces intersect when both of them remain.
293
If C1 and C2 do not intersect, as shown in Fig. 4(b)–(d), whether one cutting face (e.g., C1) lies in the abandoned part of Swait due to another cutting face (e.g., C2) should be considered and is easily judged by computing the scalar product of the unit normal vector for C2 and the vector from the center of C2 to the center of C1. There are three cases for the relative positions of C1 and C2: (1) when C1 (C2) lies in the abandoned part of Swait due to C2 (C1) while C2 (C1) does not lie in the abandoned part of Swait due to C1 (C2), C1 (C2) can be deleted, see Fig. 4(b); (2) when neither C1 nor C2 lie in the abandoned part of Swait due to the other, then both remain, see Fig. 4(c); (3) when both lie in the abandoned part of Swait due to the other, Await can be directly determined as being an internal atom, see Fig. 4(d). After reducing the number of cutting faces via this preprocessing, we must traverse each remaining pair of cutting faces. Whenever two cutting faces intersect on a line L, two tasks must be accomplished. The first is to insert a new node into the intersection-line list for Swait , with each node storing two endpoints for the intersection between L and Swait . The other is to determine whether each of these endpoints would be cut from Swait by other cutting faces. If so, the endpoint should be removed from the node storing it. After the intersection-line list for Swait has been constructed, the last task for determining whether Await is a surface atom is to check the current status of the intersection-line list. If the list is empty, as shown in Fig. 5(a), Await is a surface atom. If the list is not empty, as shown in Fig. 5(b), Await will still be a surface atom if at least one endpoint survives the previous cutting-off operations. In Fig. 5(c) and (d) ((d) is a special case of (c)), all of the endpoints are cut off from the surface of Swait ; that is, no endpoint is left in any nodes in the list. Therefore, the surface of Swait is completely covered by the extended spheres of Await 's neighbors. According to Theorem 2, Await is an internal atom. In summary, this surface atom extraction algorithm based on extended spheres can be implemented as follows: Step 1. Read the basic information for all atoms in a molecule from the PDB file, including their names, center coordinates and radii.
Fig. 3. Schematic diagram for “cutting” an extended sphere. (a) A cutting face formed by two extended spheres and (b) The shared area represents the abandoned part.
Fig. 4. The relationship of two cutting faces. (a) Two cutting faces intersect, (b) Only one cutting face is abandoned, (c) Both of the cutting faces remain and (d) Both of the cutting faces are abandoned.
Fig. 5. The four cases of Swait after being “cut”. (a) No intersection line, (b) At east one endopoint remains, (c) No surface of Swait remains and (d) Swait is all abandoned.
294
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
Step 2. Construct the extended sphere for each atom. Steps 3–6 are applied to each extended sphere to determine whether the corresponding atom is a surface atom. Step 3. Build a cutting-face list for an extended sphere. Step 4. Simplify the cutting-face list. If the atom can be determined as being internal in this step, it is marked as internal. The algorithm then returns to Step 3 to process the next extended sphere. Step 5. Build an intersection-line list for the extended sphere. Step 6. Check the existence of each endpoint. If at least one endpoint remains, the current atom is a surface atom; otherwise, it is an internal atom. If the intersection-line list for an atom is empty, it is a surface atom.
3. Calculation and triangulation of the SES 3.1. Calculating the SES based on surface atoms Because internal atoms are irrelevant to constructing the SES, the calculations for the internal atoms in Connolly's method [7] become unnecessary. In this paper, we only construct the tori [7] formed by pairs of neighboring surface atoms. Our numerical experiments show that the SES generated based on surface atoms is the same as that based on all atoms, while the computing time of the former is significantly less than that of the latter.
triangulation mesh for a saddle patch from a real model. See [17] for more details. However, the method used in [17] does not consider the size uniformity of the triangular faces of different saddle patches. For the sake of convenience, we call the arcs ðξ2 ; ξ3 Þ and ðξ1 ; ξ4 Þ traverse arcs in Fig. 6(b). According to [17], the number of subdivisions for the slicing arcs equals that of the traverse arcs. As a result, the triangular faces of some saddle patches may be quite narrow, see Fig. 7(b). To resolve this issue, a concept called the central angle of a saddle patch is presented. In Fig. 8, the dashed arcs represent the trace of a probe touching two atoms, and the solid curves show the saddle patch boundaries. The central angle θ of ðξ1 ; ξ4 Þ is defined as the central angle of the saddle patch. To maintain size uniformity, the number of subdivisions in the direction of the traverse arcs should depend on θ, while the number of subdivisions in the direction of the slicing arcs should be consistent with the two adjacent triangular concave spherical patches. Fig. 7 (c) clearly shows that our method yields better results than the one used in [17]; that is, a long-narrow saddle patch can be subdivided into a more uniform triangular mesh using our method.
3.2. Triangulating each patch of the SES Triangulating the three types of patches that compose the SES will be discussed in this section. In general, a uniform subdivision scheme is a good choice. 3.2.1. Triangulation of saddle patches Saddle patch triangulation has been discussed in detail in [17]. An example is shown in Fig. 6 [17], where the arcs ðξ1 ; ξ2 Þ, ðξ3 ; ξ4 Þ and ðm3 ; m4 Þ are named the slicing arcs [17]. Fig. 6(c) shows the
Fig. 8. Central angle of a saddle patch.
Fig. 6. A saddle patch and the evaluation of points [17]. (a) An example of a saddle patch, (b) The ovaluations of sample points on a saddle patch and (c) An example of a mesh on a saddle patch.
Fig. 7. Optimization of the number of the subdivisions of each saddle patch. (a) Two saddle patches with different central angles, (b) The same number of subdivisions in two directins [17] and (c) Adaptive number of subdivisions based on the central angle of each saddle patch.
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
3.2.2. Triangulation of triangular concave spherical patches and convex spherical patches The triangulation of triangular concave spherical patches and convex spherical patches has been discussed in detail in [17]. Fig. 9 (a) shows an example of a concave patch and Fig. 9(b) is a schematic diagram for Fig. 9(a). Fig. 9(c) shows an example triangulation mesh for a concave patch using a real model. A template triangulation is used for convex spherical patches, as shown in Fig. 10. The methods in [17] perform well for these two patch types. Therefore, we use the same methods to address the two kinds of patches. See [17] for more details.
295
Here, when triangulating the convex spherical patches, which part of the sphere surface should remain (i.e., can be reached by the probe) and which should be abandoned (i.e., cannot be reached by a probe) is not considered. This issue is discussed in Section 3.3. 3.3. Seamless connection of the three patch types Based on the SES structure, we only need consider connecting triangular concave spherical patches to saddle patches and saddle patches to convex patches.
Fig. 9. Concave patches and the midpoint evaluation [17]. (a) An example of a triangular concavo spherical patch, (b) The evaluation of midpoint between points on a great arc of concave patc and An example of triangulation mesh of a concave patch from a real model.
Fig. 10. Template triangulation for a spherical surface [17]. (a) Initial octahedron, (b) Once subdivision and (c) Twice subdivision.
Fig. 11. The overview of the process to connect saddle patches with convex spherical patches. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
296
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
3.3.1. Connection between triangular concave spherical patches and saddle patches The common boundaries between triangular concave spherical patches and saddle patches are the slicing arcs of the saddle patches. These slicing arcs are both uniformly subdivided when triangulating saddle patches and concave patches. Consequently, these two patches can be seamlessly connected as long as the number of subdivisions is the same.
3.3.2. Connection between saddle patches and convex patches The common boundaries between saddle patches and convex patches are the traverse arcs of the saddle patches. The traverse arcs are uniformly subdivided during the saddle patch triangulation, while the convex patch triangulation does not consider the location of these traverse arcs. Therefore, we use the traverse arcs as constraints to modify the triangular faces of the convex spherical patches along the traverse arcs. Fig. 11 demonstrates the process for connecting the saddle patches to the convex spherical patches. Fig. 11(a) shows a spherical surface and saddle patch with a traverse arc denoted by ðξ1 ; ξ2 Þ. The red represents the saddle patch, the green the spherical surface and the blue the torus on which the red saddle patch lies. Part of the spherical surface is clearly covered by the torus, and the dashed blue arrows point to this covered part. The main idea is to modify the spherical surface to suit the corresponding torus and saddle, as shown in Fig. 11(e). In Fig. 11(a), the circle on which ðξ1 ; ξ2 Þ lies is denoted by K. A triangular face on a spherical surface can be classified three ways based on its position relative to the torus. Faces that are not covered by the torus like T1 are kept, faces that are completely covered like T2 are abandoned, and faces that are partly covered like T3, which is shown in Fig. 11(b), should be adjusted. The covered part of faces similar to T3 should be abandoned. Fig. 12 shows such a face in detail. Fig. 12(a) shows ΔABC, a face
Fig. 12. The relationship of K and ΔABC. (a) See the relationship from side view and (b) See the relationship from top view.
like T3, intersecting the plane of the circle K at points I and IAC. - AB In Fig. 12(b), S is the center of K, and the rays S I AB and S I AC intersect K at points M and N, respectively. The covered part of the face, e.g., the quadrilateral IACIABBC in Fig. 12(a), is subsequently abandoned, and the face is replaced with a new triangle, e.g., ΔAMN in Fig. 12, or quadrilateral, as shown in Fig. 11(c). Some of the points computed for circle K, e.g., M and N in Fig. 12 (b), should be moved to fit the vertices of the saddle patch faces lying along the boundaries between the saddle and convex patches. For example, in Fig. 11(c), point P on K should be moved to the nearest vertex P′ in Fig. 11(d) along K. However, point Q does not need to be moved because it is outside the scope of the saddle patch shown in Fig. 11. Thus, the convex patches are shown in Fig. 11(d), after adjusting as detailed above. The next task is to transform the quadrilaterals in the convex patches to triangular faces that seamlessly join with the adjacent saddle patches. This transformation is easily performed as shown in Fig. 11(d) and (e). Fig. 13 shows a mesh involving a convex patch and saddle patch both before and after the adjustment. Every convex patch should be adjusted according to each adjacent torus until all corresponding tori had been iterated. Finally, a seamless triangular mesh is constructed, as shown in Fig. 14. Please note that we do not consider self-intersecting surfaces that typically form in deep grooves. Methods to address such special cases are being developed.
4. Time complexity Sections 2 and 3 show that the methods for both extracting the surface atoms and joining saddle patches to convex patches are
Fig. 14. The tight mesh for 10 atoms.
Fig. 13. Connection between a convex patch and a saddle patch. (a) a convex patch and a saddle patch before adjusting the mesh and (b) A convex patch and a saddle patch after adjusting the mesh.
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
characterized by the locality of the computation; that is, the computation only relates to an atom and the neighboring atoms. In addition, the number of atoms neighboring any atom is known to always be bounded from above by a constant value [9,10]. Therefore, once which atoms neighbor each other has been determined, a basic traversal process taking O(n2) calculations, both of the above methods require O(C n) calculations in the worst case, where C is a constant value and n is the number of atoms. For the saddle patch triangulation, only one more operation than for the method from [17] is required to calculate the central angle for each saddle patch.
5. Experiments Our algorithms were performed on a PC with an Intel Core i3 processor and 2 GB of main memory. We chose three example molecules, Escherichia coli DNA polymerase (P), myoglobin (M) and nitrogen Paracoccus cytochrome C550 (C). 5.1. Experiments for SAEES The performances of PSAGA [30], MSARA [31] and SAEES were first compared. Figs. 15(a) and (b), 16(a) and (b) and 17(a) from [31] illustrate the PSAGA and MSARA results. The yellow balls are surface atoms, and the green balls are internal atoms. The SAEES results are shown in Figs. 15(c), 16(c), and 17(b), where the red balls are surface atoms and the white balls are internal atoms.
297
Fig. 15 shows the results of the three algorithms for polymerase of E. coli DNA. Fig. 16 shows the results of the three algorithms for myoglobin. Fig. 17 shows the results of MSARA and SAEES for nitrogen Paracoccus cytochrome C550. To evaluate the accuracy of these algorithms, we use the SES calculated by Connolly [7] as the standard for whether an atom is a surface atom. Table 1 shows the number of surface atoms determined by each of the three methods and their runtime. Our experiments indicate that PSAGA cannot extract all of the surface atoms, and the extraction accuracy of MSARA is improved relative to PSAGA but is worse than SAEES, which correctly found all of the surface atoms. Furthermore, the runtimes of PSAGA and MSARA increase rapidly upon increasing the number of atoms. In contrast, SAEES had a better time efficiency than the other two methods. 5.2. Experiments for calculating the SES based on the surface atoms We compared the time required to compute the SES based all of the atoms to that based on the surface atoms. The data are shown in Table 2. According to these data, the computing time was shortened using SAEES. Furthermore, for a given molecule, increasing the proportion of internal atoms seems to improve the proportion of the saved time. Note that all of the atoms would be on the surface of a molecule with a loose atomic arrangement; therefore, the time improvement would be zero. However, compared to the time required to calculate the SES, the time required by SAEES could be
Fig. 15. Result figures for extracting surface atoms of E. coli DNA polymerase. (a) PSAGA [30], (b) MSARA [31] and SAEES. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Fig. 16. Result figures for extracting surface atoms of myoglobin. (a) PSAGA [30], (b) MSARA [31] and (c) SAEES. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
298
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
Fig. 17. Result figures for extracting surface atoms of nitrogen Paracoccus cytochrome C550. (a) MSARA [31] and (b) SAEES. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Table 1 The extracted surface atoms number by three methods for three kinds of molecules (time: s). Label of molecule
PSAGA [30] (time)
MSARA [31] (time)
SAEES (time)
Connolly [7]
Sum of atoms
Proportion of surface atoms (%)
P M C
81 (0.51 s) 676 (2.41 s) –
112 (0.55 s) 975 (2.48 s) 1246 (17.6 s)
115 (0.031 s) 932 (0.296 s) 1490 (0.719 s)
115 932 1490
135 1159 2000
85.2 80.4 74.5
Table 2 Compare the computing time of calculating the SES based on all the atoms with that based on surface atoms (time: s). Label of molecule
Calculating the SES based on all the atoms
Calculating the SES based on the surface atoms (total time ¼ time for extractionþtime for calculating the SES)
Saved time
Improving proportion (%)
P M C
0.610 29.641 126.300
0.531¼ 0.031þ 0.500 21.421 ¼0.296 þ 21.125 88.141¼ 0.719þ87.422
0.079 8.22 38.159
12.95 27.73 30.21
Table 3 The time for triangulating three types of patches and the connection (time: s). Label of Triangular concave molecule spherical patches
Saddle Convex patches patches
Connection Total time for triangulation
P M C
0.125 1.109 1.813
4.656 42.266 68.797
0.220 1.938 3.344
0.734 5.969 10.047
5.735 51.282 84.001
ignored. Therefore, we can use this technique for a variety of molecules without worrying about whether all of the atoms are on the surface. 5.3. Experiments for triangulating the SES In our experiments, the number of subdivisions for concave and convex patches, named Nc, was 3. Nc is also the number of times the saddle patches are subdivided in the direction of the slicing arcs. As for the number of subdivisions in the direction of the traverse arcs for each saddle patch, we choose a minimum integer of Nc(θ/120), where θ is the central angle for each saddle patch in degrees. Table 3 lists the triangulation time for the three types of patches in the SES and their connection.
Table 3 shows that most of the computation time is spent connecting the three types of patches. The algorithm yields good results; however, the computing efficiency still requires improvement.
6. Conclusions and future work The algorithm proposed in this paper is suitable for extracting surface atoms as a preprocess. Because determining whether an atom is a surface atom only depends on the neighboring atoms, this algorithm can be implemented in a parallelized way. Moreover, the method for calculating the SES based on the surface atoms can be parallelized by calculating patches of each surface atom in a processor. In addition, the triangulation is a highly parallel process. In future works, we will implement the above algorithms in a parallel way.
Acknowledgments This work was supported by the National Science Foundation of China (Grant no. 61272032). We thank reviewers for their constructive comments.
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
References [1] Lee B, Richards FM. The interpretation of protein structures: estimation of static accessibility. J Mol Biol 1971;55(3):379–400. [2] Richards FM. Areas, volumes, packing, and protein structure. Ann Rev Biophys Bioeng 1977;6:151–76. [3] Connolly ML. Solvent-accessible surfaces of proteins and nucleic acids. Science 1983;221(4612):709–13. [4] Connolly ML. Molecular surfaces: a review. Netw Sci. 1996 〈http://www.netsci. org/Science/Compchem/feature14.html〉. [5] Greer J, Bush BL. Macromolecular shape and surface maps by solvent exclusion. Proc Natl Acad Sci 1978;75(1):303–7. [6] Kim DS, Won CI, Bhak J. A proposal for the revision of molecular boundary typology. J Biomol Struct Dyn 2010;28(2):277–87. [7] Connolly ML. Analytical molecular surface calculation. J Appl Crystallogr 1983;16:548–58. [8] Connolly ML. Molecular surface triangulation. J Appl Crystallogr 1985;18:499–505. [9] Varshney A, Brooks FP. Fast analytical computation of Richards's smooth molecular surface. In: Proceedings of the 4th conference on visualization'93. IEEE Computer Society; 1993. p. 300–7. [10] Varshney A, Brooks Jr FP, Wright WV. Computing smooth molecular surfaces. IEEE Comput Graph Appl 1994;14(5):19–25. [11] Sanner MF, Olson AJ, Spehner JC. Reduced surface: an efficient way to compute molecular surfaces. Biopolymers 1996;38(3):305–20. [12] Kim DS, Seo J, Kim D, Ryu J, Cho CH. Three-dimensional beta shapes. ComputAided Des 2006;38(11):1179–91. [13] Ryu J, Park R, Kim DS. Molecular surfaces on proteins via beta shapes. ComputAided Des 2007;39(12):1042–57. [14] Edelsbrunner H, Mücke EP. Three-dimensional alpha shapes. ACM Trans Graph (TOG) 1994;13(1):43–72. [15] Edelsbrunner H, Kirkpatrick D, Seidel R. On the shape of a set of points in the plane. IEEE Trans Inf Theory 1983;29(4):551–9. [16] Ryu J, Park R, Seo J, Kim C, Lee HC, Kim DS. Real-time triangulation of molecular surfaces. Comput Sci Appl—ICCSA 2007;2007:55–67.
299
[17] Ryu J, Cho Y, Kim DS. Triangulation of molecular surfaces. Comput-Aided Des 2009;41(6):463–78. [18] Kim DS, Cho Y, Sugihara K. Quasi-worlds and quasi-operators on quasitriangulations. Comput-Aided Des 2010;42(10):874–88. [19] Kim DS, Cho Y, Kim JK, Ryu J. QTF: quasi-triangulation file format. ComputAided Des 2012;44(9):835–45. [20] Akkiraju N, Edelsbrunner H. Triangulating the surface of a molecule. Discrete Appl Math 1996;71(1):5–22. [21] Edelsbrunner H. The union of balls and its dual shape. Discrete Comput Geom 1995;13(1):415–40. [22] Edelsbrunner H. Deformable smooth surface design. Discrete Comput Geom 1999;21(1):87–115. [23] Cheng HL, Dey TK, Edelsbrunner H, Sullivan J. Dynamic skin triangulation. Discrete Comput Geom 2001;25(4):525–68. [24] Cheng HL, Shi X. Guaranteed quality triangulation of molecular skin. In: Proceedings of IEEE visualization. 2004; p. 481–8. [25] Bajaj C, Lee HY, Merkert R, Pascucci V. NURBS based B-rep models for macromolecules and their properties. In: Proceedings of the fourth ACM symposium on Solid modeling and applications. ACM; 1997. p. 217–28. [26] Bajaj CL, Pascucci V, Shamir A, Holt RJ, Netravali AN. Dynamic maintenance and visualization of molecular surfaces. Discrete Appl Math 2003;127 (1):23–51. [27] Dias S, Gomes A. GPU-based triangulation of the van der Waals surface. In: Proceedings of the 2010 IEEE international conference on bioinformatics and biomedicine (BIBM): 2010; p. 673–6. [28] Dias SED, Gomes AJP. Graphics processing unit-based triangulations of Blinn molecular surfaces. Concurr Comput: Pract Exp 2011;23(17):2280–91. [29] Kim B, Kim KJ, Seong JK. GPU accelerated molecular surface computing. Appl Math 2012;6(1S):185S–94S. [30] Wang Chi-she. Simulation algorithm of the protein surface. Comput Eng Appl 2007;43(034):92–4. [31] Zhang JQ, Wu Sanzi. Analysis of algorithms on extracting molecular surface atoms. In: Proceedings of the 2001 IEEE 4th international congress on image and signal processing (CISP), vol. 3; 2011. p. 1339–42.
Contents lists available at ScienceDirect
Computers & Graphics journal homepage: www.elsevier.com/locate/cag
Special Section on CAD/Graphics 2013
Triangulation of molecular surfaces based on extracting surface atoms Jingqiao Zhang, Zhe Shi n School of Computer Engineering and Science of Shanghai University, Shanghai 200444, China
art ic l e i nf o
a b s t r a c t
Article history: Received 5 August 2013 Received in revised form 7 November 2013 Accepted 7 November 2013 Available online 22 November 2013
In this paper, we propose a novel surface atom extraction algorithm. By calculating the relative positions of the extended spheres, this algorithm detects all surface atoms reachable by the probe. Our algorithm overcomes the drawbacks of current methods and can be implemented in real time. This algorithm significantly improved the efficiency of calculating the solvent-excluded surface because only the values for the surface atoms are calculated instead of for all atoms. We also construct a triangular mesh model of the molecular solvent-excluded surface in terms of the surface atoms. By subdividing the saddle patches into more uniform triangular meshes, we improved the existing methods for triangulating a molecular surface composed of triangular concave spherical patches, saddle patches and convex patches. In addition, we introduce a new method that seamlessly joins saddle patches to convex patches. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Solvent-excluded surface Extended sphere Surface atom Triangulation
1. Introduction The three-dimensional structure of a molecule is one of the most important factors determining its function. The molecular surface is a primary part of the protein structure and is crucial to protein docking and visualization. Currently, most research has focused on three molecular surface models [1–6]: the van der Waals surfaces, solvent-accessible surfaces (SASs) and solventexcluded surfaces (SESs). Solvent-excluded surfaces are also commonly known as solvent contact surfaces (SCSs), Connolly surfaces or smooth molecular surfaces. We consider the SES in this work. The SES was first proposed by Richards [2] in 1977, and it considers a probe sphere representing a solvent molecule rolling over the atoms in a molecule. SESs, which have become a wellknown definition for the molecular surface, consist of three types of patches according to Connolly [7], namely, convex spherical patches, saddle patches and triangular concave spherical patches, as shown in Fig. 1. Throughout the rest of this article the term molecular surface refers to SES unless otherwise indicated. In [7,8], Connolly showed how to calculate the three types of patches as well as the molecular area and volume and briefly discussed triangulating the molecular surface. Varshney et al. [9,10] presented a parallelizable algorithm for computing the molecular surface, which needs O(k log k) over n processors, where n is the number of atoms in the molecule and k is the average number of neighbors per atom. Note that the extended
n
Corresponding author. Tel.: þ 86 18817350768. E-mail addresses: [email protected] (J. Zhang), [email protected] (Z. Shi).
0097-8493/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cag.2013.11.007
sphere defined by Varshney et al. [9,10] is important to our algorithm for extracting molecular surface atoms. In 1996, Sanner et al. [11] proposed a reduced surface consisting of a polyhedron with faces, edges and vertices that indicate whether the probe is in contact with more than two atoms, two atoms or only one atom. This reduced surface can be computed for a molecule from a binary spatial division tree, which provides an efficient algorithm for computing the molecular surface. Sanner et al. [11] claimed that the binary spatial division tree can be constructed using O(n log n) computations in the worst case and the reduced surface can be computed using O(n log n) computations on average, where n is the number of atoms in the molecule. In 2006, Ryu et al. [12,13] introduced a β-shape, which is a generalization of the well-known α-shape [14,15] that can be used to efficiently calculate molecular surfaces. Later, Ryu et al. [16,17] discussed the uniform triangulation of the three patch types. However, constructing the β-shape is difficult. First, for a molecule with n atoms, the Voronoi diagram of atoms requires O(n3) calculations in the worst case (the average complexity is O(n)). Second, for a given Voronoi diagram, the β-shape can be computed in O(log n þk) calculations in the worst case, where k is the number of edges and faces in the resulting β-shape [13,18,19]. Nevertheless, the Voronoi diagram of the atoms and β-shape are powerful tools for other important applications such as protein docking problems. The reduced surface [11] is equivalent to one β-shape. In this paper, we improve Ryu's methods for triangulating saddle patches [17]. Our improved method generates a uniform mesh even though the saddle patch is narrow. Several other researchers have studied the calculation and visualization of various molecular surface models in other ways [20–29].
292
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
Theorem 1. The extended spheres of two atoms intersect (except when tangent) if and only if the two atoms are neighbors. Proof. The extended spheres for two given atoms, Ai ðc i ; r i Þ and Aj ðc j ; r j Þ, are Si ðci ; r i þ r p Þ and Sj ðc j ; r j þ r p Þ, respectively, and the distance between their center is D. If these two extended spheres intersect (except when tangent), we obtain D o L ¼ r i þ r p þ r j þ r p . According to Definition 2, Ai and Aj are neighboring atoms. In contrast, two neighbor atoms have D o L ¼ r i þ r p þ r j þr p ; namely, their extended spheres intersect (except when tangent).
Fig. 1. The solvent-excluded surface.
Most previous works researched various molecular surface models based on all of the atoms in the molecule. However, if the surface atoms can be effectively extracted, the scale of the problem could decrease. Because the SES is constructed independently of the internal atoms, the calculation efficiency for SES can be improved by only accounting for the surface atoms. In this paper, we define surface atoms as those that can be reached by a probe rolling over the molecule and internal atoms as those that cannot be reached. For this surface atom extraction, Wang et al. [30] proposed a geometric algorithm named the Protein Surface Atom Geometrical Algorithm (PSAGA) in 2007 to determine the molecular surface atoms according to a molecule's 3D geometric characteristics. However, the results extracted by this algorithm are not exact. Moreover, it requires the radii of all the atoms to be the same. In 2011, Zhang et al. [31] introduced a new extraction algorithm called the Molecular Surface Atom Rolling Algorithm (MSARA), which is more accurate than the PSAGA while being effective for atoms with differing radius. However, this algorithm cannot find certain surface atoms, and its efficiency is not ideal. Here, we present a new, real-time extraction algorithm that overcomes the drawbacks of the existing methods. Our algorithm consists of three steps: (1) extract the surface atoms; (2) calculate and triangulate the three molecule surface patch types; and (3) connect these patches seamlessly. This paper is organized as follows. Section 2 introduces the new extraction algorithm. Section 3 provides the calculations for the SES based on the extracted surface atoms and discusses the triangulation and connection of the three SES patches. The time complexity is analyzed in Section 4. We test our algorithms in Section 5. Conclusions and future work are presented in Section 6.
Theorem 2. If the surface of an extended sphere for one atom in a molecule is completely covered by the interior of other atoms' extended spheres, it is an internal atom; otherwise, it is a surface atom. Proof. If the surface of an atom's extended sphere is not completely covered, the probe can reach this atom without colliding with other atoms whenever the probe's center lies on the uncovered portion of the atom's extended sphere surface. Thus, the atom is a surface atom. In contrast, if the probe cannot reach this atom without colliding with other atoms, it is an internal atom. Theorem 2 provides a method for determining whether an atom is a surface atom. Fig. 2(a) and (b) shows a molecule and the extended spheres of its atoms, respectively. Because the extended sphere surface illustrated by the red dashed circle in Fig. 2(c) is completely covered, the corresponding atom in Fig. 2(d) is internal. On the other hand, because the extended sphere surfaces, which are illustrated by the black solid circles in Fig. 2(c), are not completely covered, their corresponding atoms in Fig. 2(d) are surface atoms. Moreover, according to Theorem 1, whether an atom is a surface atom only depends on it and its neighbors instead of all atoms. We propose the following method to simplify determining whether the surface of an extended sphere is completely covered.
2. Surface Atom Extraction based on Extended Spheres (SAEES) In the rest of this article, we denote r p as the radius of the probe. Definition 1. Let Aðc ; rÞ be an atom, where c is the center of the atom and r is the radius of the atom. The extended sphere of an atom is a sphere Sðc ; r þ r p Þ with the center at c and of the radius ðr þ r p Þ. The definition of an extended sphere indicates that the center of the probe cannot be inside of an atom's extended sphere; otherwise a collision would occur between the probe and the atom. Definition 2. Given two atoms, Ai ðci ; r i Þ and Aj ðc j ; r j Þ, let D be the Euclidean distance between c i and c j and L ¼ r i þ r j þ 2r p . Ai and Aj are neighbors to each other if D o L.
Fig. 2. Determine whether an atom is a surface atom based on the extended spheres. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.). (a) A molecule, (b) The molecule and the extended spheres of all the atoms, (c) All the extended spheres and (d) Surface atoms and an internal atom.
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
Let Await be an atom with an extended sphere of Swait . A1 with the extended sphere S1 neighbors Await . According to Theorem 1, Swait and S1 intersect. As shown in Fig. 3(a), the intersection between Swait 's surface and S1 's surface is a circle named the intersection circle. The planar point set surrounded by the intersection circle is called the cutting face and is denoted with C. P wait and P 1 are the centers of Swait and S1 , respectively. The vector v is the unit vector of P wait P 1 and can be used as the unit normal vector of C. Imagine Swait is a watermelon with a thin peel; the cutting face is used to cut off part of Swait and v points towards the removed parts, as shown in Fig. 3(b). Suppose Await has m neighboring atoms. There are m cutting faces, denoted by C1, C2… Cm, formed by Swait and m extended spheres from Await 's neighboring atoms. These m cutting faces are stored in a cutting-face list, with each node storing the cutting face's center and a unit normal vector pointing to the removed portion of Swait . Thus, Swait is a watermelon cut by Ci(i¼ 1,…,m). After this “cutting”, if some peel of the Swait remains, that is, the surface of Swait is not completely covered, Await is a surface atom. This method works even when the atomic radii are different and can identify the surface atoms at the boundary of internal voids. Before “cutting” Swait , a preprocess is implemented to reduce the number of cutting faces or in some cases directly determine whether Await is a surface atom. During preprocessing, the intersection line L for the plane of any two cutting faces, C1 and C2, (if they are not parallel) is calculated. If the radius of Swait is longer than the distance from the center of Swait to L, C1 and C2 intersect. Fig. 4(a) shows that two cutting faces intersect when both of them remain.
293
If C1 and C2 do not intersect, as shown in Fig. 4(b)–(d), whether one cutting face (e.g., C1) lies in the abandoned part of Swait due to another cutting face (e.g., C2) should be considered and is easily judged by computing the scalar product of the unit normal vector for C2 and the vector from the center of C2 to the center of C1. There are three cases for the relative positions of C1 and C2: (1) when C1 (C2) lies in the abandoned part of Swait due to C2 (C1) while C2 (C1) does not lie in the abandoned part of Swait due to C1 (C2), C1 (C2) can be deleted, see Fig. 4(b); (2) when neither C1 nor C2 lie in the abandoned part of Swait due to the other, then both remain, see Fig. 4(c); (3) when both lie in the abandoned part of Swait due to the other, Await can be directly determined as being an internal atom, see Fig. 4(d). After reducing the number of cutting faces via this preprocessing, we must traverse each remaining pair of cutting faces. Whenever two cutting faces intersect on a line L, two tasks must be accomplished. The first is to insert a new node into the intersection-line list for Swait , with each node storing two endpoints for the intersection between L and Swait . The other is to determine whether each of these endpoints would be cut from Swait by other cutting faces. If so, the endpoint should be removed from the node storing it. After the intersection-line list for Swait has been constructed, the last task for determining whether Await is a surface atom is to check the current status of the intersection-line list. If the list is empty, as shown in Fig. 5(a), Await is a surface atom. If the list is not empty, as shown in Fig. 5(b), Await will still be a surface atom if at least one endpoint survives the previous cutting-off operations. In Fig. 5(c) and (d) ((d) is a special case of (c)), all of the endpoints are cut off from the surface of Swait ; that is, no endpoint is left in any nodes in the list. Therefore, the surface of Swait is completely covered by the extended spheres of Await 's neighbors. According to Theorem 2, Await is an internal atom. In summary, this surface atom extraction algorithm based on extended spheres can be implemented as follows: Step 1. Read the basic information for all atoms in a molecule from the PDB file, including their names, center coordinates and radii.
Fig. 3. Schematic diagram for “cutting” an extended sphere. (a) A cutting face formed by two extended spheres and (b) The shared area represents the abandoned part.
Fig. 4. The relationship of two cutting faces. (a) Two cutting faces intersect, (b) Only one cutting face is abandoned, (c) Both of the cutting faces remain and (d) Both of the cutting faces are abandoned.
Fig. 5. The four cases of Swait after being “cut”. (a) No intersection line, (b) At east one endopoint remains, (c) No surface of Swait remains and (d) Swait is all abandoned.
294
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
Step 2. Construct the extended sphere for each atom. Steps 3–6 are applied to each extended sphere to determine whether the corresponding atom is a surface atom. Step 3. Build a cutting-face list for an extended sphere. Step 4. Simplify the cutting-face list. If the atom can be determined as being internal in this step, it is marked as internal. The algorithm then returns to Step 3 to process the next extended sphere. Step 5. Build an intersection-line list for the extended sphere. Step 6. Check the existence of each endpoint. If at least one endpoint remains, the current atom is a surface atom; otherwise, it is an internal atom. If the intersection-line list for an atom is empty, it is a surface atom.
3. Calculation and triangulation of the SES 3.1. Calculating the SES based on surface atoms Because internal atoms are irrelevant to constructing the SES, the calculations for the internal atoms in Connolly's method [7] become unnecessary. In this paper, we only construct the tori [7] formed by pairs of neighboring surface atoms. Our numerical experiments show that the SES generated based on surface atoms is the same as that based on all atoms, while the computing time of the former is significantly less than that of the latter.
triangulation mesh for a saddle patch from a real model. See [17] for more details. However, the method used in [17] does not consider the size uniformity of the triangular faces of different saddle patches. For the sake of convenience, we call the arcs ðξ2 ; ξ3 Þ and ðξ1 ; ξ4 Þ traverse arcs in Fig. 6(b). According to [17], the number of subdivisions for the slicing arcs equals that of the traverse arcs. As a result, the triangular faces of some saddle patches may be quite narrow, see Fig. 7(b). To resolve this issue, a concept called the central angle of a saddle patch is presented. In Fig. 8, the dashed arcs represent the trace of a probe touching two atoms, and the solid curves show the saddle patch boundaries. The central angle θ of ðξ1 ; ξ4 Þ is defined as the central angle of the saddle patch. To maintain size uniformity, the number of subdivisions in the direction of the traverse arcs should depend on θ, while the number of subdivisions in the direction of the slicing arcs should be consistent with the two adjacent triangular concave spherical patches. Fig. 7 (c) clearly shows that our method yields better results than the one used in [17]; that is, a long-narrow saddle patch can be subdivided into a more uniform triangular mesh using our method.
3.2. Triangulating each patch of the SES Triangulating the three types of patches that compose the SES will be discussed in this section. In general, a uniform subdivision scheme is a good choice. 3.2.1. Triangulation of saddle patches Saddle patch triangulation has been discussed in detail in [17]. An example is shown in Fig. 6 [17], where the arcs ðξ1 ; ξ2 Þ, ðξ3 ; ξ4 Þ and ðm3 ; m4 Þ are named the slicing arcs [17]. Fig. 6(c) shows the
Fig. 8. Central angle of a saddle patch.
Fig. 6. A saddle patch and the evaluation of points [17]. (a) An example of a saddle patch, (b) The ovaluations of sample points on a saddle patch and (c) An example of a mesh on a saddle patch.
Fig. 7. Optimization of the number of the subdivisions of each saddle patch. (a) Two saddle patches with different central angles, (b) The same number of subdivisions in two directins [17] and (c) Adaptive number of subdivisions based on the central angle of each saddle patch.
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
3.2.2. Triangulation of triangular concave spherical patches and convex spherical patches The triangulation of triangular concave spherical patches and convex spherical patches has been discussed in detail in [17]. Fig. 9 (a) shows an example of a concave patch and Fig. 9(b) is a schematic diagram for Fig. 9(a). Fig. 9(c) shows an example triangulation mesh for a concave patch using a real model. A template triangulation is used for convex spherical patches, as shown in Fig. 10. The methods in [17] perform well for these two patch types. Therefore, we use the same methods to address the two kinds of patches. See [17] for more details.
295
Here, when triangulating the convex spherical patches, which part of the sphere surface should remain (i.e., can be reached by the probe) and which should be abandoned (i.e., cannot be reached by a probe) is not considered. This issue is discussed in Section 3.3. 3.3. Seamless connection of the three patch types Based on the SES structure, we only need consider connecting triangular concave spherical patches to saddle patches and saddle patches to convex patches.
Fig. 9. Concave patches and the midpoint evaluation [17]. (a) An example of a triangular concavo spherical patch, (b) The evaluation of midpoint between points on a great arc of concave patc and An example of triangulation mesh of a concave patch from a real model.
Fig. 10. Template triangulation for a spherical surface [17]. (a) Initial octahedron, (b) Once subdivision and (c) Twice subdivision.
Fig. 11. The overview of the process to connect saddle patches with convex spherical patches. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
296
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
3.3.1. Connection between triangular concave spherical patches and saddle patches The common boundaries between triangular concave spherical patches and saddle patches are the slicing arcs of the saddle patches. These slicing arcs are both uniformly subdivided when triangulating saddle patches and concave patches. Consequently, these two patches can be seamlessly connected as long as the number of subdivisions is the same.
3.3.2. Connection between saddle patches and convex patches The common boundaries between saddle patches and convex patches are the traverse arcs of the saddle patches. The traverse arcs are uniformly subdivided during the saddle patch triangulation, while the convex patch triangulation does not consider the location of these traverse arcs. Therefore, we use the traverse arcs as constraints to modify the triangular faces of the convex spherical patches along the traverse arcs. Fig. 11 demonstrates the process for connecting the saddle patches to the convex spherical patches. Fig. 11(a) shows a spherical surface and saddle patch with a traverse arc denoted by ðξ1 ; ξ2 Þ. The red represents the saddle patch, the green the spherical surface and the blue the torus on which the red saddle patch lies. Part of the spherical surface is clearly covered by the torus, and the dashed blue arrows point to this covered part. The main idea is to modify the spherical surface to suit the corresponding torus and saddle, as shown in Fig. 11(e). In Fig. 11(a), the circle on which ðξ1 ; ξ2 Þ lies is denoted by K. A triangular face on a spherical surface can be classified three ways based on its position relative to the torus. Faces that are not covered by the torus like T1 are kept, faces that are completely covered like T2 are abandoned, and faces that are partly covered like T3, which is shown in Fig. 11(b), should be adjusted. The covered part of faces similar to T3 should be abandoned. Fig. 12 shows such a face in detail. Fig. 12(a) shows ΔABC, a face
Fig. 12. The relationship of K and ΔABC. (a) See the relationship from side view and (b) See the relationship from top view.
like T3, intersecting the plane of the circle K at points I and IAC. - AB In Fig. 12(b), S is the center of K, and the rays S I AB and S I AC intersect K at points M and N, respectively. The covered part of the face, e.g., the quadrilateral IACIABBC in Fig. 12(a), is subsequently abandoned, and the face is replaced with a new triangle, e.g., ΔAMN in Fig. 12, or quadrilateral, as shown in Fig. 11(c). Some of the points computed for circle K, e.g., M and N in Fig. 12 (b), should be moved to fit the vertices of the saddle patch faces lying along the boundaries between the saddle and convex patches. For example, in Fig. 11(c), point P on K should be moved to the nearest vertex P′ in Fig. 11(d) along K. However, point Q does not need to be moved because it is outside the scope of the saddle patch shown in Fig. 11. Thus, the convex patches are shown in Fig. 11(d), after adjusting as detailed above. The next task is to transform the quadrilaterals in the convex patches to triangular faces that seamlessly join with the adjacent saddle patches. This transformation is easily performed as shown in Fig. 11(d) and (e). Fig. 13 shows a mesh involving a convex patch and saddle patch both before and after the adjustment. Every convex patch should be adjusted according to each adjacent torus until all corresponding tori had been iterated. Finally, a seamless triangular mesh is constructed, as shown in Fig. 14. Please note that we do not consider self-intersecting surfaces that typically form in deep grooves. Methods to address such special cases are being developed.
4. Time complexity Sections 2 and 3 show that the methods for both extracting the surface atoms and joining saddle patches to convex patches are
Fig. 14. The tight mesh for 10 atoms.
Fig. 13. Connection between a convex patch and a saddle patch. (a) a convex patch and a saddle patch before adjusting the mesh and (b) A convex patch and a saddle patch after adjusting the mesh.
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
characterized by the locality of the computation; that is, the computation only relates to an atom and the neighboring atoms. In addition, the number of atoms neighboring any atom is known to always be bounded from above by a constant value [9,10]. Therefore, once which atoms neighbor each other has been determined, a basic traversal process taking O(n2) calculations, both of the above methods require O(C n) calculations in the worst case, where C is a constant value and n is the number of atoms. For the saddle patch triangulation, only one more operation than for the method from [17] is required to calculate the central angle for each saddle patch.
5. Experiments Our algorithms were performed on a PC with an Intel Core i3 processor and 2 GB of main memory. We chose three example molecules, Escherichia coli DNA polymerase (P), myoglobin (M) and nitrogen Paracoccus cytochrome C550 (C). 5.1. Experiments for SAEES The performances of PSAGA [30], MSARA [31] and SAEES were first compared. Figs. 15(a) and (b), 16(a) and (b) and 17(a) from [31] illustrate the PSAGA and MSARA results. The yellow balls are surface atoms, and the green balls are internal atoms. The SAEES results are shown in Figs. 15(c), 16(c), and 17(b), where the red balls are surface atoms and the white balls are internal atoms.
297
Fig. 15 shows the results of the three algorithms for polymerase of E. coli DNA. Fig. 16 shows the results of the three algorithms for myoglobin. Fig. 17 shows the results of MSARA and SAEES for nitrogen Paracoccus cytochrome C550. To evaluate the accuracy of these algorithms, we use the SES calculated by Connolly [7] as the standard for whether an atom is a surface atom. Table 1 shows the number of surface atoms determined by each of the three methods and their runtime. Our experiments indicate that PSAGA cannot extract all of the surface atoms, and the extraction accuracy of MSARA is improved relative to PSAGA but is worse than SAEES, which correctly found all of the surface atoms. Furthermore, the runtimes of PSAGA and MSARA increase rapidly upon increasing the number of atoms. In contrast, SAEES had a better time efficiency than the other two methods. 5.2. Experiments for calculating the SES based on the surface atoms We compared the time required to compute the SES based all of the atoms to that based on the surface atoms. The data are shown in Table 2. According to these data, the computing time was shortened using SAEES. Furthermore, for a given molecule, increasing the proportion of internal atoms seems to improve the proportion of the saved time. Note that all of the atoms would be on the surface of a molecule with a loose atomic arrangement; therefore, the time improvement would be zero. However, compared to the time required to calculate the SES, the time required by SAEES could be
Fig. 15. Result figures for extracting surface atoms of E. coli DNA polymerase. (a) PSAGA [30], (b) MSARA [31] and SAEES. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Fig. 16. Result figures for extracting surface atoms of myoglobin. (a) PSAGA [30], (b) MSARA [31] and (c) SAEES. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
298
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
Fig. 17. Result figures for extracting surface atoms of nitrogen Paracoccus cytochrome C550. (a) MSARA [31] and (b) SAEES. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
Table 1 The extracted surface atoms number by three methods for three kinds of molecules (time: s). Label of molecule
PSAGA [30] (time)
MSARA [31] (time)
SAEES (time)
Connolly [7]
Sum of atoms
Proportion of surface atoms (%)
P M C
81 (0.51 s) 676 (2.41 s) –
112 (0.55 s) 975 (2.48 s) 1246 (17.6 s)
115 (0.031 s) 932 (0.296 s) 1490 (0.719 s)
115 932 1490
135 1159 2000
85.2 80.4 74.5
Table 2 Compare the computing time of calculating the SES based on all the atoms with that based on surface atoms (time: s). Label of molecule
Calculating the SES based on all the atoms
Calculating the SES based on the surface atoms (total time ¼ time for extractionþtime for calculating the SES)
Saved time
Improving proportion (%)
P M C
0.610 29.641 126.300
0.531¼ 0.031þ 0.500 21.421 ¼0.296 þ 21.125 88.141¼ 0.719þ87.422
0.079 8.22 38.159
12.95 27.73 30.21
Table 3 The time for triangulating three types of patches and the connection (time: s). Label of Triangular concave molecule spherical patches
Saddle Convex patches patches
Connection Total time for triangulation
P M C
0.125 1.109 1.813
4.656 42.266 68.797
0.220 1.938 3.344
0.734 5.969 10.047
5.735 51.282 84.001
ignored. Therefore, we can use this technique for a variety of molecules without worrying about whether all of the atoms are on the surface. 5.3. Experiments for triangulating the SES In our experiments, the number of subdivisions for concave and convex patches, named Nc, was 3. Nc is also the number of times the saddle patches are subdivided in the direction of the slicing arcs. As for the number of subdivisions in the direction of the traverse arcs for each saddle patch, we choose a minimum integer of Nc(θ/120), where θ is the central angle for each saddle patch in degrees. Table 3 lists the triangulation time for the three types of patches in the SES and their connection.
Table 3 shows that most of the computation time is spent connecting the three types of patches. The algorithm yields good results; however, the computing efficiency still requires improvement.
6. Conclusions and future work The algorithm proposed in this paper is suitable for extracting surface atoms as a preprocess. Because determining whether an atom is a surface atom only depends on the neighboring atoms, this algorithm can be implemented in a parallelized way. Moreover, the method for calculating the SES based on the surface atoms can be parallelized by calculating patches of each surface atom in a processor. In addition, the triangulation is a highly parallel process. In future works, we will implement the above algorithms in a parallel way.
Acknowledgments This work was supported by the National Science Foundation of China (Grant no. 61272032). We thank reviewers for their constructive comments.
J. Zhang, Z. Shi / Computers & Graphics 38 (2014) 291–299
References [1] Lee B, Richards FM. The interpretation of protein structures: estimation of static accessibility. J Mol Biol 1971;55(3):379–400. [2] Richards FM. Areas, volumes, packing, and protein structure. Ann Rev Biophys Bioeng 1977;6:151–76. [3] Connolly ML. Solvent-accessible surfaces of proteins and nucleic acids. Science 1983;221(4612):709–13. [4] Connolly ML. Molecular surfaces: a review. Netw Sci. 1996 〈http://www.netsci. org/Science/Compchem/feature14.html〉. [5] Greer J, Bush BL. Macromolecular shape and surface maps by solvent exclusion. Proc Natl Acad Sci 1978;75(1):303–7. [6] Kim DS, Won CI, Bhak J. A proposal for the revision of molecular boundary typology. J Biomol Struct Dyn 2010;28(2):277–87. [7] Connolly ML. Analytical molecular surface calculation. J Appl Crystallogr 1983;16:548–58. [8] Connolly ML. Molecular surface triangulation. J Appl Crystallogr 1985;18:499–505. [9] Varshney A, Brooks FP. Fast analytical computation of Richards's smooth molecular surface. In: Proceedings of the 4th conference on visualization'93. IEEE Computer Society; 1993. p. 300–7. [10] Varshney A, Brooks Jr FP, Wright WV. Computing smooth molecular surfaces. IEEE Comput Graph Appl 1994;14(5):19–25. [11] Sanner MF, Olson AJ, Spehner JC. Reduced surface: an efficient way to compute molecular surfaces. Biopolymers 1996;38(3):305–20. [12] Kim DS, Seo J, Kim D, Ryu J, Cho CH. Three-dimensional beta shapes. ComputAided Des 2006;38(11):1179–91. [13] Ryu J, Park R, Kim DS. Molecular surfaces on proteins via beta shapes. ComputAided Des 2007;39(12):1042–57. [14] Edelsbrunner H, Mücke EP. Three-dimensional alpha shapes. ACM Trans Graph (TOG) 1994;13(1):43–72. [15] Edelsbrunner H, Kirkpatrick D, Seidel R. On the shape of a set of points in the plane. IEEE Trans Inf Theory 1983;29(4):551–9. [16] Ryu J, Park R, Seo J, Kim C, Lee HC, Kim DS. Real-time triangulation of molecular surfaces. Comput Sci Appl—ICCSA 2007;2007:55–67.
299
[17] Ryu J, Cho Y, Kim DS. Triangulation of molecular surfaces. Comput-Aided Des 2009;41(6):463–78. [18] Kim DS, Cho Y, Sugihara K. Quasi-worlds and quasi-operators on quasitriangulations. Comput-Aided Des 2010;42(10):874–88. [19] Kim DS, Cho Y, Kim JK, Ryu J. QTF: quasi-triangulation file format. ComputAided Des 2012;44(9):835–45. [20] Akkiraju N, Edelsbrunner H. Triangulating the surface of a molecule. Discrete Appl Math 1996;71(1):5–22. [21] Edelsbrunner H. The union of balls and its dual shape. Discrete Comput Geom 1995;13(1):415–40. [22] Edelsbrunner H. Deformable smooth surface design. Discrete Comput Geom 1999;21(1):87–115. [23] Cheng HL, Dey TK, Edelsbrunner H, Sullivan J. Dynamic skin triangulation. Discrete Comput Geom 2001;25(4):525–68. [24] Cheng HL, Shi X. Guaranteed quality triangulation of molecular skin. In: Proceedings of IEEE visualization. 2004; p. 481–8. [25] Bajaj C, Lee HY, Merkert R, Pascucci V. NURBS based B-rep models for macromolecules and their properties. In: Proceedings of the fourth ACM symposium on Solid modeling and applications. ACM; 1997. p. 217–28. [26] Bajaj CL, Pascucci V, Shamir A, Holt RJ, Netravali AN. Dynamic maintenance and visualization of molecular surfaces. Discrete Appl Math 2003;127 (1):23–51. [27] Dias S, Gomes A. GPU-based triangulation of the van der Waals surface. In: Proceedings of the 2010 IEEE international conference on bioinformatics and biomedicine (BIBM): 2010; p. 673–6. [28] Dias SED, Gomes AJP. Graphics processing unit-based triangulations of Blinn molecular surfaces. Concurr Comput: Pract Exp 2011;23(17):2280–91. [29] Kim B, Kim KJ, Seong JK. GPU accelerated molecular surface computing. Appl Math 2012;6(1S):185S–94S. [30] Wang Chi-she. Simulation algorithm of the protein surface. Comput Eng Appl 2007;43(034):92–4. [31] Zhang JQ, Wu Sanzi. Analysis of algorithms on extracting molecular surface atoms. In: Proceedings of the 2001 IEEE 4th international congress on image and signal processing (CISP), vol. 3; 2011. p. 1339–42.