The main aim of this study was to develop and implement an algorithm for the rapid, accurate and automated identification of paths leading from buried protein clefts, pockets and cavities in dynamic and static protein structures to the outside solvent.
The algorithm to perform a skeleton search was based on a reciprocal distance function grid that was developed and implemented for the CAVER program. The program identifies and visualizes routes from the interior of the protein to the bulk solvent. CAVER was primarily developed for proteins, but the algorithm is sufficiently robust to allow the analysis of any molecular system, including nucleic acids or inorganic material. Calculations can be performed using discrete structures from crystallographic analysis and NMR experiments as well as with trajectories from molecular dynamics simulations. The fully functional program is available as a stand-alone version and as plug-in for the molecular modeling program PyMol. Additionally, selected functions are accessible in an online version.
The algorithm developed automatically finds the path from a starting point located within the interior of a protein. The algorithm is sufficiently rapid and robust to enable routine analysis of molecular dynamics trajectories containing thousands of snapshots. The algorithm is based on reciprocal metrics and provides an easy method to find a centerline, i.e. the spine, of complicated objects such as a protein tunnel. It can also be applied to many other molecules. CAVER is freely available from the web site http://loschmidt.chemi.muni.cz/caver/ webcite.
The shape of a protein is complicated by its many clefts, pockets, protrusions, channels and cavities. Protein concavities offer a unique microenvironment for biological functions, such as ligand binding or enzymatic catalysis. Protein shape is of great interest to medicinal chemists working in the drug discovery industry and looking for inhibitors, enzymologists interested in identifying substrate molecules based on the well known "lock and key" mechanism and protein chemists studying protein-protein or protein-DNA interactions. The identification of protein pockets and cavities has been the focus of several studies [1-4] and various algorithms have been developed for the calculation of protein volume and surface area. A large number of enzymes possess buried active sites that are connected to the external solvent environment by access routes (tunnels or channels). A catalytic step must always be preceded by the formation of an enzyme-substrate complex, which may require passage of the substrate through these routes. The size and shape of the access routes may become an important determinant of enzyme substrate specificity . Changes in the diameter of the access tunnels during the dynamic movement of a protein play an important biological role, such as that described for acetylcholinesterase . Two narrow active site gorges are positioned deep inside the protein core and movement of the residues making up the gorge walls is necessary to allow ligands access to the active site. A method based on molecular surface was used for the calculation of the gorge diameter in acetylcholinesterase. The diameter was defined as the maximum probe size that produces a continuous molecular surface between an active site and a solvent. Calculation of one diameter in this approach requires the generation of several molecular surfaces using a series of probes of increasing size . A more effective method is implemented in the CAST program, which utilizes the alpha shape theory. CAST computation of pockets and their openings does not require direct human interaction. The required inputs are atomic coordinates, van der Waals radii, and the radius of the probe sphere . The program VOIDOO, a component of O package utilizes a grid-based algorithm for detection, delineation, and measurement of protein cavities and solvent accessible pockets. The VOIDOO algorithm suffers from crude grid spacing and the "can-of-worms" phenomenon . The central problem in the analysis of tunnels in protein structures is the identification of the centerline, i.e. spine, of a 3D object. Algorithms dealing with centerlines have been applied to medical procedures, for example in virtual colonoscopy and bronchoscopy [8-11].
The aim of this study was to develop a rapid and accurate algorithm for the identification of routes from buried active sites to the external solvent in static protein structures. We aimed to produce an algorithm that could also be applied to molecular dynamic trajectories. Further, the algorithm was intended to allow changes in the radius of a channel gorge with time to be monitored and the most probable access routes to be identified. Several other requirements were taken into consideration during development of the algorithm and its implementation: (i) speed, thus enabling rapid analysis of an entire trajectory from a molecular dynamic simulation, i.e. thousands of snapshots, in a few hours; (ii) easy identification of a starting point for the calculation; (iii) that the algorithm functions independently of the probe radius; (iv) storage of paths in PDB format; and (v) intuitive visualization.
The most accessible path from the protein cavity to the bulk solvent has to be found by systematic exploration of the protein interior, in order to calculate the access route gorge radius (Fig. 1). In our model, a protein consists of hard sphere atoms with appropriate van der Waals radii. The protein body is modeled on a discrete three-dimensional grid space and all grid nodes are clustered into two classes: nodes located in the interior of the protein body (inside atomic vdW radii) and nodes located outside the protein body. Outer nodes can lie in the cavity, access tunnels or in the external environment of the protein, e.g. a bulk solvent. The convex approximation of the protein, termed the 'convex hull', is used to distinguish nodes that lie either in the interior or exterior of the protein (Fig. 1). Nodes that are located outside of the convex hull, are eliminated and not used in further calculations.
Figure 1. Sketch of the method implemented in CAVER. The black bold circle represents the starting point. The protein is visualized by gray circles with van der Walls atom radii mapped on a discrete grid (black dots). The solid line represents the boundary between the protein (convex hull) interior and its surroundings. Empty circles represent the maximally inscribed balls on the probable route (dashed line).
Attention is paid to nodes that lie on a boundary of the modeled convex hull. These nodes are potential end-stops of the grid search algorithm because each boundary node can be treated as a putative outfall of the channel. The mathematical object, which is called a vertex-weighted graph, is constructed and the algorithm applied to identify the shortest low-cost path. Each possible path from the active site to the exterior is evaluated as a positive value. This value represents the relative cost to navigate each path in what could be described as a "highway-toll". Long and complicated paths are more "expensive", while the short direct paths are "cheaper" (Fig. 2). More formally, the cost function C(P) is defined (Eq. 1) for the given path P as the sum of node-price-function values calculated for all nodes forming the path P. Let N(P) be the set of all points form the path P, then the C(P) is expressed as
Figure 2. Evaluation of grid nodes by cost function. A grid point evaluation using the cost function (Eq. 3). The line represents the optimal centerline.
This cost function depends on the number of nodes in each sum and, as such, this function is not suitable for the purposes of comparison. A normalized cost function is defined (Eq. 2) to avoid cost function dependence on a summand number:
where N is the number of the summand. Next, the single node cost function c(x) must fulfill two requirements. First, it must provide a positive value for each node and a low-value for nodes that are surrounded by empty space. Second, it must identify preferred nodes that are surrounded by sufficient empty space to allow a hypothetical substrate to pass through a channel without risk of collision. These low-weighted nodes are preferentially selected by the search-algorithm. In our case, the cost function c(x) for a single node was chosen (Eq. 3)
where function rmax (x) is equal to the maximal radius of a hypothetical ball that can be inserted into node x just touching the protein surface. The small constant ε is here only for technical purposes to get rid of a singularity of the function in points where rmax (x) equals zero.
The graph-searching algorithm then establishes the lowest cost path from the active site to the external environment. The calculated path can be visualized (Fig. 3 and Fig. 4) using the rmax radii for each node of the path. The smallest radius represents the channel gorge and as such the point coordinates can be determined together with the gorge radius rgorge.
Figure 3. Access path visualized by pyMol. Visualization of the access route using the PyMol plug-in. Wires represent the protein surface, balls are nodes and the surface represents the export route.
Figure 4. Path profile convergence. Convergence of the path profile found on the grid with increasing precision, i.e. decreasing node distance d. Calculation of one path takes approximately 11 sec in case of d = 0.7°A (Athlon 2600+, 2GB RAM, NetBSD 1.6.1) but increases ten-fold in a substrate where d = 0.3°A.
The method was implemented in the CAVER program (1 and web page http://loschmidt.chemi.muni.cz/caver/ webcite. The program uses the publicly-licensed software qhull  that quickly (in O(N log N) time) computes a convex hull for a given set of N points in three dimensions. The result of qhull was used to eliminate nodes located outside the convex hull. All points were investigated regardless of whether they lay inside or outside the convex hull. This can be achieved by traversing all facets of the convex hull and testing whether a point lies in the same halfspace as another convex hull point. This process takes O(NK) time, where N is the number of points inside the convex hull, and K is the number of facets forming the convex hull.
Additional File 1. Source codes of CAVER. Additional file (gzipped tar archive caver_unix_v0.99.4.tar.gz) contains source codes of CAVER ver. 0.99.4 for UNIX platforms and pyMol plug-in. Instructions for the installation and updates are provided at web page http://loschmidt.chemi.muni.cz/caver/ webcite.
Format: GZ Size: 1.1MB Download file
Based on graph theory, several methods for the shortest path problem have been described. The most widely used algorithms are Dijkstra's , Bellman-Ford's [14,15], A* search algorithm, and the Floyd-Warshall algorithm . In our case, a positively vertex-weighted graph is plotted on a three dimensional grid, where the source vertex is known while the destination vertex is not. We used a modified form of the Dijkstra's algorithm.
The Dijkstra's algorithm effectively solves the problem of the shortest path from a single source vertex to the destination one. It was originally designed for edge-weighted graphs but its vertex-weighted variation is easy to implement. The algorithm can be used even if the destination vertex is unknown. In the main loop of the algorithm, the shortest path to the closest available vertex (measured from the source vertex) is determined. Then, estimates of the shortest path for all the adjacent vertices are updated. This means searching can be terminated if the nearest available vertex is the boundary vertex indicating that the shortest path has been identified. To speed up the algorithm, a modification related to the cost function was implemented. The single node cost function is evaluated as part of the main loop of the Dijkstra's algorithm. The cost function is evaluated only at nodes where it is required.
The identified path can be easily visualized since the program writes a PDB file containing the coordinates of the path nodes accompanied by the maximum probe radius that touches the vdW protein surface (Fig. 3). A user-friendly GUI was implemented in the graphic program PyMol .
An algorithm to perform a skeleton search on a defined grid was developed and implemented in the CAVER program as described in the Methods section. The algorithm developed automatically finds the easiest path from the starting point, typically located inside the molecule, to the exterior of the molecule. The identified path resembles a tunnel that connects protein clefts, pockets or cavities with the surrounding bulk solvent. The tunnel characteristics, e.g. length, mean radius and gorge radius are determined and can be further analyzed. In molecular dynamics trajectories it is possible to analyze time fluctuations of tunnel characteristics and construct a dynamic picture of tunnel behavior.
The tunnel gorge radius rgorge is one of the most important tunnel characteristics because the tunnel gorge can form a bottleneck for substrate access or product release to and from the active site of a protein. The radius rgorge as estimated by the algorithm is always underestimated in finite grids. The maximal error "max" of an rgorge estimation is expressed by the equation (Eq. 4):
where d is equal to the length of the grid cell edge. The probability of εmax realization is equal to zero and this error is overestimated, therefore the mean error should be defined. The mean error of rgorge determination is equal to (Eq. 5):
and its variance and deviation equal to (Eq. 6)
The rgorge estimation should be corrected by adding 0.48d to the rgorge. The corrected rgorge estimate has a mean error value με = 0 and the variance of error is var(ε) = 0.019d2. In the case of a globular protein (50 × 50 × 50 Å3) the εmax (rgorge) costs 0.43 Å for a grid with d = 0.5 Å, however, the mean ε(rgorge) equals 0.24 Å, and σε = 0.07 Å. The results of the tests (depicted in Fig. 4) focus on the convergence of the identified paths with d decreasing from 1.0 to 0.3 Å.
Performance of the method is given as the tunnel volume i.e. number of vertexes searched rather than the number of atoms. In the case of haloalkane dehalogenases (the active site volume ~200Å3) the typical calculation of one tunnel takes about 10–12 sec. but in the case of cytochrome P450 2C9 or 3A4 (which has larger active sites ~500–600Å3) the calculation takes about 20–25 sec. In case of very large cavities (e.g. RNA) calculation may take several minutes at low resolution (d = 1–2Å). The program performance was tested on Pentium IV 3.2 GHz machine with 2 GB RAM running on Windows XP Professional operating system.
Haloalkane dehalogenases (EC 18.104.22.168) are microbial enzymes that cleave a carbon-halogen bond in a broad range of halogenated compounds . The molecular structures of three different haloalkane dehalogenases are known: DhlA from Xanthobacter autotrophicus GJ10 [19-26], DhaA from Rhodococcus sp.  and LinB from Sphingomonas paucimobilis UT26 [28-31]. The overall shape of haloalkane dehalogenases is globular, with the active site buried between the main domain of an α/β hydrolase fold and a cap domain with a uteroglobin fold. There are three access routes connecting the protein surface with the active site, denoted the main tunnel, the upper tunnel and a slot (Fig. 5). The three proteins differ in the number of routes that provide access to their active site: LinB has the most available active site, accessible through three tunnels. The active site in DhaA is accessible through the upper tunnel and the slot, while DhlA is believed to have a single accessible route via the main tunnel . Here we used CAVER to conduct a thorough analysis of the access paths using all available X-ray structures and pre-calculated molecular dynamics trajectories.
Figure 5. Scheme of haloalkane dehalogenases' tunnels. Schematic representation of access paths for DhlA (A), DhaA (B) and LinB (C) identified by protein crystallography [19–31] and molecular dynamic simulations . The slot in DhlA was described in this study.
Analysis of X-ray structures
Analysis of fifteen available crystal structures of the DhlA enzyme (1B6G,
Figure 6. Tunnels found by CAVER. Accessible paths identified by CAVER in DhlA (A), DhaA (B) and LinB (C). The most accessible tunnel in every structure (main in DhlA, upper in DhaA and lower in LinB) is colored in red. The upper tunnel of LinB is in blue. Slots are highlighted in green for each structure.
Table 1. Characteristics of tunnels. Characteristics of the four most accessible paths connecting the active sites of DhlA, DhaA and LinB with the bulk solvent.
The three crystal structures available for the DhaA enzyme (1BN6, 1BN7 and 1CQW) were analyzed for the presence of access routes. CAVER detected one clearly preferred route, which corresponded to the upper tunnel. The gorge of this tunnel is made up of W141, F144, A145, F159, V172 and C176 (numbered according to dhaA sequence). The additional paths were located in the slot (Fig. 6B). The slots of the three structures studied showed a slightly different spatial position and variable size with respect to the mean gorge radius, mainly due to repositioning of the side-chain R133. The cost function of these routes was almost twice that of the upper tunnel (Table 1), but still comparable with the main tunnel and the slot of DhlA, which are known to be involved in substrate binding and product release. Exchange of solvent molecules between the active site and the slot was observed during a 1ns molecular dynamic simulation .
Analysis of eleven available LinB crystal structures (
Analysis of molecular dynamics trajectories
A molecular dynamics simulation of the DhlA enzyme was analyzed using CAVER to determine the easiest routes from the active site. About 400 snapshots taken at 5 ps intervals were analyzed, the tunnels identified and their corresponding gorgeswere further analyzed by visual inspection (Fig. 7A). CAVER identified two clusters of gorges that correspond to two different paths. The most populated tunnel gorges (78%) are located in the main tunnel and the other remaining gorge clusters are located in the slot.
Figure 7. Tunnels' gorges as found by CAVER in MD trajectories. Analysis of snapshots taken from a molecular dynamics simulation of DhlA (A), DhaA (B) and LinB (C). Two clusters of tunnel gorges were identified in DhlA simulation (denoted by the numbers 1 and 2). Three clusters were identified in the DhaA simulation (denoted by the numbers 1–3) and one cluster was identified in the LinB simulation. Tunnel gorges are represented by small balls. Ball color correlates with gorge radius, balls representing narrow gorges are red and wide gorges are blue. The mesh represents the protein surface.
As in the previous case, the molecular dynamics simulation of the DhaA enzyme was analyzed using the CAVER program. Two main clusters were identified by CAVER, one having two subclusters (Fig. 7B) of gorges resulting in three different paths,. The most populated access gorges (64%) were located in the upper tunnel and the two remaining gorge clusters were positioned in the slot. The two subclusters in the slot (Fig. 7B) have populations equal to 26% (cluster 2) and 10% (cluster 3), respectively. Analysis of the LinB molecular dynamics trajectory identifies the main tunnel as the most accessible route to the active site (Fig. 7C).
A new algorithm for the identification of tunnels in large molecular systems was developed and implemented in the CAVER program, which is available within the public domain. The algorithm automatically explores a grid, which is constructed over the molecule and stripped to its convex hull. Nodes are evaluated using a cost function, which determines the amount of free space around the node. The grid search algorithm is used to find the lowest-cost centerline path between a given starting point and the surface of the molecule. The user needs only to provide the molecular geometry, atomic van der Waals radii and the designated starting point, to enable the analysis of any molecular system be it protein, nucleic acid or inorganic material. The algorithm is sufficiently rapid and robust for the routine analysis of molecular dynamics trajectories that contain thousands of snapshots. The program is also available as a plug-in for PyMol and, additionally a Web-based version of the program offers analysis of static protein structures online.
Availability and requirements
Project name: CAVER;
Operating systems: UNIX/Linux, Windows;
Programming language: C++;
Licence: CAVER licence;
Any restrictions to use by non-academics: licence needed.
MP developed the algorithm for search of access paths, wrote and tested the software, prepared initial draft of the manuscript and prepared web pages; MO developed concept of the software, conducted performance tests, wrote parts of the manuscript and prepared graphic material for the web; PB contributed ideas on algorithm and conducted statistical analyses; PK conducted performance tests; JK financially supported MP; JD contributed fundamental biochemical concept and interpreted data, wrote parts of the manuscript, contributed ideas on web pages and coordinated project.
The authors express their thanks to Petr Kulhánek (NCBR, Brno, Czech Republic) for his valuable suggestions. This work was supported by the Czech Ministry of Education and the Grant Agency of the Czech Republic (MSM6198959216, MSM0021622413, LC512 and 305/06/P139). J.D. is funded by an EMBO/HMMI Young Investigator grant.
Chaloupkova R, Sykorova J, Prokop Z, Jesenska A, Monincovaa M, Pavlova M, Tsuda M, Nagata Y, Damborsky J: Modification of activity and specificity of haloalkane dehalogenase from Sphingomonas paucimobilis UT26 by engineering of its entrance tunnel.
Transactions on Visualization and Computer Graphics 2001, 7:195-206. Publisher Full Text
Transactions on Medical Imaging 2002, 21:1450-1460. Publisher Full Text
Communications of the ACM 2005, 48:37-41. Publisher Full Text
ACM Transactions on Mathematical Software 1996, 22:469-483. Publisher Full Text
Communications of the ACM 1962, 5:345. Publisher Full Text
Environmental Toxicology and Chemistry 2001, 20:2681-2689. Publisher Full Text
Verschueren KHG, Kingma J, Rozeboom HJ, Kalk KH, Janssen DB, Dijkstra BW: Crystallographic and fluorescence studies of the interaction of haloalkane dehalogenase with halide-ions – studies with halide compounds reveal a halide binding-site in the active-site.
Schanstra JP, Ridder IS, Heimeriks GJ, Rink R, Poelarends GJ, Kalk KH, Dijkstra BW, Janssen DB: Kinetic characterization and X-ray structure of a mutant of haloalkane dehalogenase with higher catalytic activity and modified substrate range.
Oakley AJ, Prokop Z, Bohac M, Kmunicek J, Jedlicka T, Monincova M, Kuta-Smatanova I, Nagata Y, Damborsky J, Wilce MCJ: Exploring the structure and activity of haloalkane dehalogenase from Sphingomonas paucimobilis UT26: Evidence for product- and water-mediated inhibition.
Oakley AJ, Klvana M, Otyepka M, Nagata Y, Wilce MCJ, Damborsky J: Crystal structure of haloalkane dehalogenase LinB from Sphingomonas paucimobilis UT26 at 0.95 Angstrom resolution: Dynamics of catalytic residues.
Streltsov VA, Prokop Z, Damborsky J, Nagata Y, Oakley A, Wilce MCJ: Haloalkane dehalogenase LinB from Sphingomonas paucimobilis UT26: X-ray crystallographic studies of dehalogenation of brominated substrates.