PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Comput Sci Technol. Author manuscript; available in PMC 2010 May 20.
Published in final edited form as:
J Comput Sci Technol. 2008 November 1; 23(6): 1026–1036.
doi:  10.1007/s11390-008-9184-1
PMCID: PMC2873780
NIHMSID: NIHMS87564

Bio-molecule Surfaces Construction via a Higher-Order Level-Set Method

Abstract

We present a general framework for a higher-order spline level-set (HLS) method and apply this to bio-molecule surfaces construction. Starting from a first order energy functional, we obtain a general level set formulation of geometric partial differential equation, and provide an efficient approach to solve this partial differential equation using a C2 spline basis. We also present a fast cubic spline interpolation algorithm based on convolution and the Z-transform, which exploits the local relationship of interpolatory cubic spline coefficients with respect to given function data values. One example of our HLS method is demonstrated, which is the construction of bio-molecule surfaces (an implicit solvation interface) with their individual atomic coordinates and solvated radii as prerequisite.

Keywords: higher-order spline level-set, geometric partial differential equation, bio-molecular surfaces

1 Introduction

Given a non-negative function g(x) over a domain Ω [subset or is implied by] R3, find a surface Γ in Ω to minimize the energy functional

E(Γ)=Γg(x)dx+εΓh(x,n)dx,
(1)

where x and n are surface point and normal to the surface, respectively. Furthermore, h(x,n) is another non-negative function defined over R3 × R3\{0} chosen for regularizing the constructed smooth surface. Constant ε ≥ 0 is to control the regularization effect. Many solid and physical modeling problems, such as surface (solid boundary) reconstruction, and physics-based simulation of deformable interfaces could be formulated as minimizing one kind of energy in the form of (1). By calculating the variation of this energy functional, a partial differential equation (PDE) in the level-set formulation can be generated. In this paper, we propose a higher-order spline level-set method to solve the obtained PDE. To verify the effectiveness of this method, smooth surface constructions experiments are carried out. As illustrative examples, Figure 1 shows the smooth solvent excluded molecular surface construction.

Figure 1
(a) shows the van der Waals surface of a molecule. (b) shows the corresponding solvent excluded molecular surface constructed using our C2 tri-cubic spline level-set method. (c) illustrates that the smooth solvent excluded surface constructed encloses ...

Why Use a Level-Set Method?

In shape deformation simulations, topology changes may occur. This topology change makes parametric form surface tracking difficult even though the hard problem of global parametrization of a discrete surface could be successfully solved. However, implicit form surface (level-set surface) deformation could overcome this difficulty. Thus, level-set method has been a leading subject in many areas (the interested readers are referred to two monographs [1] and [2] for details and references) since the seminal work of Osher and Sethian ([3]) appeared. The level-set method allows one to dynamically deform and track an implicit surface using a governing PDE, which describes various laws of motion depending on geometry, external forces, or a desired energy minimization. Furthermore, the underlining data structure is simple with the computation being restricted to a narrow band surrounding the level-set.

Why Use a Higher-Order Spline Method?

The level-set surfaces obtained from classical level-set methods are generally bumpy due to the use of piecewise tri-linear interpolation from the discrete function data computed on a rectilinear grid. To produce a higher quality surface, a denser grid needs to be used. However, the increased grid resolution substantially increases the computation costs. Another drawback of using discrete data over grids is the non-trivial requirement of estimating derivatives for smooth interpolation. In many surface construction problems, such as the construction of molecular surface, the underlining surface is at least C1 smooth. Therefore, a smooth level-set function is highly desirable. In this paper, second-order geometric partial differential equations are required to be solved with mean curvature involved. Therefore, we use C2 tri-cubic spline as the level-set function basis. Note that tri-cubic is the lowest order of B-spline that could achieve C2 in 3D. The advantages using C2 spline function basis include:

  1. Since the level-set function is C2, the level-set surface is G2. There does exist a finite number of critical level-set values where the level-set surface may have a finite number of isolated singular points (i.e., the gradient of the level-set function vanishes). However working in a finite precision numerical domain one can automatically avoid this finite set of critical level-set values.
  2. Derivatives up to the second order and curvatures, which appear in the governing geometric partial differential equations, can be easily and accurately computed from the C2 level-set function.

Main Contributions

We derive and describe an efficient approach of higher-order local level-set method to solve PDEs. Using a first order energy functional, a second order geometric PDE in the level-set formulation is obtained. To achieve higher efficiency in using tri-cubic spline functions for solving PDEs, a fast and exact cubic spline interpolation algorithm is presented. Furthermore, the local property of cubic spline function is analyzed. Based on this analysis, a local spline level-set method is implemented. Though we use tri-cubic spline functions here, we could easily use splines of any other orders as well. We apply our method to the smooth construction of bio-molecular solvation interfaces. Our construction method is simple, efficient and error bounded.

