Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2873780

Formats

Article sections

- Abstract
- 1 Introduction
- 2 Algorithm Outline
- 3 Fast Cubic Spline Interpolation
- 4 Solving the PDE by a Local Level-set Method
- 5 Smooth Bio-molecular Surfaces Construction
- 6 Conclusions
- References

Authors

Related links

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-1PMCID: PMC2873780

NIHMSID: NIHMS87564

See other articles in PMC that cite the published article.

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 *C*^{2} 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.

Given a non-negative function *g*(**x**) over a domain Ω ^{3}, find a surface Γ in Ω to minimize the energy functional

$$E(\Gamma )={\int}_{\Gamma}g\left(\mathbf{\text{x}}\right)d\mathbf{\text{x}}+\epsilon {\int}_{\Gamma}h\left(\mathbf{\text{x}},\text{n}\right)d\mathbf{\text{x}},$$

(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 ^{3} × ^{3}\{**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.

(a) shows the van der Waals surface of a molecule. (b) shows the corresponding solvent excluded molecular surface constructed using our *C*^{2} tri-cubic spline level-set method. (c) illustrates that the smooth solvent excluded surface constructed encloses **...**

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.

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 *C*^{1} 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 *C*^{2} tri-cubic spline as the level-set function basis. Note that tri-cubic is the lowest order of B-spline that could achieve *C*^{2} in 3D. The advantages using *C*^{2} spline function basis include:

- Since the level-set function is
*C*^{2}, the level-set surface is*G*^{2}. 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. - 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
*C*^{2}level-set function.

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.

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

$$\begin{array}{l}\frac{\partial \phi}{\partial t}=\left(g+\epsilon h\right)\text{div}\phantom{\rule{0.2em}{0ex}}\left(\frac{\nabla \phi}{\Vert \nabla \phi \Vert}\right)\phantom{\rule{0.2em}{0ex}}\Vert \nabla \phi \Vert \hfill \\ \phantom{\rule{2.8em}{0ex}}+\epsilon \phantom{\rule{0.2em}{0ex}}\text{div}\phantom{\rule{0.1em}{0ex}}\left(\text{P}{\nabla}_{n}h\right)\Vert \nabla \phi \Vert +2{\left[\nabla \left(g+\epsilon h\right)\right]}^{T}\nabla \phi \hfill \\ \phantom{\rule{1.7em}{0ex}}=\phantom{\rule{0.2em}{0ex}}L\left(\phi \right)+H\left(\nabla \phi \right),\hfill \end{array}$$

(2)

where

$$\begin{array}{l}\phantom{\rule{1.2em}{0ex}}L\left(\phi \right)=\left(g+\epsilon h\right)\text{div}\phantom{\rule{0.2em}{0ex}}\left(\frac{\nabla \phi}{\Vert \nabla \phi \Vert}\right)\phantom{\rule{0.2em}{0ex}}\Vert \nabla \phi \Vert ,\hfill \\ H\left(\nabla \phi \right)=\epsilon \phantom{\rule{0.2em}{0ex}}\mathit{\text{div}}\left(\text{P}{\nabla}_{n}h\right)\Vert \nabla \phi \Vert +2{\left[\nabla \left(g+\epsilon h\right)\right]}^{T}\nabla \phi ,\hfill \end{array}$$

$\text{P}=\text{I}-\frac{\nabla \phi}{\Vert \nabla \phi \Vert}\otimes \frac{\nabla \phi}{\Vert \nabla \phi \Vert}$ is a projection operator onto the tangent space and I indicates the identity mapping. and * _{n}* denote the usual gradient operator with respect to

Consider the solution of equation (2) over the domain Ω = [*a,b*] × [*c,d*] × [*e,f*] ^{3}. For simplicity, we assume *b − a = d − c = f − e* > 0 and the domain Ω is uniformly partitioned with vertices
${G}_{0}={\{{\mathbf{\text{x}}}_{\mathit{\text{ijk}}}\}}_{\mathit{\text{ijk}}=0}^{n}:={\{{x}_{i}\}}_{i=0}^{n}\times {\{{y}_{j}\}}_{j=0}^{n}\times {\{{z}_{k}\}}_{k=0}^{n}$, where

$${x}_{i}=a+i\Delta x,{y}_{j}=c+j\Delta y,{z}_{k}=e+k\Delta z,$$

and Δ*x* = Δ*y* = Δ*z* = (*b − a*)/*n*. Let G* _{l}* be the set of grid points generated by bisecting

**Initialization.**Given an initial Γ, construct a piecewise tri-linear level-set function over the grid*G*. If necessary, apply a reinitialization step to set to be a signed distance function to Γ (see section 4.2 for details). Convert to Φ (see Section 3)._{l}**Evolution.**Resample Φ to obtain a new over the grid*G*. Compute_{l}*L*(Φ) and*H*() in the thin shell (traditionally called a narrow band for curve evolution)*N*= {(*x*)_{i},y_{j}, Z_{k}*G*: |(_{l}*x*)| < }. Update in_{i},y_{j},z_{k}*N*for one time step to get by an ODE time stepping method (see section 4.3).**Re-initialization.**Apply re-initialization step to in the shellto get a new (see section 4.4). Convert to Φ (see Section 3). Go back to step 2 if the termination condition is not satisfied.$$\begin{array}{l}N=\{({x}_{i},{y}_{j},{z}_{k})\in {G}_{l}:\\ \underset{-1\le \mu ,\nu ,\lambda \le 1}{\text{min}}|{\phi}_{i+\mu ,j+\nu ,k+\lambda}|\le \mathscr{H}\}\end{array}$$**Iso-contouring.**Extract 3-sided or 4-sided iso-surface patches (vertices with normals) of Φ =*c*, where*c*is a given iso-value.*G*^{1}approximate these patches by parametric surfaces.

For the problem of molecular surface construction, the grid size *G*_{0} 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.

The aim of using *l* > 0 is to make 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 *G*_{0} for denoising.

Let

$${\beta}^{3}\left(x\right)=\{\begin{array}{cc}\frac{2}{3}-{x}^{2}+\frac{1}{2}{\left|x\right|}^{3},\hfill & 0\le \left|x\right|<1,\hfill \\ \frac{1}{6}{\left(2-\left|x\right|\right)}^{3},\hfill & 1\le \left|x\right|<2,\hfill \\ 0,\hfill & 2\le \left|x\right|.\hfill \end{array}$$

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

$$s\left(x\right)=\sum _{i=1}^{n-1}{c}_{i}{\beta}^{3}\left(x-i\right)$$

(3)

be a cubic spline function. Then the spline interpolation problem is to determine the coefficients $\text{c}={\{{c}_{i}\}}_{i=1}^{n-1}$ such that the following interpolation conditions

$$s\left(k\right)=\sum _{i=1}^{n-1}{c}_{i}{\beta}^{3}\left(k-i\right),k=1,\dots ,n-1$$

(4)

are satisfied for any given function values
${\left\{s\left(k\right)\right\}}_{k=1}^{n-1}$. Such a problem could be easily solved by solving the linear system (4) directly for the unknowns *c*_{1}, …,*c _{n}*

$$s\left(k\right)=\sum _{i=1}^{n-1}{c}_{i}{\beta}^{3}\left(k-i\right),k=0,\dots ,n.$$

Using the initial values

$${c}_{0}^{+}=-\frac{1}{1-{z}_{1}^{2n}}\sum _{k=1}^{n-1}s\left(k\right){z}_{1}^{k}(1-{z}_{1}^{2\left(n-k\right)}),\phantom{\rule{0.2em}{0ex}}{c}_{n}^{-}=0,$$

(5)

the recursive process

$$\begin{array}{l}{c}_{k}^{+}\phantom{\rule{0.4em}{0ex}}=\phantom{\rule{0.4em}{0ex}}s\left(k\right)+{z}_{1}{c}_{k-1}^{+},\phantom{\rule{0.2em}{0ex}}k=1,2,\dots ,n-1,\\ {c}_{k}^{-}\phantom{\rule{0.4em}{0ex}}=\phantom{\rule{0.4em}{0ex}}{z}_{1}({c}_{k+1}^{-}-{c}_{k}^{+}),\phantom{\rule{0.2em}{0ex}}k=n-1,n-2,\dots ,1,\end{array}$$

yields exact B-spline coefficients ${c}_{i}=6{c}_{i}^{-}$ of the interpolation problem, where ${\mathit{\text{z}}}_{1}=\sqrt{3}-2$.

The computation of
${\mathit{\text{c}}}_{0}^{+}$ can be accelerated by neglecting small terms. Given an error control tolerance *ε*, we compute

$$K=\left[\frac{\text{log}\left(\epsilon \right)}{\text{log}(|{z}_{1}|)}\right].$$

If *n* > *K*, we replace
${\mathit{\text{c}}}_{0}^{+}$ in (5) with

$${c}_{0}^{+}=-\frac{1}{1-{z}_{1}^{2n}}\sum _{k=1}^{n-1}s\left(k\right){z}_{1}^{k},$$

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

Since the basis function *β*^{3} (*x*) is locally supported, the *k*th coefficient *c _{k}* of the spline function (3) merely has influence on

$$s\left(i\right)=\{\begin{array}{cc}0,& i\ne k,\\ 1,& i=k,\end{array}$$

define a set of coefficients {*c _{i}*} with none of them zero. However, the coefficients

$${c}_{k}^{-}=-{z}_{1}(1+{z}_{1}^{2}+{z}_{1}^{4}+\cdots )=-\frac{{z}_{1}}{1-{z}_{1}^{2}}=\frac{\sqrt{3}}{6}.$$

Similarly,

$${c}_{k\pm l}^{-}=-{z}_{1}^{l+1}(1+{z}_{1}^{2}+{z}_{1}^{4}+\cdots )={z}_{1}^{l}\phantom{\rule{0.2em}{0ex}}\frac{\sqrt{3}}{6},$$

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

$${c}_{k\pm l}={z}_{1}^{l}\sqrt{3},l=0,1,2,\mathrm{....}$$

Note that
${\mathit{\text{z}}}_{1}=\sqrt{3}-2\approx -0.2679491924$ and c_{k}_{±}* _{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

$${c}_{k\pm l},\phantom{\rule{0.3em}{0ex}}l=0,1,2,3,4.$$

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

$${|{z}_{1}|}^{l}\sqrt{3}<\Delta {x}^{2}\epsilon ,$$

(6)

where Δ*x* is the spacing of the knots.

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

Let * _{l}* be a piecewise tri-linear function over the grid

The method we use in this paper is to convert * _{l}* to

In this section, we describe a local level-set method, using a tri-cubic spline level-set function Φ over the grid *G*_{0} and a piecewise tri-linear level-set function over the grid *G _{l}*, 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.

During the evolution process, we confine the movement of the level-set surface to no more than one grid size Δ*x* of *G _{l}* for each time step. Since ‖‖ = 1 before the movement, we require that function values, and up to the second order partial derivatives of , are accurately computed (within the given error tolerance) in the

$$N=\left\{\mathbf{\text{x}}\in {G}_{l}:\left|\phi \left(\mathbf{\text{x}}\right)\right|\le \mathscr{H}\right\},$$

where

$$\mathscr{H}=\left[\frac{\text{log}{\left(\Delta x\right)}^{2}\epsilon -\text{log}\sqrt{3}}{\text{log}|{z}_{1}|}+2\right]\Delta x+|c|.$$

(7)

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

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 in magnitude.

Let Γ^{0} be a given initial closed surface with interior ^{3}. We define an initial level-set function ^{0}(**x**) in ^{3} satisfying

$$\{\begin{array}{cc}{\phi}^{0}\left(\mathbf{\text{x}}\right)<0\hfill & \text{in}\phantom{\rule{0.3em}{0ex}}\mathcal{D},\hfill \\ {\phi}^{0}\left(\mathbf{\text{x}}\right)=0\hfill & \text{on}\phantom{\rule{0.3em}{0ex}}{\Gamma}^{0},\hfill \\ {\phi}^{0}\left(\mathbf{\text{x}}\right)>0\hfill & \text{elsewhere}.\hfill \end{array}$$

(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 ^{0} needs to be computed differently.

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 . This is done by the following algorithms.

- Set initial value for each grid point to be .
- 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 . 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.

- 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. - Compute distance from the grid points to the (perturbed) surface mesh using Algorithm 4.1.
- First we associate each grid point $\mathbf{\text{q}}\in {\left\{{\mathbf{\text{q}}}_{jk}\right\}}_{jk=0}^{n}:={\left\{{y}_{j}\right\}}_{j=0}^{n}\times {\left\{{z}_{k}\right\}}_{k=0}^{n}$ on the
*yz*-plane with a set of numbers, called**q**-pool. This pool is empty at beginning. Let Δbe the projection of the triangle Δ of the surface mesh onto the_{x}*yz*-plane. For each grid point**q**Δ, compute the intersection of the line segment [_{x}*x*,**q**]^{T}^{3},*x*[*a,b*] with the triangle Δ. Put the*x*-component of the intersection point into**q**-pool. - 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*t*_{1}<*t*_{2}< … <*t*be the sequence in the_{k}**q**-pool. Then if*a*+*j*Δ*x*is in the interval [*t*_{i},t_{i}_{+1}], the sign to the grid point [*a*+*j*Δ*x*,**q**]is set to be (−1)^{T}, where^{i}*t*_{0}=*a, t*_{k}_{+1}=*b*.

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}*,

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

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 *d*(**x**) ≠ **0**, **x** Γ^{0}. Then a simple way is to use *d*(**x**) *− c* or *c − d*(**x**) as the initial ^{0}. For example, in the problem of solvent excluded surface construction, *d*(**x**) is a Gaussian map, *c*=1 (see Section 5).

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

$${N}^{0}=\{\mathbf{\text{x}}\in {G}_{l}:|{\phi}^{0}\left(\mathbf{\text{x}}\right)|<\mathscr{H}\}.$$

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 *C*^{1} Hermite interpolating the boundary data on [_{0},]

$$c\left(x\right)=\{\begin{array}{cc}1\hfill & \text{if}\phantom{\rule{0.3em}{0ex}}\left|x\right|\le {\mathscr{H}}_{0,}\hfill \\ \frac{{\left(\left|x\right|-\mathscr{H}\right)}^{2}\left(2\left|x\right|+\mathscr{H}-3{\mathscr{H}}_{0}\right)}{{\left(\mathscr{H}-{\mathscr{H}}_{0}\right)}^{3}}\hfill & \text{if}\phantom{\rule{0.3em}{0ex}}{\mathscr{H}}_{0}<\left|x\right|\le \hfill \\ 0\hfill & \text{if}\phantom{\rule{0.3em}{0ex}}\left|x\right|>\mathscr{H},\end{array}\mathscr{H},$$

where _{0} is chosen to be 0.5. Our experiments show that linear cut-off function works equally well. We update ^{0} by solving the following equation,

$$\{\begin{array}{c}\frac{\partial \phi}{\partial t}=c\left(\phi \right)\left[L\left(\phi \right)+\mathscr{H}\left(\nabla \phi \right)\right],\hfill \\ \phi \left(\mathbf{\text{x}},0\right)={\phi}^{0}\left(\mathbf{\text{x}}\right),\hfill \end{array}$$

(9)

on *N*^{0} for one time step and get ^{1}(**x**). The time step is chosen such that the interface moves less than one grid size Δ*x*. At each grid point **x*** _{ijk}* in the thin shell

$$\tau =\text{min}\phantom{\rule{0.2em}{0ex}}\left\{\Delta x,\Delta x/\underset{{\mathbf{\text{x}}}_{\mathit{\text{ijk}}}\in {N}^{0}}{\text{max}}\left|{\nu}^{0}\left({\mathbf{\text{x}}}_{\mathit{\text{ijk}}}\right)\right|\right\}.$$

Then update ^{0} by the explicit Euler scheme

$${\phi}^{1}({\mathbf{\text{x}}}_{\mathit{\text{ijk}}})={\phi}^{0}({\mathbf{\text{x}}}_{\mathit{\text{ijk}}})+\tau {\nu}^{0}({\mathbf{\text{x}}}_{\mathit{\text{ijk}}}),\phantom{\rule{0.2em}{0ex}}{\mathbf{\text{x}}}_{\mathit{\text{ijk}}}\in {N}^{0}.$$

Since ^{1} is no longer a signed distance function, a re-initialization step is required (see section 4.4) to get a new ^{1} and a new thin shell *N*^{1}. The process from ^{0} to ^{1} described above is repeated to get a sequence {* ^{m}*}

$$\underset{{\mathbf{\text{x}}}_{\mathit{\text{ijk}}}\in {N}^{m}}{\text{max}}|{\nu}^{m}({\mathbf{\text{x}}}_{\mathit{\text{ijk}}})|<\epsilon \phantom{\rule{0.3em}{0ex}}\text{and}\phantom{\rule{0.3em}{0ex}}m<\text{M}$$

(10)

are satisfied. We choose *ε* = 0.001. *M* is a given upper bound of the iteration number. We choose *M = n*, where *n*^{3} is the number of grid points in *G _{l}*.

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*(* ^{m}*). The higher order term

To simplify the notation, let * ^{m}*(

Re-initialization here is a process of replacing (**x**) by another function (**x**) that is the signed distance function to the zero level-set of and then taking this new function (**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**,*τ*)

$$\{\begin{array}{c}\frac{\partial d}{\partial \tau}+S\left(d\right)\left(\Vert \nabla d\Vert -1\right)=0,\hfill \\ d\left(\mathbf{\text{x}},0\right)=\phi \left(\mathbf{\text{x}}\right),\hfill \end{array}$$

(11)

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

$$\{\begin{array}{c}\frac{\partial d}{\partial \tau}+S\left(d\right)(\nabla {d}^{T}\nabla d-1)=0,\hfill \\ d\left(\mathbf{\text{x}},0\right)=\phi \left(\mathbf{\text{x}}\right).\hfill \end{array}$$

(12)

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

$$s=\frac{d}{\sqrt{{d}^{2}+{\left|Dd\right|}^{2}\Delta {x}^{2}}},$$

(13)

where *Dd* is a discrete approximation of *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.

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)

$$\begin{array}{l}{d}_{\mathit{\text{ijk}}}^{n+1}\phantom{\rule{0.4em}{0ex}}=\phantom{\rule{0.4em}{0ex}}{d}_{\mathit{\text{ijk}}}^{n}-\Delta \tau {s}_{\mathit{\text{ijk}}}^{+}(\text{max}[{\left({a}^{+}\right)}^{2},{\left({b}^{-}\right)}^{2}]\hfill \\ \phantom{\rule{1em}{0ex}}+\text{max}[{\left({c}^{+}\right)}^{2},{\left({d}^{-}\right)}^{2}]+\text{max}[{\left({e}^{+}\right)}^{2},{\left({f}^{-}\right)}^{2}]-1)\hfill \\ \phantom{\rule{1em}{0ex}}-\Delta \tau {s}_{\mathit{\text{ijk}}}^{-}(\text{max}[{\left({a}^{-}\right)}^{2},{\left({b}^{+}\right)}^{2}]+\text{max}[{\left({c}^{-}\right)}^{2},{\left({d}^{+}\right)}^{2}]\hfill \\ \phantom{\rule{1em}{0ex}}+\text{max}[{\left({e}^{-}\right)}^{2},{\left({f}^{+}\right)}^{2}]-1),\hfill \end{array}$$

(14)

where *S _{ijk}* is the approximation to

$$\begin{array}{ccc}a={D}_{x}^{-}{d}_{\mathit{\text{ijk}}}^{n},& c={D}_{y}^{-}{d}_{\mathit{\text{ijk}}}^{n},& e={D}_{z}^{-}{d}_{\mathit{\text{ijk}}}^{n},\\ b={D}_{x}^{+}{d}_{\mathit{\text{ijk}}}^{n},& d={D}_{y}^{+}{d}_{\mathit{\text{ijk}}}^{n},& f={D}_{z}^{+}{d}_{\mathit{\text{ijk}}}^{n},\end{array}$$

respectively, where
${D}_{x}^{\pm}{d}_{\mathit{\text{ijk}}},{D}_{y}^{\pm}{d}_{\mathit{\text{ijk}}}$ and
${D}_{z}^{\pm}{d}_{\mathit{\text{ijk}}}$ denote the one-sided divided differences at the grid point **x*** _{ijk}* in the

$$\frac{\Delta \tau}{\Delta x}|{s}_{\mathit{\text{ijk}}}|\le \frac{1}{3}.$$

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

During the process of re-initialization, function values around the level-set change, and so the thin shell defined by {**x** *G _{l}* : |(

In this section, we use our higher order level-set method with *C*^{2} 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
${\sum}_{i=1}^{m}{e}^{-d({\Vert \mathbf{\text{x}}-{\mathbf{\text{x}}}_{i}\Vert}^{2}-{r}_{i}^{2})/{r}^{2}}$ 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
${\{{\mathbf{\text{x}}}_{i}\}}_{i=1}^{m}$ and radii
${\{{r}_{i}\}}_{i=1}^{m}$ (see Figure 3(a)) is retrievable from the Protein Data Bank (PDB). To construct the SES for *M*, we minimize the energy functional

(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(\Gamma )={\int}_{\Gamma}g{\left(\mathbf{\text{x}}\right)}^{2}d\mathbf{\text{x}}+\epsilon {\int}_{\Gamma}d\mathbf{\text{x}},$$

(15)

where

$$g\left(\mathbf{\text{x}}\right)=1-\sum _{i=1}^{m}{e}^{-C({\Vert \mathbf{\text{x}}-{\mathbf{\text{x}}}_{i}\Vert}^{2}-{\stackrel{\sim}{r}}_{i}^{2})},\phantom{\rule{0.5em}{0ex}}{\stackrel{\sim}{r}}_{i}={r}_{i}+{r}_{b},$$

with *r _{b}* the probe radius, which is usually 1.4 Å. The constant

$$\begin{array}{l}\frac{\partial \phi}{\partial t}=({g}^{2}+\epsilon )\text{div}\phantom{\rule{0.2em}{0ex}}\left(\frac{\nabla \phi}{\Vert \nabla \phi \Vert}\right)\phantom{\rule{0.2em}{0ex}}\Vert \nabla \phi \Vert +4g{\left(\nabla g\right)}^{T}\nabla \phi \hfill \\ \phantom{\rule{1.6em}{0ex}}=L\left(\phi \right)+H\left(\nabla \phi \right),\hfill \end{array}$$

(16)

where

$$\begin{array}{l}\phantom{\rule{1em}{0ex}}L\left(\phi \right)=({g}^{2}+\epsilon )\text{div}\phantom{\rule{0.2em}{0ex}}\left(\frac{\nabla \phi}{\Vert \nabla \phi \Vert}\right)\phantom{\rule{0.2em}{0ex}}\Vert \nabla \phi \Vert ,\hfill \\ H\left(\nabla \phi \right)=4g{\left(\nabla g\right)}^{T}\nabla \phi .\hfill \end{array}$$

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

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

In step 2, (**x**) may change rapidly, we may slow down the changes of by dividing it with a constant before taking the re-initialization step. Since the gradient of
${e}^{-C[{\Vert \mathbf{\text{x}}-{\mathbf{\text{x}}}_{i}\Vert}^{2}-{\stackrel{\sim}{r}}_{i}^{2}]}$ is
$-2C{e}^{-C[{\Vert \mathbf{\text{x}}-{\mathbf{\text{x}}}_{i}\Vert}^{2}-{\stackrel{\sim}{r}}_{i}^{2}]}{[x-{x}_{i},y-{y}_{i},z-{z}_{i}]}^{T}$, the length of the gradient is approximately
$2C{\stackrel{\sim}{r}}_{i}$ at the sphere surface
$\{\mathbf{\text{x}}:\Vert \mathbf{\text{x}}-{\mathbf{\text{x}}}_{i}\Vert ={\stackrel{\sim}{r}}_{i}\}$. Hence we normalize the function by dividing it with
$2C\ast {\text{max}}_{i}{\stackrel{\sim}{r}}_{i}$.

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

$${e}^{-C[{\Vert \mathbf{\text{x}}-{\mathbf{\text{x}}}_{i}\Vert}^{2}-{\stackrel{\sim}{r}}_{i}^{2}]}={e}^{-C[{({\stackrel{\sim}{r}}_{i}+\epsilon )}^{2}-{\stackrel{\sim}{r}}_{i}^{2}]}=\frac{1}{k}.$$

From this we cobtain

$$C\approx \frac{\text{ln}\phantom{\rule{0.1em}{0ex}}k}{2{\stackrel{\sim}{r}}_{i}\epsilon}.$$

For instance, if *k* is 4, we choose *ε* = 0.01, and *C* is thus chosen to be
$\frac{50\phantom{\rule{0.1em}{0ex}}\text{ln}\phantom{\rule{0.1em}{0ex}}4}{{r}_{b}+{\text{min}}_{i}{r}_{i}}$.

Since each term
${e}^{-C[{\Vert \mathbf{\text{x}}-{\mathbf{\text{x}}}_{i}\Vert}^{2}-{\stackrel{\sim}{r}}_{i}^{2}]}$ in *g*(**x**) decreases rapidly as ‖**x** − **x*** _{i}*‖ increases, we evaluate
${e}^{-C[{\Vert \mathbf{\text{x}}-{\mathbf{\text{x}}}_{i}\Vert}^{2}-{\stackrel{\sim}{r}}_{i}^{2}]}$ locally around

$${e}^{-C[{\Vert \mathbf{\text{x}}-{\mathbf{\text{x}}}_{i}\Vert}^{2}-{\stackrel{\sim}{r}}_{i}^{2}]}>\epsilon ,$$

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

$${\Vert \mathbf{\text{x}}-{\mathbf{\text{x}}}_{i}\Vert}^{2}<{\stackrel{\sim}{r}}_{i}^{2}-\frac{\text{ln}\phantom{\rule{0.1em}{0ex}}\epsilon}{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.

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.

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 = −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** ^{3} : *g*(**x**) = 1}, where

$$g\left(\mathbf{\text{x}}\right)=\sum _{i=1}^{n}{e}^{-\frac{d}{{r}_{i}^{2}}[{\Vert \mathbf{\text{x}}-{\mathbf{\text{x}}}_{i}\Vert}^{2}-{r}_{i}^{2}]}.$$

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)).

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 *C*^{2} 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 *C*^{1} and higher-order smoothness is both necessary and/or desirable.

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)

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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |