Internal:GPCD

Introduction

Protein Cavities

Cavities are structural features of proteins and they are often important for drug discovery and design. There are different kinds of cavities, including pockets, voids, tunnels, and channels. The types of cavities are listed in more detail in figure 1.

Figure 1: The basic types of protein cavities [1]

1. Void: A cavity that looks like internal "hole" in a protein with no connection to the surface.
2. Cleft: A pocket that looks like a "dent" in the surface of a protein. Can also be found as part of a tunnel or channel.
3. Invagination: A pocket that looks like a "cove". It has a small opening on the surface of a protein and leads to a larger interior. Can also be found as part of a tunnel or after a cleft.
4. Tunnel: A pocket that has an opening on the surface of a protein and leads to somewhere in the interior. Has no opening on the other side.
5. Channel: A cavity that has an opening on one side of the surface of a protein and leads to another opening in another location on the surface of the protein.

Of course, the types of cavities listed in figure 1 are only the basic ones. Combinations of these are also possible as shown in figure 2.

Figure 2: Combination types of protein cavities [1]

Protein Cavity Detection Methods

There are several methods for detecting protein cavities. There are sphere based methods, grid based methods, grid-and-sphere based methods, surface based methods, grid-and-surface based methods, time-varying methods, GPU based methods, consensus based methods, and finally tessellation based methods. These algorithms are shown in more detail in figure 3. Tessellation based methods mainly rely on Voronoi tessellation. Some methods rely on weighted Voronoi tessellation because of its ability to more accurately represent the volume of atoms. The Voronoi tessellation in these methods is computed from a set of points representing the atomic centers of the underlying protein.

Figure 3: A hierarchy of protein cavity detection algorithms [1]

Literature Review

MOLE

Method

The MOLE algorithm is specifically designed to identify tunnels and channels. Finding superficial cavities is a secondary concern and does so as part of the surface approximation step. MOLE utilizes the surface obtained by the refinement of the convex hull to find superficial cavities such as clefts and invaginations. In essence, the approximation of a more accurate surface and the detection of these superficial cavities are closely related. As tetrahedra are removed from the convex hull to get a more accurate surface approximation, a parameter called interior threshold controls how many superficial cavities are found. For example, if the number of tetrahedra removed from a specific part of the molecule exceeds this threshold, the algorithm marks these tetrahedra as part of a superficial cavity.

Figure 4: Shows in 2D the Voronoi diagram (thin lines) and the convex hull (dotted line) of molecule represented by atoms (circles). Each Voronoi edge is evaluated by the cost function (numbers). The thick line represents the optimal path found by the Dijkstra's algorithm from the given starting point (black circle), i.e., centerline of found tunnel. [1]

For the more complicated types of cavities, tunnels and channels, MOLE uses the Voronoi diagram to identify them. Start and end points are defined as the centers of the deepest tetrahedra in each cavity. The depth of the tetrahedron is defined as the number of Voronoi edges from the closest boundary tetrahedron (on the surface). Tunnels/channels are computed as the shortest paths between all pairs of start and end points in the same cavity diagram.

Figure 5: Shows the weight function used to estimate the weight of each Voronoi edge for Dijkstra's algorihtm. l(e) is the length of the edge, d(e) is the distance of the edge to the closest atom Van der Waals sphere and ϵ is a small number to avoid division by zero [1]

Complexity

• Worst-case: O(N^2 log N), where N is the number of atoms in the molecule.
• Average: O(M log M), where M is the number of vertices in the Voronoi diagram (in the worst case, M = N^2, but in most cases M = O(N)). As a result of the use of the incremental algorithm and Hilbert curve ordering, the complexity of calculating the Delaunay triangulation of most molecular structures is O(N log N). Finally, the complexity of all the remaining steps of the algorithm is at most O(M log M).

MolAxis

Method

The MolAxis algorithm relies on Voronoi and Dijkstra's algorithm to find channels, similarly to MOLE. The algorithm consists of the following high-level steps:

