# Molecular Surface Mesh Generation by Filtering Electron Density Map

Authors: Joachim Giard; Benoît Macq
Journal: International Journal of Biomedical Imaging
 (2010)

---

## Abstract

Bioinformatics applied to macromolecules are now widely spread and in continuous expansion. In this context, representing external molecular surface such as the Van der Waals Surface or the Solvent Excluded Surface can be useful for several applications. We propose a fast and parameterizable algorithm giving good visual quality meshes representing molecular surfaces. It is obtained by isosurfacing a filtered electron density map. The density map is the result of the maximum of Gaussian functions placed around atom centers. This map is filtered by an ideal low-pass filter applied on the Fourier Transform of the density map. Applying the marching cubes algorithm on the inverse transform provides a mesh representation of the molecular surface.


---

## Body

## 1. Introduction

The geometric structure of macromolecules, such as proteins or nucleic acids, is directly related to their function [1–3]. Consequently, studying this structure is of capital importance in the understanding and simulation of numerous life processes. It allows researchers to save a lot of time and money for various applications such as drug design [4, 5] or mutation effect prediction [6, 7]. In this context, working with molecules external surface can be useful, for instance, to predict the geometrical complementarity between two molecules [8] or to visualize them [9]. The prediction of geometric complementarity is one of the keystones of molecular docking [10–13], the modeling of interactions between molecules. The localization of potential binding sites of molecules [14–16] is a frequently used tool for docking and generally requires a good description of the protein surface [17–19]. The estimation of the surface area may also be related to the stability of a particular molecule 3D conformation [20].At first, the external surface of a molecule has to be defined. Indeed, molecules are made of atoms which have no real surface. The most frequent molecular surface representations are the Van der Waals Surface (VdWS), the Solvent Accessible Surface (SAS), and the Solvent Excluded Surface (SES) [21, 22]. In the case of the VdWS, the electron clouds around atoms are approximated by rigid spheres with radii corresponding to the Van der Waals (VdW) radii of the atoms. The SAS (resp., SES) is the inner surface of the volume filled by the possible positions of the center (resp., exterior surface) of a ball representing a molecule of solvent, for example, water (see Figure 1).Figure 1
A cutaway view of a small molecule. The most frequent molecular surface representations are detailed. Gray discs depict the Van der Waals volumes, the gray outer line depicts the Solvent Accessible Surface (SAS) and the continuous black line depicts the Solvent Excluded Surface (SES). The SAS (resp., SES) is the limit surface for the solvent molecule center (resp., external surface). The solvent molecule is represented as a black dashed circle.Efficient tools to represent such surfaces are the polygonal meshes, which are collection of points related by edges and faces that approximate the considered surfaces. A lot of methods have been proposed in the last few years for the generation of a molecular surface meshes.However, the computational time remains generally high for quality meshes, and it can be a problem when there is a great amount of data to treat. In this paper, we introduce the Filtered Density Map (FDM) algorithm, which is a fast and parameterizable algorithm to generate smooth molecular surface meshes. The generated mesh is the isosurface of frequency filtered electron density map.This paper is organized as follows. First, some other works related to molecular surface generation are succinctly described in Section2. Then, the FDM method is described in details in Section 3. And finally, results and comparisons with other methods are presented and discussed in Section 4.


## 2. Related Work

In the last few years, a lot of methods have been developed for the generation of molecular surface meshes. In 1983, Connolly [23] proposed an analytical algorithm in which points were strategically placed around the molecule with a specific analytical role (maximum, minimum or saddle point) depending on the number of atoms present in the neighborhood. In 2003, Bajaj et al. [24] introduced another analytical method based on NURBS that offers the advantage to be parameterizable without recalculation. In 2002, Laug and Borouchaki [25] used a parametric representation of intersecting spheres to create the surface mesh. MSMS, developed by Sanner et al. [26] is based on alpha-shapes [27] of molecules. This algorithm is widely used because it is time efficient. However, the generated mesh is not a manifold and is composed of very irregular triangles. The beta-shapes [28] are a generalization of the alpha-shapes and were used by Ryu et al. [29] in 2007 to design a similar algorithm. Another vertex based method was used by Cheng and Shi [30]. In this method, molecular surfaces are generated with the help of restricted union of balls. Finally, some methods based on volumetric computation exist, such as the one of Zhang et al. [31] in which the solvent accessible surface is seen as the isosurface of Gaussian shaped electron density maps, and the algorithm of Can et al. [32] (the LSMS) which is based on a front propagation from atom center and on level-sets.Comparisons between the FDM method and the methods mentioned in this section are shown in Section4.2.


## 3. Method