The rest of the paper is organized as follows. Section 2 outlines the main algorithm steps for solving the geometric PDE. The details of these steps are considered in the sections that follow. In Section 3, we describe a fast spline interpolation algorithm. The derivation of the interpolation algorithm is given in our technical report [4]. The initialization, evolution and re-initialization steps in the main algorithm are discussed in Section 4. Application to bio-molecular surfaces construction is considered in Sections 5. Section 6 concludes the paper.

2 Algorithm Outline

By calculating the variation of energy functional (1), we obtain the following evolution equation (see [5])

φt=(g+εh)div(φφ)φ+εdiv(Pnh)φ+2[(g+εh)]Tφ=L(φ)+H(φ),
(2)

where

L(φ)=(g+εh)div(φφ)φ,H(φ)=εdiv(Pnh)φ+2[(g+εh)]Tφ,

P=Iφφφφ is a projection operator onto the tangent space and I indicates the identity mapping. [nabla] and [nabla]n denote the usual gradient operator with respect to x and n, respectively. Note that L([var phi]) is a parabolic term and H([nabla][var phi]) is a hyperbolic term. Hence, in solving equation (2) in the following, the first order term H([nabla][var phi]) is computed using an upwind scheme (see [6] for the reason) over a finer grid, and the higher order term L([var phi]) is computed using a spline presentation defined on a coarser grid. In the experiments carried out in this paper, we always choose h(x,n) to be 1 for simplicity.

Consider the solution of equation (2) over the domain Ω = [a,b] × [c,d] × [e,f] [subset or is implied by] R3. For simplicity, we assume b − a = d − c = f − e > 0 and the domain Ω is uniformly partitioned with vertices G0={xijk}ijk=0n:={xi}i=0n×{yj}j=0n×{zk}k=0n, where

xi=a+iΔx,yj=c+jΔy,zk=e+kΔz,

and Δx = Δy = Δz = (b − a)/n. Let Gl be the set of grid points generated by bisecting G0 l times. Let [var phi] be a piecewise tri-linear level-set function over the grid Gl, Φ be a tri-cubic spline approximation of [var phi] over the grid G0. In general, l is chosen to be 0 or 1 or 2. In our implementation, we take l=1. If l = 0, Φ and [var phi] are defined on the same grid G0, which is the simplest case. The aim of the following algorithm is to compute the spline level set function Φ.

Algorithm Steps

  1. Initialization. Given an initial Γ, construct a piecewise tri-linear level-set function [var phi] over the grid Gl. If necessary, apply a reinitialization step to set [var phi] to be a signed distance function to Γ (see section 4.2 for details). Convert [var phi] to Φ (see Section 3).
  2. Evolution. Resample Φ to obtain a new [var phi] over the grid Gl. Compute L(Φ) and H([nabla][var phi]) in the thin shell (traditionally called a narrow band for curve evolution) N = {(xi,yj, Zk) [set membership] Gl : |[var phi](xi,yj,zk)| < H}. Update [var phi] in N for one time step to get [phi] by an ODE time stepping method (see section 4.3).
  3. Re-initialization. Apply re-initialization step to [phi] in the shell
    N={(xi,yj,zk)Gl:min1μ,ν,λ1|φi+μ,j+ν,k+λ|}
    to get a new [var phi] (see section 4.4). Convert [var phi] to Φ (see Section 3). Go back to step 2 if the termination condition is not satisfied.
  4. Iso-contouring. Extract 3-sided or 4-sided iso-surface patches (vertices with normals) of Φ = c, where c is a given iso-value. G1 approximate these patches by parametric surfaces.

Remark 2.1

For the problem of molecular surface construction, the grid size G0 should be less than the radii of atoms so that atoms are distinguishable from the level set surface. In our implementation, the grid size is chosen to be one-half of the minimal value of the atom radii.

Remark 2.2

The aim of using l > 0 is to make [var phi] a more accurate approximation of the signed distance function. The larger the value of l we use, the better approximation of the signed distance function we have. Since the scanned data to be approximated in general suffers from noise, we use the approximation Φ over a coarse grid G0 for denoising.

3 Fast Cubic Spline Interpolation

3.1 Algorithm

Let