1. Approximate the molecule using a collection of fixed-size spheres. The number of balls needed depends on the required quality of the approximation. In the extreme case, where the user is interested in a high resolution, each atom is approximated using a single sphere.
2. Compute the Voronoi diagram of the centers of the fixed-sized spheres. The diagram is intersected with a user-defined bounding sphere. This defines the area of interest.
3. Discard the parts of the Voronoi diagram that intersect the molecule.
4. Define a graph, using the Voronoi vertices and edges. A weight is given to each edge, which depends on its length and on its distance from the surface of the molecule.
5. Define a source point (root vertex). The user can either supply a source point or allow MolAxis to automatically choose the source point to be the center of the largest cavity in the molecule (basically detect channels inside-out).
6. Apply a shortest path optimization algorithm (Dijktra's algorithm) on the graph, computing a tree called the "corridor tree".
7. The desired channels are paths in the corridor tree. Duplicate channels are avoided by using the split distance parameter. The split distance bounds the distance to the least common ancestor of the channel and all previously reported channels.

Figure 6: Shows an example of fixed-size circles (light-blue) approximating varying-size circles and the pathway diagram of their centers. The discarded portion of the Voronoi diagram of their centers is depicted using dotted lines [2]. This figure also helps illustrate the usefulness of the Voronoi tessellation in computing pathways.

Protein Energetics

Metrics

1. Hydrophobicity: Physical property of a molecule that is seemingly repelled from a mass of water. Hydrophilic molecules are attracted to water. Hydrophobicity plays a major role in protein folding. Hydrophobic molecules also tend to be nonpolar and prefer other neutral molecules and nonpolar solvents. Hydrophobicity is the most widely used metric.
2. Hydrogen Bonds: Weak bonds between two molecules resulting from an electrostatic attraction between a proton in one molecule and an electronegative atom in the other. We often see these bonds between oxygen and hydrogen.
3. Van der Waals Interactions: Also known as Van der Waals forces, they are forces driven by induced electrical interactions between two or more atoms or molecules that are very close to each other. Van der Waals interaction is the weakest of all intermolecular attractions between molecules. These forces are distance-dependent.
4. Salt Bridges: Non-covalent (no sharing of eletrons; electromagnetic) interaction between two ionized sites. It has two components: a hydrogen bond and an electrostatic interaction.
5. Disulfide Bridges: Strong bonds that can form between two cysteines (an amino acid). The strength of disulfide bridges helps stabilize a protein. Disulfide bridges are especially common in proteins that are secreted from cells.
6. Charge: Plays a role in strength of bonds. A higher charge means that an atom will have stronger forces in an electromagnetic field.
7. Polarity: Separation of electric charge leading to a molecule or its chemical groups having an electric dipole moment, with a negatively charged end and a positively charged end. Indicates the existence of polar bonds, which are covalent (sharing of electrons) and are quite strong. About bond polarity.

Candidate Method: Protein Cavity Identification using Triangulation, unsigned Distance, and Clustering (PCITDC)

The PCITDC method employs Delaunay triangulation to generate a coarse and fine surface of a protein (the fine surface is enclosed within the coarse surface). Then, the unsigned distance of the fine surface from the coarse surface is computed. This is done by iterating through every vertex (atom) in the fine surface and computing its distance from the closest point on the coarse surface. Finally, three-level clustering is employed to cluster atoms based on their distance values.

The algorithm contains the following high-level steps:

(Spiros/Αngelos) Compute Inaccurate and Accurate surfaces

1. Solvent Excluded Surface (SES) image
   1. Create a space filling image with spheres for each atom where their radius is: Van der Waals radius + probe radius [DONE]
      1. Note: Shared memory parallelism can be employed [DONE]
   2. Perform Euclidean Distance Transform on the space filling image [DONE]
      1. Note: Shared memory parallelism can be employed [DONE]
   3. Extract Solvent excluded surface from the EDT result by setting: [DONE]
      1. All pixel values < probe-radius = 0
      2. All pixel values >= probe-radius = 1
      3. Note: Shared memory parallelism can be employed [DONE]
2. Solvent Excluded Surface (SES): [DONE]
   1. Computed using:
      1. PODM3D
   2. Advantages:
      1. Can control the number of elements on the surface
      2. the edge lengths bounds are relatively small which helps in the clustering.
   3. Disadvantages:
      1. Doesn't have a direct relation to the molecule
         1. But can be fixed using a KD-tree [DONE]
      2. Surface may need to be smoothed
   4. Note: Shared memory parallelism can be employed [DONE]

(Angelos/Spiros) Compute ambient occlusion of the SES

1. Computed using:
   1. Using Embree library (of Intel) and Eigen [DONE]
      1. Note: Shared memory parallelism can be used [DONE]
      2. Note: GPU parallelism can be used [TODO]

(Spiros) Surface Statistics

1. Compute min/mean/median/max edge lengths of the surface triangles [DONE]
   1. Note: Shared memory parallelism can be used [DONE]

(Angelos) Compute atom clusters based on ambient occlusion values:

1. A modified version of DBSCAN that finds neighbors based on not only their density but also their ambient occlusion values.
   1. Variations:
      1. One-level clustering [DONE]
      2. Multi-level clustering to capture all the types of cavities
         1. E.g. Three-level clustering. The first level consists of large surface cavities. The second level consists of small surface cavities. The third level consists of deep cavities (i.e. channels). [DONE]
   2. Can be improved by:
      1. Parameters Tuning:
         1. Utilize min/mean/median/max edge lengths to define eps / search-radius [DONE]
         2. Utilize distibution statistics of distances to define minimum points per cluster. [DONE]
   3. Note: parallelism can be employed, see A Fast Parallel Clustering Algorithm for Large Spatial Databases PDBSCAN [TODO]
2. A Novel Spatial Clustering Algorithm Based on Delaunay Triangulation [TODO]
   1. Note: parallelism can be explored [TODO]
3. OPTICS: ordering points to identify the clustering structure [TODO]

(Spiros) Classify atom types of surface points

1. Clasify atom types of surface points: [DONE]
   1. Using a KD-tree, find the closest atom (from pdb) of every surface point
      1. if a point lies within the radius of the atom, it's assigned its atom type
      2. if a point does not lie with the radius of the atom, it's assigned a value which is considered re-entrace surface point
   2. Note: Shared memory parallelism can be used [DONE]

(Angelos) Classify clusters based on the type of cavity they are:

1. Algorithms:
   1. Ambient occlusion values based classification using mean distance cluster value [TODO]
      1. Note: Shared memory parallelism can be used [TODO]
   2. Analyze the shape of clusters to identify what type of cavity they are. [TODO]
      1. This can potentially be done using machine learning.
      2. Note: Shared memory parallelism can be used [TODO]

Output:

1. Labeled fine surface [DONE]
2. Labeled protein [TODO]
   1. Note: a Paraview plugin may be developed for visualization purposes.

PCITDC vs MOLE

Feature	PCITDC	MOLE
Detect Surface Cavities	YES (more)	YES (less)
Detect Interior Cavities	YES	YES
Visualize Exact Channel Shape	NO	YES
Works on Very Large Molecules	YES	NO (very slow)
Works on Multi-Model Files	YES	NO
Parallelization	YES	NO
Relative Performance	VERY FAST	SLOW
Intuitive Visualization	YES	YES
False Positives	YES (more)	YES (less)
Quantify Major Structural Features	YES	NO

Resources

• ParaView vertex colors

Resources

Software

• UCSF ChimeraX: Next-generation molecular visualization program from the Resource for Biocomputing, Visualization, and Informatics (RBVI).
• VMD: Molecular visualization program for displaying, animating, and analyzing large biomolecular systems using 3D graphics and built-in scripting. Used in conjunction with NAMD for visualizing molecular dynamics simulations.
• NAMD: Parallel molecular dynamics code designed for high-performance simulation of large biomolecular systems. Scales to hundreds of cores for typical simulations and beyond 500,000 cores for the largest simulations.
• MOLE: Universal toolkit for rapid location and characterization of channels and pores in biomacromolecules.
• CAVER: Software tool for the analysis and visualization of tunnels and channels in protein structures.
• RCSB PDB Protein Data Bank: The largest and most popular database of protein structures.

Papers

• Survey: Geometric Protein Cavity Detection Algorithms: Paper surveying the current geometric algorithms for protein cavity detection. This is a great resource for understanding this project.
• MOLE 2.0: Paper describing the MOLE 2.0 protein cavity detection software. Contains a good algorithm that can be used for extracting a good approximation of the surface of a point cloud of atoms.
• MOLE 1.0: Paper describing the MOLE 1.0 algorithm.
• MolAxis: Paper describing the MolAxis protein cavity detection software.
• CAVER 3.0: Paper describing the CAVER 3.0 protein cavity detection software. This algorithm specializes in identifying tunnels.
• Travel Depth: Paper describing the "travel depth" algorithm for protein cavity visualization.
• GPU Protein Cavity Detection using Ambient Occlusion: Paper describing a GPU-based protein cavity detection algorithm that utilizes Ambient Occlusion.
• DBSCAN: Paper describing the DBSCAN clustering algorithm.
• PDBSCAN: Paper describing a parallel version of DBSCAN.

Organizations/Groups

• University of Illinois at Urbana-Champaign: Theoretical and Computational Biophysics Group: Perhaps the most important group in molecular dynamics. Developers of VMD and NAMD. Ran the simulation for the 200 million atom model of COVID-19.
• Oak Ridge National Laboratory (ORNL): Major national lab with active involvement in drug discovery and design using supercomputers.
• Resource for Biocomputing, Visualization, and Informatics (RBVI): Developers of UCSF Chimera and UCSF ChimeraX.
• IST Austria: Edelsbrunner Group: Important group in algorithms, computational geometry, and computational topology. Involved with protein cavity detection in the past and potentially involved now. The head of the group, Edelsbrunner created alpha shapes.
• Texas Advanced Computing Center: Maintainers of the Frontera supercomputer, the fifth fastest supercomputer in the world (2020). The simulation for the 200 million atom model of COVID-19 was executed on this supercomputer.
• WebChemistry: Developed MOLE and many other software tools related to molecular analysis.

Important Links

• COVID-19 200 million atom NAMD simulation
• ORNL COVID-19 spike drug identification

Proteins

• Large Ribosomal Subunit (98,693 atoms)
• COVID-19 Spike Glycoprotein, Closed State (23,694 atoms)
• COVID-19 Main Protease (5,370 atoms)
• COVID-19 NSP15 Endoribonuclease (24,036 atoms)

References

[1] Simões, T., Lopes, D., Dias, S., Fernandes, F., Pereira, J., Jorge, J., Bajaj, C., & Gomes, A. (2017). Geometric Detection Algorithms for Cavities on Protein Surfaces in Molecular Graphics: A Survey. Computer graphics forum : journal of the European Association for Computer Graphics, 36(8), 643–683. https://doi.org/10.1111/cgf.13158

[2] Yaffe, Fishelovitch, Wolfson, Halperin, & Nussinov. (2008). MolAxis: Efficient and accurate identification of channels in macromolecules. Proteins: Structure, Function, and Bioinformatics, 73(1), 72-86.

[3] Coleman, R. G., & Sharp, K. A. (2006). Travel depth, a new shape descriptor for macromolecules: application to ligand binding. Journal of molecular biology, 362(3), 441–458. https://doi.org/10.1016/j.jmb.2006.07.022