The FDM method is based on volumetric electron density and a frequency filtering. Each atom is seen as a Gaussian electron cloud, the dimensions of which are depending on the VdW radius. Then, the electron density map is created by taking the local maximum value of these clouds. After a Fourier Transform, it is filtered by an ideal low pass filter, in order to remove frequencies corresponding to a spatial element smaller than a solvent molecule. Finally, a marching cubes [33] algorithm is used on the inverse Fourier Transform to find an isosurface. A refinement of the final mesh constitutes an optional step of the method. The whole algorithm was implemented in C++ with vtk (Visual ToolKit) (http://www.vtk.org/).


### 3.1. Electron Density Map

A Gaussian function is constructed around each atom. The value of this function at a pointi for an atom a is:(1)Ga(i)=te-(ra2-∥ai∥2/r2),
where t is a threshold parameter, r is a radius parameter, ∥ai∥ is the Euclidean distance between the center of a and the point i, ra is the VdW radius of the atom a. So, the isosurface for the threshold t is the VdW sphere because if i is located on the VdWS of a, ∥ai∥=ra and Ga(i)=t. In this work, r is set to 3 Å  because this value is suitable for an ideal low-pass filter (see Section 3.2).For the implementation, the three-dimensional space is divided into voxels. The spacing (Te) between voxels is an adaptable parameter. The more Te is small, the more the surface approximation is fine.The density map of the whole molecule for a point in the space is defined as the maximal value of all the Gaussian functions at this point. The maximum of the Gaussian functions is chosen instead of the summation because it is not possible to evaluate the SES using the isosurface of a summation of Gaussian functions. It can be shown by the following counterexample, in which the Gaussian affected to the atomsa and b must have contradictory properties depending on the situation.In the first situation, the space betweena and b is just small enough to block the way to a solvent molecule (see Figure 2(a)). The other atoms are considered to be too far to have an influence. Thus, the SES has a concave shape at this place and the value of the density map at the “center” of the concavity c1 must be influenced by the fields of a and b: Ga(c1)>0, Gb(c1)>0. We can state that ∥ac1∥>ra+rs because s, the center of the solvent molecule, can be very close to the ab axis.Two situations in which the Gaussian functions would have contradictory properties. (a) Two atoms of the molecule,a and b, are just close enough to block the way to a solvent molecule s. The point c1 is on the SES and must be influenced by the filed generated by a, Ga. (b) a and b are close. The point c2 is on the SES but also on the VdWS. So, it should not be influenced by Ga. It is possible to show that in some (frequent) configurations, ∥ac1∥>∥ac2∥.
(a)(b)In the second situation,∥ab∥<ra+rs/2, which is often the case for covalent bonds. Let c2 be a point belonging both to the SES and to the VdWS of b with the necessary condition ∥bc2∥=rb (see Figure 2(b)). The point c2 should not be influenced by the field of a, so Ga(c2)=0. If ∥ac2∥>∥ac1∥, this condition is in contradiction with Ga(c1)>0, because the Gaussian function is strictly decreasing in the positive domain. Using the Al-Kashi theorem, we know that(2)∥ac2∥2=rb2+∥ab∥2-2rb∥ab∥cosβ,∥ac1∥2>(ra+rs)2=(rb+rs)2+∥ab∥2-2(rb+rs)∥ab∥cosβ,
where β is the abf̂ angle. Thus, ∥ac1∥>∥ac2∥ if(3)rs2+2rbrs-2rs∥ab∥cosβ>0,
what is verified by the hypothesis: ∥ab∥<rb+rs/2 because cosβ≤1.In order to avoid interferences, the maximum is preferred to the summation of Gaussian functions. Isosurfacing this density map returns the VdWS. This surface is not smooth and in order to compute the SES, the density map must first be filtered.


### 3.2. Fourier Transform and Filtering

The Fourier Transform of this electron density map is computed using the FFT algorithm [34]. The frequency representation of the function is filtered by an ideal low pass filter in order to eliminate frequencies corresponding to elements smaller than a solvent molecule, for example, inflexion points between two VdW spheres.The cutoff frequency is fc=1/4rs, where rs is the radius of the sphere approximating the solvent molecule (typically 1.4 Å  for water). The wavelength must be four times longer than rs because a molecule solvent diameter has to fit in a half wavelength (see Figure 3).Figure 3
Minimal wavelength allowed on the SES. It corresponds to four times the solvent radiusrs and determines the cutoff frequency for the filtering.Gaussian functions are preferred to balls in the spatial domain because an ideal low-pass filter makes the Gibb's phenomenon appear on sharp edges. An ideal filter is used because the cutoff frequency is exactly known and because it is numerically possible. An ideal low-pass filter in the frequency domain is equivalent to a convolution product with asinc function in the space domain. Let X=(x,y,z) be the space variable in ℝ3 and ℳ(X) the initial electron density map. Then, the filtered density map is:(4)ℳ̃(X)=(ℳ(•)*2fcsinc(2fc•))(X),
where • represents the variable. The parameter r of the Gaussian functions in (1) is related to the width of the function. To keep the isosurface at the same place, a wider function has a smaller maximum. If this maximum is too high, that is, if r is too small, the secondary ripples of the sinc function take too much importance when they are in phase with this maximum. It makes oscillations appear in the final density map, which can lead to the apparition of unwanted surfaces after isosurfacing. Simulations with several 1D and 2D functions were performed to verify the effect of r. The three main conditions to verify are the follwing:(5)ℳ̃(Xc1)=tforc1lyingontheSESandontheVdWS,ℳ̃(Xc2)=tforc2lyingontheSESbutnotontheVdWS,ℳ̃(Xc3)<tforc3lyingoutsidethemolecule.
Here is a 2D example:(6)ℳ(x)=max(Ga(x),Gb(x)),
with the atom a centered in (xa,ya)=(0,-2.5) and the atom b in (xb,yb)=(0,2.5), the threshold t=1, and ra=rb=1.8 (Figure 4). The values of the density maps ℳ(X) and ℳ̃(X) for c1=(0,xa-ra), c2=(xa2+(ra+rs)2-rs,0), and c3=(x3,y3) are plotted as a function of r in Figure 5. (x3,y3) is the position of the maximal value of the density map outside the “molecule” for r=1. When the parameter r=3, conditions (5) are verified and the Gaussian functions are not too wide, what leads to shorter execution times. The 2D density maps ℳ(X) and ℳ̃(X) with r=3 as well as ℳ̃(X) with r=1 are shown in Figure 6. The isocontours, representing the VdWS or the SES, are depicted in white and we can see the artifacts appearing for too small values of r.Figure 4
In this example, the atomsa and b are 5 Å  far from each other and have both a VdW radius of 1.8 Å. The radius of the solvent molecule s is 1.4 Å. The conditions to verify for the density map before (ℳ) and after (ℳ̃) filtering are: ℳ(c1)=ℳ̃(c1)=t, ℳ(c2)<ℳ̃(c2)=t, and ℳ(c3)<tandℳ̃(c3)<t, with t, the threshold value for the isosurfacing.Value of the density map before (continuous gray line) and after (dashed black line) filtering at a point belonging to both the SES and the VdWS,c1 (a) a point belonging to the SES but not the VdWS, c2 (b) and a point outside the molecule, c1 (c) (see Figure 4). t=1, so, when r=3, conditions (5) are satisfied.
(a)(b)(c)Value of the density map for the example of Figure4. x and y are the two spatial dimensions (Å) and the colors represent the value of the density map before filtering and with r=3 (a) after filtering and with r=3 (b) and after filtering with r=1 (c). The contour for t=1 is depicted with a white line. It represents the VdWS (a), the SES (b), and the SES with artifacts (c).
(a)(b)(c)It is important to notice that if the spacing (Te) for the spatial sampling is too large, there would be no filtering. Normally, the sampling frequency, fe=1/Te, has to verify the Nyquist-Shannon theorem: fe>2fmax, where fmax is the higher frequency with a nonzero coefficient in the original signal. However, for this application the errors resulting from a subsampling are not too important and the sampling frequency is chosen such that fe>2fc, that is, Te<2rs. In this situation, the filtering is always possible.


### 3.3. Isosurfacing

The final triangular mesh is an approximation of the isosurface of the filtered electron density map. The most popular technique to extract an isosurface from a 3D image is the marching cubes algorithm [33]. In this algorithm, the voxels are screened by group of eight sharing a same point. Mesh vertices, faces, and edges are added depending on the value of these eight voxels. There are 256 (28) possibilities that can be reduced to 15 situations thanks to symmetries and complementarities.


### 3.4. Refinement

The visual appearance of the final mesh can be improved by magnifying the number of vertices. The number of vertices is increased using a smooth interpolation scheme such as the piecewise smooth surface reconstruction of Hoppe et al. [35], or the algorithm based on the butterfly scheme proposed by Zorin et al. [36].


## 3.1. Electron Density Map

A Gaussian function is constructed around each atom. The value of this function at a pointi for an atom a is:(1)Ga(i)=te-(ra2-∥ai∥2/r2),
where t is a threshold parameter, r is a radius parameter, ∥ai∥ is the Euclidean distance between the center of a and the point i, ra is the VdW radius of the atom a. So, the isosurface for the threshold t is the VdW sphere because if i is located on the VdWS of a, ∥ai∥=ra and Ga(i)=t. In this work, r is set to 3 Å  because this value is suitable for an ideal low-pass filter (see Section 3.2).For the implementation, the three-dimensional space is divided into voxels. The spacing (Te) between voxels is an adaptable parameter. The more Te is small, the more the surface approximation is fine.The density map of the whole molecule for a point in the space is defined as the maximal value of all the Gaussian functions at this point. The maximum of the Gaussian functions is chosen instead of the summation because it is not possible to evaluate the SES using the isosurface of a summation of Gaussian functions. It can be shown by the following counterexample, in which the Gaussian affected to the atomsa and b must have contradictory properties depending on the situation.In the first situation, the space betweena and b is just small enough to block the way to a solvent molecule (see Figure 2(a)). The other atoms are considered to be too far to have an influence. Thus, the SES has a concave shape at this place and the value of the density map at the “center” of the concavity c1 must be influenced by the fields of a and b: Ga(c1)>0, Gb(c1)>0. We can state that ∥ac1∥>ra+rs because s, the center of the solvent molecule, can be very close to the ab axis.Two situations in which the Gaussian functions would have contradictory properties. (a) Two atoms of the molecule,a and b, are just close enough to block the way to a solvent molecule s. The point c1 is on the SES and must be influenced by the filed generated by a, Ga. (b) a and b are close. The point c2 is on the SES but also on the VdWS. So, it should not be influenced by Ga. It is possible to show that in some (frequent) configurations, ∥ac1∥>∥ac2∥.
(a)(b)In the second situation,∥ab∥<ra+rs/2, which is often the case for covalent bonds. Let c2 be a point belonging both to the SES and to the VdWS of b with the necessary condition ∥bc2∥=rb (see Figure 2(b)). The point c2 should not be influenced by the field of a, so Ga(c2)=0. If ∥ac2∥>∥ac1∥, this condition is in contradiction with Ga(c1)>0, because the Gaussian function is strictly decreasing in the positive domain. Using the Al-Kashi theorem, we know that(2)∥ac2∥2=rb2+∥ab∥2-2rb∥ab∥cosβ,∥ac1∥2>(ra+rs)2=(rb+rs)2+∥ab∥2-2(rb+rs)∥ab∥cosβ,
where β is the abf̂ angle. Thus, ∥ac1∥>∥ac2∥ if(3)rs2+2rbrs-2rs∥ab∥cosβ>0,
what is verified by the hypothesis: ∥ab∥<rb+rs/2 because cosβ≤1.In order to avoid interferences, the maximum is preferred to the summation of Gaussian functions. Isosurfacing this density map returns the VdWS. This surface is not smooth and in order to compute the SES, the density map must first be filtered.


## 3.2. Fourier Transform and Filtering

The Fourier Transform of this electron density map is computed using the FFT algorithm [34]. The frequency representation of the function is filtered by an ideal low pass filter in order to eliminate frequencies corresponding to elements smaller than a solvent molecule, for example, inflexion points between two VdW spheres.The cutoff frequency is fc=1/4rs, where rs is the radius of the sphere approximating the solvent molecule (typically 1.4 Å  for water). The wavelength must be four times longer than rs because a molecule solvent diameter has to fit in a half wavelength (see Figure 3).Figure 3
Minimal wavelength allowed on the SES. It corresponds to four times the solvent radiusrs and determines the cutoff frequency for the filtering.Gaussian functions are preferred to balls in the spatial domain because an ideal low-pass filter makes the Gibb's phenomenon appear on sharp edges. An ideal filter is used because the cutoff frequency is exactly known and because it is numerically possible. An ideal low-pass filter in the frequency domain is equivalent to a convolution product with asinc function in the space domain. Let X=(x,y,z) be the space variable in ℝ3 and ℳ(X) the initial electron density map. Then, the filtered density map is:(4)ℳ̃(X)=(ℳ(•)*2fcsinc(2fc•))(X),
where • represents the variable. The parameter r of the Gaussian functions in (1) is related to the width of the function. To keep the isosurface at the same place, a wider function has a smaller maximum. If this maximum is too high, that is, if r is too small, the secondary ripples of the sinc function take too much importance when they are in phase with this maximum. It makes oscillations appear in the final density map, which can lead to the apparition of unwanted surfaces after isosurfacing. Simulations with several 1D and 2D functions were performed to verify the effect of r. The three main conditions to verify are the follwing:(5)ℳ̃(Xc1)=tforc1lyingontheSESandontheVdWS,ℳ̃(Xc2)=tforc2lyingontheSESbutnotontheVdWS,ℳ̃(Xc3)<tforc3lyingoutsidethemolecule.
Here is a 2D example:(6)ℳ(x)=max(Ga(x),Gb(x)),
with the atom a centered in (xa,ya)=(0,-2.5) and the atom b in (xb,yb)=(0,2.5), the threshold t=1, and ra=rb=1.8 (Figure 4). The values of the density maps ℳ(X) and ℳ̃(X) for c1=(0,xa-ra), c2=(xa2+(ra+rs)2-rs,0), and c3=(x3,y3) are plotted as a function of r in Figure 5. (x3,y3) is the position of the maximal value of the density map outside the “molecule” for r=1. When the parameter r=3, conditions (5) are verified and the Gaussian functions are not too wide, what leads to shorter execution times. The 2D density maps ℳ(X) and ℳ̃(X) with r=3 as well as ℳ̃(X) with r=1 are shown in Figure 6. The isocontours, representing the VdWS or the SES, are depicted in white and we can see the artifacts appearing for too small values of r.Figure 4
In this example, the atomsa and b are 5 Å  far from each other and have both a VdW radius of 1.8 Å. The radius of the solvent molecule s is 1.4 Å. The conditions to verify for the density map before (ℳ) and after (ℳ̃) filtering are: ℳ(c1)=ℳ̃(c1)=t, ℳ(c2)<ℳ̃(c2)=t, and ℳ(c3)<tandℳ̃(c3)<t, with t, the threshold value for the isosurfacing.Value of the density map before (continuous gray line) and after (dashed black line) filtering at a point belonging to both the SES and the VdWS,c1 (a) a point belonging to the SES but not the VdWS, c2 (b) and a point outside the molecule, c1 (c) (see Figure 4). t=1, so, when r=3, conditions (5) are satisfied.
(a)(b)(c)Value of the density map for the example of Figure4. x and y are the two spatial dimensions (Å) and the colors represent the value of the density map before filtering and with r=3 (a) after filtering and with r=3 (b) and after filtering with r=1 (c). The contour for t=1 is depicted with a white line. It represents the VdWS (a), the SES (b), and the SES with artifacts (c).
(a)(b)(c)It is important to notice that if the spacing (Te) for the spatial sampling is too large, there would be no filtering. Normally, the sampling frequency, fe=1/Te, has to verify the Nyquist-Shannon theorem: fe>2fmax, where fmax is the higher frequency with a nonzero coefficient in the original signal. However, for this application the errors resulting from a subsampling are not too important and the sampling frequency is chosen such that fe>2fc, that is, Te<2rs. In this situation, the filtering is always possible.


## 3.3. Isosurfacing

The final triangular mesh is an approximation of the isosurface of the filtered electron density map. The most popular technique to extract an isosurface from a 3D image is the marching cubes algorithm [33]. In this algorithm, the voxels are screened by group of eight sharing a same point. Mesh vertices, faces, and edges are added depending on the value of these eight voxels. There are 256 (28) possibilities that can be reduced to 15 situations thanks to symmetries and complementarities.


## 3.4. Refinement

The visual appearance of the final mesh can be improved by magnifying the number of vertices. The number of vertices is increased using a smooth interpolation scheme such as the piecewise smooth surface reconstruction of Hoppe et al. [35], or the algorithm based on the butterfly scheme proposed by Zorin et al. [36].


## 4. Results and Discussion

Some numerical results pointing out advantages and drawbacks of the FDM are shown in this section. The main characteristics to be observed are the computation time and the quality of the generated mesh. The section is divided into three parts: the analysis of the effects of the different parameters of the FDM, the results of computation time comparisons with other existing methods, and a quality measurement of the generated meshes.


### 4.1. Parameters

There are three main parameters modifiable by the user. First, the spatial spacingTe, that is, the distance between two neighbor voxels center, which determines the total number of voxels. Second, the cutoff frequency fc, which determines the smoothness of the final mesh. And third, the refinement rate kr, that is, the number of new points in a triangle for the final mesh magnification. In this section, the effect of these parameters on the visual quality, on the computation time and on the memory space, are discussed.


#### 4.1.1. Spatial Spacing

With a small spatial spacing, it is possible to represent fine details. However, it drastically increases the memory space needed as well as the computation time. Indeed, reducingTe by a factor α increases the number of voxels by α3. The parts of the method depending on the number of voxels are the creation of the density map (time: 𝒪(nv) and space: 𝒪(nv), with nv the number of voxels) and the Fast Fourier Transform (time: 𝒪(nvlognv) and space: 𝒪(nv)). A visual comparison between meshes generated with Te=1.9rs, Te=rs and Te=rs/3 is shown in Figure 7. In these examples, fc=1/4rs and kr=0, that is, the meshes represent the SES without final refinement.Meshes generated from electron density map at a spatial spacing ofTe=1.9rs (a), Te=rs (b), and Te=rs/3 (c) (PDB code: 3EBZ).
(a)(b)(c)


#### 4.1.2. Cutoff Frequency

In order to generate a mesh representing the SES,fc is set to 1/4rs (see Section 3.2). However, depending on the application, the surface could be other than the SES. For instance, if fc>fe/2, there is no actual filtering, and so, the generated mesh represents the VdWS. On the other hand, to obtain a smooth approximation of the molecule shape, fc can be reduced. It is equivalent to consider a bigger solvent molecule. Changing this parameter does not have any effect neither on the computation time nor on the memory space needed. A visual comparison between meshes generated with fc>fe/2 (VdWS), fc=1/4rs, and fc=1/16rs (shape approximation) is shown in Figure 8. In these examples, Te=rs/3 and kr=0, that is, a solvent molecule diameter takes 6 voxels and there is no final refinement.Meshes generated from electron density map filtered at a cutoff frequency offc>fe/2 (a) to give the VdWS, fc=1/4rs (b) to give the SES, and fc=1/16rs (c) to give an approximation of the general shape. (PDB code: 3EBZ).
(a)(b)(c)


#### 4.1.3. Refinement Factor

The final mesh refinement gives foremost an esthetic advantage. The memory space needed does not increase a lot because the number of voxels remains the same. Only the size of the mesh changes and this is negligible in comparison with the space needed by the voxels representation. The computation time slightly increases but, whenkr=1 (which gives a good visual improving), this is negligible in comparison with the voxels operations computation time. A visual comparison between meshes generated with kr=0, kr=1 and kr=2 is shown in Figure 9. In these examples, Te=rs and fc=1/4rs, that is, a solvent molecule diameter takes 2 voxels and the meshes represent the SES.Meshes refined with a factor ofkr=0 (a), kr=1 (b), and kr=2 (c) (PDB code: 3EBZ).
(a)(b)(c)


### 4.2. Time Comparisons

In this section, computation times are compared between the FDM algorithm and algorithms found in the literature for equivalent qualities. When available, the algorithms were run on the same computer, when not, the computation times were the ones announced in the original paper. Can et al. made a comparison of their method computation time with three molecular visualization tools: UCSF Chimera [37], Swiss-PDBViewer [38], and PyMol [9]. We added MSMS [26] to this set of methods. These programs, as well as the one of Can (LSMS), are available for free, so, the computation time could be measured on the same computer than for our method, except for Chimera that was not supported by the system. Thus, the computation times mentioned here for Chimera are the one announced in the paper of Can et al. [32]. For the LSMS method, the grid size is set to 256 × 256 × 256. In this condition, Te≃1 with the tested molecules. Other programs are run with defaults settings, that is, Te≃1. Two tests are made for the FDM method. In the first one, the parameters are set to Te=1.9rs and kr=0 to be as fast as possible while keeping a correct solution. In the second one, the parameters are set to Te=4/3rs and kr=1 which gives a correct mesh with a good appearance (see Figure 10). The computation times are shown in Table 1. The computation times reported only include the mesh generation time, that is, it does not take the loading time into account. In addition, for three molecules, computation times for the method of Cheng and Shi [30] are reported from their paper. The computation times comparison is shown in Table 2. All the tests were performed on a AMD Athlon(tm) 64 X2 Dual Core Processor 3800+, 2 gigabytes RAM. The computers were more or less equivalent in the cited papers.Table 1
Computation times (in s) for different methods. ((b) from [32].)
PDB codeNo. atomsFDM1.9rs,0FDM4/3rs,1LSMSPyMolSwiss-PDBVChimeraMSMS1A8R264000.652.615.5610.526.3816.36(b)0.951H2I323180.721.836.5011.375.2540.04(b)3.031GTP347400.512.206.9813.154.7567.04(b)9.021FKA349770.752.777.8926.297.3677.25(b)4.501GT7427000.952.437.3216.106.5054.39(b)3.321GAV433350.552.547.0528.867.7178.35(b)4.221G3I455280.542.858.1819.456.21—7.671PMA458920.401.978.2318.676.72—12.901FJG519950.712.888.0125.188.05—15.111AON588700.622.648.8326.368.91—10.871J0B601440.693.076.8732.667.92—5.611OTZ686200.672.288.4630.149.56—9.031IR2770880.652.887.0929.319.5593.87(b)9.49Table 2
Computation times (in s) for different methods. ((a) from [30].)
PDB codeNo. atoms≃ No. trianglesFDMrs/2,1ChengMSMS200D23265 k0.551.35(a)0.331FG1873100 k0.852.41(a)0.653EBZ1651200 k1.3715.43(a)0.97Figure 10
Mesh of the SES of a molecule withTe=4/3rs and kr=1. (PDB Code: 1H2I).It appears in Tables1 and 2 that all the molecule surfaces in this data set are computed faster with the FDM than with any other method and for different values of the parameters.


### 4.3. Quality Results

In order to validate the quality of the results, different generated surfaces (SESs) were compared with references surfaces. These reference surfaces were generated by isosurfacing a field composed of a union of VdW balls at good resolution (spatial spacing ofrs/8) after morphological closing with a structuring element of the size of the solvent molecule. This approach, similar to [32], directly follows the definition of the SES [22], because morphologically close this volume is equivalent to make a solvent molecule roll on the VdW balls and to consider unaccessible parts to be inside the molecule. The molecules tested were 200D, 1FG1, and 3EBZ, because they are small enough to generate a good reference surface with the available memory space. The mean weighted root mean square deviation (RMSD) for three spatial spacing is reported in Table 3. The weighted RMSD isTable 3
RMSD and percentage of big differences with reference surfaces.
Spacing (Å)RMSDw (Å)∥·∥>1.9rs∥·∥>rs∥·∥>rs/41.9rs=2.660.780.41%6.89%67.04%rs=1.40.260%0%20.20%rs/4=0.350.200%0%8.71%(7)RMSDw=∑n=1N∥pnpn′∥2snS,
where N is the number of vertices in the reference mesh, pn is a vertex of this mesh, pn′ is the closest point on the other mesh (not necessary a vertex, it can be on an edge or on a face), sn is the mean area of the faces pn belongs to, and S is the total surface area. The percentages of the surface for which ∥pnpn′∥ is greater than the spacings are shown in the right-hand-side columns of Table 3. These error indexes are not completely correct because the reference surfaces is not a ground truth. However, it shows that the FDM algorithm can provide a surface with a quality comparable to robust methods.A visual comparison between the SES of 1FG1 computed with the FDM algorithm and the reference SES is shown in Figure11.Visual comparison between the 1FG1 SES generated with the FDM algorithm with a spacing ofrs/4 and the reference SES generated with a spacing of rs/8.
(a)(b)


## 4.1. Parameters

There are three main parameters modifiable by the user. First, the spatial spacingTe, that is, the distance between two neighbor voxels center, which determines the total number of voxels. Second, the cutoff frequency fc, which determines the smoothness of the final mesh. And third, the refinement rate kr, that is, the number of new points in a triangle for the final mesh magnification. In this section, the effect of these parameters on the visual quality, on the computation time and on the memory space, are discussed.


### 4.1.1. Spatial Spacing

With a small spatial spacing, it is possible to represent fine details. However, it drastically increases the memory space needed as well as the computation time. Indeed, reducingTe by a factor α increases the number of voxels by α3. The parts of the method depending on the number of voxels are the creation of the density map (time: 𝒪(nv) and space: 𝒪(nv), with nv the number of voxels) and the Fast Fourier Transform (time: 𝒪(nvlognv) and space: 𝒪(nv)). A visual comparison between meshes generated with Te=1.9rs, Te=rs and Te=rs/3 is shown in Figure 7. In these examples, fc=1/4rs and kr=0, that is, the meshes represent the SES without final refinement.Meshes generated from electron density map at a spatial spacing ofTe=1.9rs (a), Te=rs (b), and Te=rs/3 (c) (PDB code: 3EBZ).
(a)(b)(c)


### 4.1.2. Cutoff Frequency

In order to generate a mesh representing the SES,fc is set to 1/4rs (see Section 3.2). However, depending on the application, the surface could be other than the SES. For instance, if fc>fe/2, there is no actual filtering, and so, the generated mesh represents the VdWS. On the other hand, to obtain a smooth approximation of the molecule shape, fc can be reduced. It is equivalent to consider a bigger solvent molecule. Changing this parameter does not have any effect neither on the computation time nor on the memory space needed. A visual comparison between meshes generated with fc>fe/2 (VdWS), fc=1/4rs, and fc=1/16rs (shape approximation) is shown in Figure 8. In these examples, Te=rs/3 and kr=0, that is, a solvent molecule diameter takes 6 voxels and there is no final refinement.Meshes generated from electron density map filtered at a cutoff frequency offc>fe/2 (a) to give the VdWS, fc=1/4rs (b) to give the SES, and fc=1/16rs (c) to give an approximation of the general shape. (PDB code: 3EBZ).
(a)(b)(c)


### 4.1.3. Refinement Factor

The final mesh refinement gives foremost an esthetic advantage. The memory space needed does not increase a lot because the number of voxels remains the same. Only the size of the mesh changes and this is negligible in comparison with the space needed by the voxels representation. The computation time slightly increases but, whenkr=1 (which gives a good visual improving), this is negligible in comparison with the voxels operations computation time. A visual comparison between meshes generated with kr=0, kr=1 and kr=2 is shown in Figure 9. In these examples, Te=rs and fc=1/4rs, that is, a solvent molecule diameter takes 2 voxels and the meshes represent the SES.Meshes refined with a factor ofkr=0 (a), kr=1 (b), and kr=2 (c) (PDB code: 3EBZ).
(a)(b)(c)


## 4.1.1. Spatial Spacing

With a small spatial spacing, it is possible to represent fine details. However, it drastically increases the memory space needed as well as the computation time. Indeed, reducingTe by a factor α increases the number of voxels by α3. The parts of the method depending on the number of voxels are the creation of the density map (time: 𝒪(nv) and space: 𝒪(nv), with nv the number of voxels) and the Fast Fourier Transform (time: 𝒪(nvlognv) and space: 𝒪(nv)). A visual comparison between meshes generated with Te=1.9rs, Te=rs and Te=rs/3 is shown in Figure 7. In these examples, fc=1/4rs and kr=0, that is, the meshes represent the SES without final refinement.Meshes generated from electron density map at a spatial spacing ofTe=1.9rs (a), Te=rs (b), and Te=rs/3 (c) (PDB code: 3EBZ).
(a)(b)(c)


## 4.1.2. Cutoff Frequency

In order to generate a mesh representing the SES,fc is set to 1/4rs (see Section 3.2). However, depending on the application, the surface could be other than the SES. For instance, if fc>fe/2, there is no actual filtering, and so, the generated mesh represents the VdWS. On the other hand, to obtain a smooth approximation of the molecule shape, fc can be reduced. It is equivalent to consider a bigger solvent molecule. Changing this parameter does not have any effect neither on the computation time nor on the memory space needed. A visual comparison between meshes generated with fc>fe/2 (VdWS), fc=1/4rs, and fc=1/16rs (shape approximation) is shown in Figure 8. In these examples, Te=rs/3 and kr=0, that is, a solvent molecule diameter takes 6 voxels and there is no final refinement.Meshes generated from electron density map filtered at a cutoff frequency offc>fe/2 (a) to give the VdWS, fc=1/4rs (b) to give the SES, and fc=1/16rs (c) to give an approximation of the general shape. (PDB code: 3EBZ).
(a)(b)(c)


## 4.1.3. Refinement Factor

The final mesh refinement gives foremost an esthetic advantage. The memory space needed does not increase a lot because the number of voxels remains the same. Only the size of the mesh changes and this is negligible in comparison with the space needed by the voxels representation. The computation time slightly increases but, whenkr=1 (which gives a good visual improving), this is negligible in comparison with the voxels operations computation time. A visual comparison between meshes generated with kr=0, kr=1 and kr=2 is shown in Figure 9. In these examples, Te=rs and fc=1/4rs, that is, a solvent molecule diameter takes 2 voxels and the meshes represent the SES.Meshes refined with a factor ofkr=0 (a), kr=1 (b), and kr=2 (c) (PDB code: 3EBZ).
(a)(b)(c)


## 4.2. Time Comparisons

In this section, computation times are compared between the FDM algorithm and algorithms found in the literature for equivalent qualities. When available, the algorithms were run on the same computer, when not, the computation times were the ones announced in the original paper. Can et al. made a comparison of their method computation time with three molecular visualization tools: UCSF Chimera [37], Swiss-PDBViewer [38], and PyMol [9]. We added MSMS [26] to this set of methods. These programs, as well as the one of Can (LSMS), are available for free, so, the computation time could be measured on the same computer than for our method, except for Chimera that was not supported by the system. Thus, the computation times mentioned here for Chimera are the one announced in the paper of Can et al. [32]. For the LSMS method, the grid size is set to 256 × 256 × 256. In this condition, Te≃1 with the tested molecules. Other programs are run with defaults settings, that is, Te≃1. Two tests are made for the FDM method. In the first one, the parameters are set to Te=1.9rs and kr=0 to be as fast as possible while keeping a correct solution. In the second one, the parameters are set to Te=4/3rs and kr=1 which gives a correct mesh with a good appearance (see Figure 10). The computation times are shown in Table 1. The computation times reported only include the mesh generation time, that is, it does not take the loading time into account. In addition, for three molecules, computation times for the method of Cheng and Shi [30] are reported from their paper. The computation times comparison is shown in Table 2. All the tests were performed on a AMD Athlon(tm) 64 X2 Dual Core Processor 3800+, 2 gigabytes RAM. The computers were more or less equivalent in the cited papers.Table 1
Computation times (in s) for different methods. ((b) from [32].)
PDB codeNo. atomsFDM1.9rs,0FDM4/3rs,1LSMSPyMolSwiss-PDBVChimeraMSMS1A8R264000.652.615.5610.526.3816.36(b)0.951H2I323180.721.836.5011.375.2540.04(b)3.031GTP347400.512.206.9813.154.7567.04(b)9.021FKA349770.752.777.8926.297.3677.25(b)4.501GT7427000.952.437.3216.106.5054.39(b)3.321GAV433350.552.547.0528.867.7178.35(b)4.221G3I455280.542.858.1819.456.21—7.671PMA458920.401.978.2318.676.72—12.901FJG519950.712.888.0125.188.05—15.111AON588700.622.648.8326.368.91—10.871J0B601440.693.076.8732.667.92—5.611OTZ686200.672.288.4630.149.56—9.031IR2770880.652.887.0929.319.5593.87(b)9.49Table 2
Computation times (in s) for different methods. ((a) from [30].)
PDB codeNo. atoms≃ No. trianglesFDMrs/2,1ChengMSMS200D23265 k0.551.35(a)0.331FG1873100 k0.852.41(a)0.653EBZ1651200 k1.3715.43(a)0.97Figure 10
Mesh of the SES of a molecule withTe=4/3rs and kr=1. (PDB Code: 1H2I).It appears in Tables1 and 2 that all the molecule surfaces in this data set are computed faster with the FDM than with any other method and for different values of the parameters.


## 4.3. Quality Results

In order to validate the quality of the results, different generated surfaces (SESs) were compared with references surfaces. These reference surfaces were generated by isosurfacing a field composed of a union of VdW balls at good resolution (spatial spacing ofrs/8) after morphological closing with a structuring element of the size of the solvent molecule. This approach, similar to [32], directly follows the definition of the SES [22], because morphologically close this volume is equivalent to make a solvent molecule roll on the VdW balls and to consider unaccessible parts to be inside the molecule. The molecules tested were 200D, 1FG1, and 3EBZ, because they are small enough to generate a good reference surface with the available memory space. The mean weighted root mean square deviation (RMSD) for three spatial spacing is reported in Table 3. The weighted RMSD isTable 3
RMSD and percentage of big differences with reference surfaces.
Spacing (Å)RMSDw (Å)∥·∥>1.9rs∥·∥>rs∥·∥>rs/41.9rs=2.660.780.41%6.89%67.04%rs=1.40.260%0%20.20%rs/4=0.350.200%0%8.71%(7)RMSDw=∑n=1N∥pnpn′∥2snS,
where N is the number of vertices in the reference mesh, pn is a vertex of this mesh, pn′ is the closest point on the other mesh (not necessary a vertex, it can be on an edge or on a face), sn is the mean area of the faces pn belongs to, and S is the total surface area. The percentages of the surface for which ∥pnpn′∥ is greater than the spacings are shown in the right-hand-side columns of Table 3. These error indexes are not completely correct because the reference surfaces is not a ground truth. However, it shows that the FDM algorithm can provide a surface with a quality comparable to robust methods.A visual comparison between the SES of 1FG1 computed with the FDM algorithm and the reference SES is shown in Figure11.Visual comparison between the 1FG1 SES generated with the FDM algorithm with a spacing ofrs/4 and the reference SES generated with a spacing of rs/8.
(a)(b)


## 5. Conclusion

In this paper, we introduced an algorithm to compute molecular surface meshes (the FDM algorithm). It is constructed as an isosurface of a filtered electron density map (FDM). This algorithm is faster than other algorithm tested in equivalent conditions. It is slower than the MSMS algorithm for small molecules (<30000 atoms) but it returns a smooth manifold surface, which is not the case with MSMS. It makes possible to compute a precise representation of the surface with a limited number of voxels, so that the computation time and the memory space needed are reduced. Moreover, it is parameterizable on the spatial resolution, the refinement of the final mesh, and the size of the solvent molecule. Thus, the spatial resolution can be improved for a finer result but with an important computation time increase. Similarly, a smoother result can be obtained with a final refinement with a small influence on the computation time but with less precise results than reducing the spacing. Finally, the solvent molecule size can be chosen without influence on the computation time.The refinement could be improved to be specific to molecular surface. It would enable coarse meshes to be generated rapidly and to be improved by a priori knowledge about local geometry of molecule surfaces, such that the curvature deduced from the closest atom radius. In future works, this algorithm will be used in surface-based method to detect protein hot spots [19].



---

*Source: 923780-2010-04-12.xml*