β3(x)={23x2+12|x|3,0|x|<1,16(2|x|)3,1|x|<2,0,2|x|.

be the cubic B-spline base function over the knots {−2,−1,0,1,2}, and let

s(x)=i=1n1ciβ3(xi)
(3)

be a cubic spline function. Then the spline interpolation problem is to determine the coefficients c={ci}i=1n1 such that the following interpolation conditions

s(k)=i=1n1ciβ3(ki),k=1,,n1
(4)

are satisfied for any given function values {s(k)}k=1n1. Such a problem could be easily solved by solving the linear system (4) directly for the unknowns c1, …,cn−1. The algorithm derived here is in O(n) complexity. Given two function values s(0) and s(n), system (4) could be augmented as

s(k)=i=1n1ciβ3(ki),k=0,,n.

Recursive Algorithm

Using the initial values

c0+=11z12nk=1n1s(k)z1k(1z12(nk)),cn=0,
(5)

the recursive process

ck+=s(k)+z1ck1+,k=1,2,,n1,ck=z1(ck+1ck+),k=n1,n2,,1,

yields exact B-spline coefficients ci=6ci of the interpolation problem, where z1=32.

Remark 3.1

The computation of c0+ can be accelerated by neglecting small terms. Given an error control tolerance ε, we compute

K=[log(ε)log(|z1|)].

If n > K, we replace c0+ in (5) with

c0+=11z12nk=1n1s(k)z1k,

which can be computed by the Horner scheme for evaluating a polynomial.

3.2 Locality of Spline Interpolation

Since the basis function β3 (x) is locally supported, the kth coefficient ck of the spline function (3) merely has influence on s(x) within the interval (k − 2, k + 2). A set of function values {s(i)}

s(i)={0,ik,1,i=k,

define a set of coefficients {ci} with none of them zero. However, the coefficients ck±l are decreasing with respect to l. To see this, let us compute ck. Suppose that the interpolation problem considered is over (−∞,+∞). From the definition of ck, we have

ck=z1(1+z12+z14+)=z11z12=36.

Similarly,

ck±l=z1l+1(1+z12+z14+)=z1l36,

for l = 0,1,2, …. Therefore,

ck±l=z1l3,l=0,1,2,....

Note that z1=320.2679491924 and ck±l will decrease to zero rapidly. The first ten terms are 1.732, −0.464, 0.124, −0.0333, 0.00893, −0.00239, 0.000641, −0.000172, 0.0000460, −0.0000123. Given a threshold ε, the coefficients ck±l with |ck±l| < ε could be ignored. We take ε = 0.005. Therefore, we keep the terms

ck±l,l=0,1,2,3,4.

If we require that the second order derivatives have accuracy ε, we need to determine l, such that

|z1|l3<Δx2ε,
(6)

where Δx is the spacing of the knots.

Remark 3.2

Higher dimensional spline interpolation can be recursively computed by using lower dimensional spline computations. We do not exposit this here.

3.3 Conversion of Piecewise Tri-linear Functions to Tri-cubic Splines

Let [var phi]l be a piecewise tri-linear function over the grid Gl. Now we intend to convert it approximately to a trivariate tri-cubic spline function Φ over the grid G0. If l = 0, the conversion algorithm is described in section 3.1. Now we assume l > 0. The conversion could be done in several ways. One simple way is to use Φ to interpolate [var phi]l over the grid Gl in the least square sense. This requires us to solve a linear system with size n3, which is the number of grid in G0. Another algorithm is to convert [var phi]l to a spline function Φl over the grid Gl firstly, and then use the knots removal algorithm of spline function to obtain Φl−1, Φl−2, …, Φ0 repeatedly. Since Φk may not be represented exactly by Φk−1, a least square approximation has to be used. This also leads to solving a linear system.

The method we use in this paper is to convert [var phi]l to [var phi]l−1, [var phi]l−2, …, [var phi]0, where [var phi]k is a piecewise tri-linear function over the grid Gk. Finally, we convert [var phi]0 to Φ using the conversion algorithm described in section 3.1. For converting [var phi]k to [var phi]k−1, we use a local tri-linear function at each grid point of Gk−1 to interpolate [var phi]k at 27 neighbor grid points of Gk in the least square sense. This leads to an 8×8 linear system of equations for each of the grid point of Gk−1. Note that all these systems have the same coefficient matrix. Hence, we only need to inverse the matrix once. Furthermore, this conversion are conducted only in a neighborhood of the level-set surface (see Section 4). Therefore it is very efficient.

4 Solving the PDE by a Local Level-set Method

In this section, we describe a local level-set method, using a tri-cubic spline level-set function Φ over the grid G0 and a piecewise tri-linear level-set function [var phi] over the grid Gl, which represents a signed distance function. We first define a thin shell around the interface in which the level-set function is updated. Then we determine the initial level-set function followed by the evolution of the level-set function. Finally we explain the reinitialization step.

4.1 Thin Shell of a Tri-cubic Spline

During the evolution process, we confine the movement of the level-set surface to no more than one grid size Δx of Gl for each time step. Since ‖[nabla][var phi]‖ = 1 before the movement, we require that function values, and up to the second order partial derivatives of [var phi], are accurately computed (within the given error tolerance) in the c + Δx neighborhood of the interface, where c is the iso-value used for extracting the level-set (see the last step of main algorithm in Section 2). From (6), we define the thin shell as

N={xGl:|φ(x)|},

where

=[log(Δx)2εlog3log|z1|+2]Δx+|c|.
(7)

As mentioned before, we take ε = 0.005. Function values that are beyond H are truncated by H, meaning if [var phi] (x) > H, then set [var phi] (x) as H, if [var phi](x) < H, then set [var phi](x) as H. The evolution of [var phi] is performed in this thin shell. The time step is so determined that the function values of [var phi] are changed not more than Δx. By using an explicit Euler method, this is easy to achieve. After the evolution step, the function [var phi] is truncated by H.

In the re-initialization step, the thin shell is updated by introducing new grid points that are the neighbor points of the previous thin shell, and furthermore by removing the points whose values are greater than H in magnitude.

4.2 Initialization

Let Γ0 be a given initial closed surface with interior An external file that holds a picture, illustration, etc.
Object name is nihms87564ig1.jpg [subset or is implied by] R3. We define an initial level-set function [var phi]0(x) in R3 satisfying

{φ0(x)<0inD,φ0(x)=0onΓ0,φ0(x)>0elsewhere.
(8)

Γ0 may be given as a close triangular surface mesh, or a level-set of a given function. For each of these two cases, the initial [var phi]0 needs to be computed differently.

(a) Γ0 is a surface mesh

If Γ0 is a surface mesh, we first compute the Euclidean distance from the grid points to the surface mesh using a local fast method (see Algorithm 4.1), and then we specify negative/positive signs to the distance values for the grids inside/outside the surface (see Algorithm 4.2). Since the function values need to be exact only in the thin shell, we compute distance values for the grid points locally around each triangle within the distance H. This is done by the following algorithms.

Algorithm 4.1. Distance Computation
  1. Set initial value for each grid point to be H.
  2. For each triangle of the mesh, find a containing cube such that the distance of the cube boundary to the triangle is no less than H. For each grid point in the cube, compute the distance from the grid point to the triangle. If this distance is less than the existing value previously computed, then replace it with the newly computed one.

Algorithm 4.2. Sign Designation
  1. Perturb the vertices of the given triangular surface mesh to meet the following two requirements. (i) The x-components of all the face normals are not zero; (ii) The projections of all the edges of the mesh onto the yz-plane do not pass through any grid point on this plane.
  2. Compute distance from the grid points to the (perturbed) surface mesh using Algorithm 4.1.
  3. First we associate each grid point q{qjk}jk=0n:={yj}j=0n×{zk}k=0n on the yz-plane with a set of numbers, called q-pool. This pool is empty at beginning. Let Δx be the projection of the triangle Δ of the surface mesh onto the yz-plane. For each grid point q [set membership] Δx, compute the intersection of the line segment [x,q]T [set membership] R3, x [set membership] [a,b] with the triangle Δ. Put the x-component of the intersection point into q-pool.
  4. After the computation above, we have a pool for each grid point in the yz-plane. Then we sort the numbers in each of the pool in the increasing order. Let t1 < t2 < … < tk be the sequence in the q-pool. Then if a + jΔx is in the interval [ti,ti+1], the sign to the grid point [a+jΔx, q]T is set to be (−1)i, where t0 = a, tk+1 = b.

Remark 4.1

The aim of the perturbation in the first step of Algorithm 4.2 is to make the intersection of the line segment [x,q]T, x [set membership] [a,b] with Δ to be unique and there is no duplicated numbers in each of q-pool. This perturbation is not harmful to the surface construction algorithm, since Γ0 is served as an initial surface.

Remark 4.2

The time complexities of both Algorithm 4.1 and 4.2 are linear with respect to the number of triangles of the given mesh.

(b) Γ0 is a level-set

If Γ0 is a level-set {x : d(x) = c} of a smooth function d(x), a naive method is to first extract the level-set surface and then compute a signed distance function using the method described earlier. Since extracting iso-surface is not a trivial task, we seek a more efficient method. Suppose that [nabla]d(x) ≠ 0, x [set membership] Γ0. Then a simple way is to use d(x) − c or c − d(x) as the initial [var phi]0. For example, in the problem of solvent excluded surface construction, d(x) is a Gaussian map, c=1 (see Section 5).

4.3 Evolution

Suppose that we have an initial [var phi]0 satisfying (8). First we apply the local re-initialization step, if necessary, to set [var phi]0 to be the signed distance function of Γ0 (see section 4.4). Then we define a thin shell around Γ0 by

N0={xGl:|φ0(x)|<}.

To prevent numerical oscillations at the boundary of the thin shell, we introduce a cubic cut-off function c(x) in (9), which is defined by C1 Hermite interpolating the boundary data on [H0,H]

c(x)={1if|x|0,(|x|)2(2|x|+30)(0)3if0<|x|0if|x|>,,

where H0 is chosen to be 0.5H. Our experiments show that linear cut-off function works equally well. We update [var phi]0 by solving the following equation,

{φt=c(φ)[L(φ)+(φ)],φ(x,0)=φ0(x),
(9)

on N0 for one time step and get [var phi]1(x). The time step is chosen such that the interface moves less than one grid size Δx. At each grid point xijk in the thin shell N0, compute ν0(xijk) = c([var phi]0(xijk))[L0(xijk)) + H([nabla][var phi]0(xijk))]. Let

τ=min{Δx,Δx/maxxijkN0|ν0(xijk)|}.

Then update [var phi]0 by the explicit Euler scheme

φ1(xijk)=φ0(xijk)+τν0(xijk),xijkN0.

Since [var phi]1 is no longer a signed distance function, a re-initialization step is required (see section 4.4) to get a new [var phi]1 and a new thin shell N1. The process from [var phi]0 to [var phi]1 described above is repeated to get a sequence {[var phi]m}m≥0 of [var phi], and a sequence of thin shells {Nm}m≥0, until the following termination conditions

maxxijkNm|νm(xijk)|<εandm<M
(10)

are satisfied. We choose ε = 0.001. M is a given upper bound of the iteration number. We choose M = n, where n3 is the number of grid points in Gl.

Remark 4.3

The solution to equation (2) may have discontinuous derivatives, even if the initial data is smooth (see [6]). Using simple central difference to approximate the first order term is not appropriate. Hence, we use monotone and upwind scheme for Hamilton-Jacobi equations as developed in [6] for computing H([nabla][var phi]m). The higher order term L([var phi]m) is computed from the spline function Φm directly. Four monotone discretization schemes (Lax-Friedrichs scheme, Godunov type scheme, Local Lax-Friedrichs scheme and Roe with LLF entropy correction) have been reviewed in [6]. We use Godunov's scheme because it is also upwinding.

4.4 Adaptive Re-initialization

To simplify the notation, let [var phi]m(x) be denoted by [var phi] (x). In general, it is impossible to prevent [var phi] (x) from deviating away from a signed distance function. Flat and/or steep regions will develop around the interface, making further computation and contour plotting highly inaccurate. Hence a re-initialization step to reset the level-set function [var phi] (x) to be a signed distance function to the interface {x : [var phi] (x) = 0} is necessary.

Re-initialization here is a process of replacing [var phi] (x) by another function [phi] (x) that is the signed distance function to the zero level-set of [var phi] and then taking this new function [phi] (x) as the initial data to use until the next round of re-initialization. To achieve this goal, people usually solve the following equation (see [7]) for the unknown function d (x,τ)

{dτ+S(d)(d1)=0,d(x,0)=φ(x),
(11)

where S(d) is a sign function. We propose to solve the following Hamilton-Jacobi equation,

{dτ+S(d)(dTd1)=0,d(x,0)=φ(x).
(12)

for its steady state solution, with S(d) approximated by

s=dd2+|Dd|2Δx2,
(13)

where Dd is a discrete approximation of [nabla]d. The advantages of using (12) instead of (11) is that the right hand side of (12) is a smooth function of d, and it also facilitates the usage of semi-implicit temporal discretization scheme.

Discretization

Again, we use the Godunov scheme (see [6]) for the spatial discretization because it is monotone and upwinding. The scheme is described abstractly in [6] for a 2D problem. Here we present the details for our 3D problem (12)

dijkn+1=dijknΔτsijk+(max[(a+)2,(b)2]+max[(c+)2,(d)2]+max[(e+)2,(f)2]1)Δτsijk(max[(a)2,(b+)2]+max[(c)2,(d+)2]+max[(e)2,(f+)2]1),
(14)

where Sijk is the approximation to S(dijk) with (13), (x)+ = max(x,0), (x) = min(x,0); a, …, f are defined by

a=Dxdijkn,c=Dydijkn,e=Dzdijkn,b=Dx+dijkn,d=Dy+dijkn,f=Dz+dijkn,

respectively, where Dx±dijk,Dy±dijk and Dz±dijk denote the one-sided divided differences at the grid point xijk in the x, y and z directions, separately. These one-sided differences are computed by ENO (essentially non-oscillatory) interpolation (see [6]). Scheme (14) is easily seen to be monotone; i.e., the right hand side of (14) is a nondecreasing function of all the dijkn if

ΔτΔx|sijk|13.

With the same time step restriction, this scheme is also upwind.

Adaptive Updating the Thin Shell

During the process of re-initialization, function values around the level-set change, and so the thin shell defined by {x [set membership] Gl : |[var phi](x)| < H} should be updated correspondingly. This is automatically achieved by including new grid points with at least one of their one-ring neighbor grid points with value less than H. Additionally we throw away the old grid points in the current thin shell where the function values and function values of all their one-ring neighbor grid points are all greater than H.

5 Smooth Bio-molecular Surfaces Construction

In this section, we use our higher order level-set method with C2 tri-cubic spline functions to construct bio-molecular interfaces. As is well known, there are typically three types of molecular surfaces (e.g. [8]): the van der Waals surface (VWS), the solvent-accessible surface (SAS), and the solvent-excluded surface (SES) or Connolly surface ([9]). The van der Waals surface is defined from the van der Waals radii of the atoms, which is the boundary of the region formed by the union of all the spheres (atoms). The SAS introduced by Lee-Richards is defined to be the locus of the center of the rolling sphere (always water molecule) which makes contact with the VWS surface ([10]). Hence the SAS is an inflated VWS with a probe radius of 1.4 Angstroms. The solvent-excluded surface is the solvent (sphere) surface inside of which the probe sphere never intrudes. In other words, SES is the offset surface of SAS in the inward direction with the solvent probe radius as the offset radius. This kind of surface can be represented by Alpha shapes ([11, 12, 13]). The Alpha shape technique has been extensively used for molecular surface modeling (e.g. [14]), area and volume computation and cavity and pocket recognition (e.g. [11, 15]). This technique requires a complex geometric data structure.

Gaussian blurred maps i=1med(xxi2ri2)/r2 have been frequently used in the molecular surface construction (see [16, 17, 18, 19]). Both VWS and SAS can be easily approximated by the level-set of the Gaussian map within any given tolerance by properly choosing the parameter d. However, the Gaussian map method cannot generate an accurate approximation to the SES (see Figure 5(a)). It is well known that SES is an important molecular surface (dielectric interface) for many applications, such as electrostatic energy computation (see [20]), various simulations (see [21]). Hence, we focus our attention in this section on the computation of SES using our HLS method.

Structural model of molecule M consisting of a collection of atoms with centers {xi}i=1m and radii {ri}i=1m (see Figure 3(a)) is retrievable from the Protein Data Bank (PDB). To construct the SES for M, we minimize the energy functional

Figure 3
(a) and (d) show the van der Waals surface (VWS). (b) and (e) are the corresponding solvent excluded surface (SES) constructed. (c) and (f) illustrate the tight enclosure of the smooth SES to the VWS.
E(Γ)=Γg(x)2dx+εΓdx,
(15)

where

g(x)=1i=1meC(xxi2ri2),ri=ri+rb,

with rb the probe radius, which is usually 1.4 Å. The constant C > 0 is so determined that g(x) = 0 is an approximation of the solvent accessible surface within a given tolerance (see Figure 3(b) and section 5.1). The second term of (15) is used to regularize the constructed surface, where ε is a small number. In the examples provided in the following, we choose ε as 0.1. The corresponding level-set formulation of the evolution equation for the energy (15) is

φt=(g2+ε)div(φφ)φ+4g(g)Tφ=L(φ)+H(φ),
(16)

where

L(φ)=(g2+ε)div(φφ)φ,H(φ)=4g(g)Tφ.

As before, the first order term H([nabla][var phi]) is computed using an upwind scheme over a finer grid, the higher order term L([var phi]) is computed using a spline presentation defined on a coarser grid. If [var phi] is a signed distance function and a steady solution of equation of (16), then the iso-surface [var phi] = − rb is an approximation of SES (see Figure 3(c)). The algorithm in Section 2 could be specialized as follows.

Algorithm 5.1. Solvent Excluded Surface Construction

  1. Compute g(x) over the grid Gl (see section 5.2 for the fast computation of g(x)).
  2. Compute an initial [var phi] by taking [var phi](x) = g(x) and then update it with a re-initialization step, such that ‖[nabla][var phi]‖ = 1 (see section 4.4).
  3. Update [var phi] by solving equation (16) for one time step using a divided difference method (see Section 4).
  4. Re-initialize [var phi], and then go back to the previous step if the stopping criterion (10) is not satisfied.
  5. Generate iso-surface {x : Φ(x) = −rb}, which is the required approximation of SES.

Remark 5.1

In step 2, [var phi](x) may change rapidly, we may slow down the changes of [var phi] by dividing it with a constant before taking the re-initialization step. Since the gradient of eC[xxi2ri2] is 2CeC[xxi2ri2][xxi,yyi,zzi]T, the length of the gradient is approximately 2Cri at the sphere surface {x:xxi=ri}. Hence we normalize the function [var phi] by dividing it with 2Cmaxiri.

5.1 Error Bounded Approximation

Let ε be a given error tolerance for the constructed SES. We determine the constant C to satisfy this error constraint. Assume that at each surface point of SAS, at most k atomic spheres intersect. Each of these spheres thus has a maximum contribution of 1/k to the Gauss map g(x) at such surface points of the SAS. That is

eC[xxi2ri2]=eC[(ri+ε)2ri2]=1k.

From this we cobtain

Clnk2riε.

For instance, if k is 4, we choose ε = 0.01, and C is thus chosen to be 50ln4rb+miniri.

5.2 Fast Computation of Gaussian Map

Since each term eC[xxi2ri2] in g(x) decreases rapidly as ‖xxi‖ increases, we evaluate eC[xxi2ri2] locally around xi. Suppose that we evaluate it for x such that

eC[xxi2ri2]>ε,

where ε is a given threshold value (e.g., we take ε = 10−5). We have

xxi2<ri2lnεC.
(17)

Hence, we only evaluate the term over the grids within the spherical range defined by (17). For the simplicity of implementation, we determine a minimal bounding cube containing the sphere defined by (17), and evaluate the function within the cube. The time complexity of this algorithm is obviously linear with respect to the number of atoms.

5.3 Illustrative Examples

We have implemented our higher-order level-set algorithm in the molecular visualization software package TexMol ([22]). We present a few examples to demonstrate the quality of the SES for molecules from the PDB. In these examples, the solvent probe (water molecule) radius is chosen to be 1.4. Figure 4 shows two molecular surface constructions, where figures (a) and (d) show the VWS, respectively, (b) and (e) are the corresponding SES. (c) and (f) illustrate how the smooth SES encloses the VWS tightly and how the SES transits between the atoms (spheres) smoothly.

Figure 4
Comparing the results of Gaussian blurring and our new method. (a) and (b) show the result of Gaussian blurring with d = −2.3. (c) and (d) are the level-set surfaces of [var phi] = −1.4 of the level set method. (b) and (d) show cross-section ...

Figure 5 shows the difference between our method with Gaussian blurring method (e.g. [19]) for an enzyme with 1UDI as its PDB ID. The surface of Gaussian blurring is defined by {x [set membership] R3 : g(x) = 1}, where

g(x)=i=1nedri2[xxi2ri2].

The surface generated by Gaussian blurring reveals creases (see Figure 5(a)). These creases should be covered by the rolling solvent probe spheres. The cutoff results further show that a lot of redundant interior structures exist (see Figure 5(b)). Our method gives more accurate results to the final SES (see Figure 5(c) and (d)).

6 Conclusions

We have proposed a general framework of higher-order spline level-set method for surface construction that could be used to solve smooth surface construction problems. We applied our method to the construction of C2 smooth solvent excluded surfaces of bio-molecules. Compared with Gaussian blurring technique, our method yields even smoother solvent excluded surfaces with a tighter enclosure of the atomic sphere models. Our HLS's generalization of the linear order level set method, is of course applicable to several other smooth surface construction applications [1] where C1 and higher-order smoothness is both necessary and/or desirable.

Figure 2
Schematic illustration of the construction algorithm. Figure (a) shows the van der Waals surface for a simple molecule of four atoms. (b) shows the solvent accessible surface defined by [var phi] = 0. (c) shows the solvent excluded surface defined by ...

Acknowledgments

We wish to thank Jianguang Sun and Vinay Siddavanahalli for all their help with the the software package TexMol, and Albert Chen for providing us with the PDB/PQR data of enzymes. A substantial part of this work in this paper was done when Guoliang Xu was visiting C. L. Bajaj at UT-CVC. His visit was supported by the J. T. Oden ICES visitor fellowship.

* Bajaj is supported in part by NSF grants IIS-0325550, CNS-0540033 and NIH contracts P20-RR020647, R01-GM074258, R01-GM07308. Xu and Zhang are supported by NSFC grant 10371130 and National Key Basic Research Project of China (2004CB318000)

References

1. Osher S, Fedkiw R. Level Set Method and Dynamic Implicit Surfaces. Springer; New York: 2003.
2. Sethian JA. Cambridge Monographs on Applied and Computational Mathematical. second. Cambridge University Press; 1999. Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Visison, and Materials Science.
3. Osher S, Sethian J. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics. 1988;79:12–49.
4. Bajaj CL, Xu G, Zhang Q. ICES Report 06-18. Institute for Computational and Engineering Sciences, The University of Texas at Austin; 2006. Smooth Surface Constructions via a Higher-Order Level-Set Method.
5. Xu G, Zhang Q. Construction of geometric partial differential equations in computational geometry. Mathematica Numerica Sinica. 2006;28(4):337–356.
6. Osher S, Shu CW. High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations. SIAM Journal of Numerical Analysis. 1991;28(4):907–922.
7. Peng D, Merriman B, Osher S, Zhao H, Kang M. A PDE-based fast local level set method. Journal of Computational Physics. 1999;155:410–438.
8. Sanner M, Olson A, Spehner J. Reduced surface: An efficient way to compute molecular surfaces. Biopolymers. 1996;38:305–320. [PubMed]
9. Connolly M. Analytical molecular surface calculation. J Appl Cryst. 1983;16:548–558.
10. Richards FM. Areas, volumes, packing, and protein structure. Ann Rev Biophys Bioeng. 1977;6:151–176. [PubMed]
11. Liang J, Edelsbrunner H, Fu P, Sudhakar PV, Subramaniam S. Analytical shape computation of macromolecules: I. molecular area and volume through Alpha shape. Proteins: Structure, Function, and Genetics. 1998;33:1–17. [PubMed]
12. Cazals F, Proust F. Revisiting the description of Protein-Protein interfaces. Part I: Algorithms. Research Report 5346, INRIA. 2004
13. Cazals F, Proust F. Revisting the description of Protein-Protein interfaces. Part II: Experimental study. Research Report 5501, INRIA. 2005
14. Akkiraju N, Edelsbrunner H. Triangulating the surface of a molecule. Discr Appl Math. 1996;71:5–22.
15. Liang J, Edelsbrunner H, Fu P, Sudhakar PV, Subramaniam S. Analytical shape computation of macromolecules: Ii. inaccessible cavities in protein. Proteins: Structure, Function, and Genetics. 1998;33:18–29. [PubMed]
16. Blinn J. A generalization of algebraic surface drawing. ACM Transactions on Graphics. 1982;1(3):235–256.
17. Duncan BS, Olson AJ. Shpae analysis of molecular surfaces. Biopolymers. 1993;33:231–238. [PubMed]
18. Grant J, Pickup B. A Gaussian description of molecular shape. Journal of Phys Chem. 1995;99:3503–3510.
19. Zhang Y, Xu G, Bajaj C. Quality meshing of implicit solvation models of biomolecular structures. Computer Aided Geometric Design. 2006;23(6):510–530. [PMC free article] [PubMed]
20. Baker N, Sept D, Joseph S, Holst M, McCammon J. Proc Natl Acad Sci. USA: 2001. Electrostatics of nanosystems: application to microtubules and the ribosome; pp. 10037–10041. [PubMed]
21. Holst M, Saied F. Multigrid solution of the Poisson-Boltzmann equation. J Comput Chem. 1993;14:105–113.
22. Bajaj CL, Djeu P, Siddavanahalli V, Thane A. Texmol: interactive visual exploration of large flexible multi-component molecular complexes. Proceedings of the Annual IEEE Visualization Cnference'04; Austin, Texax, USA. 2004. pp. 243–250.
23. Zhao H, Osher S, Merriman B, Kang M. Implicit nonparametric shape reconstruction from unorganized points using a variational level set method. Computer Vision and Image Understanding. 2000;80(3):295–319.
24. Zhao H, Osher S, Fedkiw R. Fast surface reconstruction using the level set method. Proceedings of the IEEE Workshop on Variational and Level Set Methods in Computer Vision (VLSM2001); 2001. pp. 194–201.