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: Geometric Protein Cavities Detection (GPCD)

The Pipeline contains the following high-level algorithmic steps:

(Angelos/Spiros) Read protein

1. Input: PDB file
2. Output: A Polydata of the protein
3. Algorithmic Step:
   1. Using vtkPDBReader, read the PDB file that included the protein. [DONE]
4. Challenges:
   1. Upgrade vtkPDBReader to be able to manage multimodel proteins [DONE]
   2. Fix VTKPDBReader error of setting wrong atom numbers [DONE]
   3. Set up vtkPDBReader to utilize The Blue Obelisk Dataset for the Covalent radius and Van-der-Waals radius [DONE]
5. Ideas that might be used:
   1. Improve The Blue Obelisk Dataset by using: [TODO]
      1. updated Van-der-Waals radius values provided by Los Alamos National Laboratory
      2. updated Covalent radius values provided by Cambridge Structural Database

(Spiros/Angelos) Remove Het Atoms

1. Input: A Polydata of the protein
2. Output: A Polydata of the protein without Het atoms
3. Algorithmic step:
   1. Using vtkSelectionFilter, we extract all the non-Het atoms. [DONE]

(Spiros/Angelos) Sort Atoms By Radius

1. Input: A Polydata of the protein without Het atoms
2. Output: A Polydata of the protein without Het atoms sorted by radius in descending order
3. Algorithmic step:
   1. Using vtkSortFieldData we sort the non-Het atoms by radius in descending order. [DONE]
      1. Note: Shared memory parallelism can be employed [DONE]

(Spiros) Extract protein Chains

1. Input: A Polydata of the protein without Het atoms sorted by radius in descending order
2. Output: A vector of Polydata of the protein chains without Het atoms sorted by radius in descending order
3. Algorithmic step:
   1. Using vtkSelectionFilter, we extract all the protein-chains. [DONE]
      1. Note: This way we can mesh each chain individually

(Spiros/Αngelos) Compute the Solvent Excluded Image

1. Input: A Polydata of a chain of the protein without Het atoms sorted by radius in descending order, probe radius, spacing
2. Output: Solvent Excluded Image of chain
3. Algorithmic steps:
   1. Create the Solvent Accessible image by creating spheres for each atom that 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 Solvent Accessible image. [DONE]
      1. Note: Shared memory parallelism can be employed [DONE]
   3. Extract Solvent Excluded image from the EDT result by setting: [DONE]
      1. All pixel values < probe-radius - truncationError[atomType] - paddingDistance = 0
      2. All pixel values >= probe-radius - truncationError[atomType] - paddingDistance = 1
      3. Note: Shared memory parallelism can be employed [DONE]
4. Challenges:
   1. Handle the calculation of truncationError[atomType] in a thread-safe way by sorting the atoms by radius in descending order [DONE]
      1. Truncation error is caused because of the spacing of the image.
   2. Select a paddingDistance value to help PODM3D produce an accurate surface [DONE]
   3. Close holes that should not exist (based on other surfaces) [DONE]

(Spiros) Compute the Solvent Excluded Surface (SES)

1. Input: Solvent Excluded Image of chain
2. Output: Solvent Excluded Surface of chain
3. Algorithmic steps:
   1. Using PODM3D, produce a volumetric mesh of the image. [DONE]
      1. Note: Shared memory parallelism can be employed [DONE]
   2. Using PODM3D, extract the surface of the volumetric mesh. [DONE]
      1. Note: Shared memory parallelism can be employed [DONE]
4. Advantages:
   1. Can control the number of elements on the surface
   2. the edge lengths bounds are relatively small which helps in the clustering.
5. Challenges:
   1. Using vtkDataSetSurfaceFilter, fix the issue of surface triangles being oriented the wrong way. [DONE]
   2. Select a default delta value. [DONE]
   3. Surface needs to be smoothed. [TODO]

(Spiros/Angelos) Merge Protein Chain SESs

1. Input: A vector of the Solvent Excluded Surface of chains
2. Output: Solvent Excluded Surface
3. Algorithmic step:
   1. Using vtkAppendPolyData, we merge all the protein-chains surfaces polydata. [DONE]

(Angelos/Spiros) Compute ambient occlusion of the SES

1. Input: Solvent Excluded Surface
2. Output: Ambient Occlusion values of surface points, normals of surface points
3. Algorithmic step:
   1. Using IGL Library, which uses Embree and Eigen, calculate normals of surface points and then calculate Ambient occlusion values for all the surface points [DONE]
      1. Note: Shared memory parallelism can be used [DONE]
      2. Note: GPU parallelism can be used [TODO]
4. Ideas that might be used:
   1. Assign the ambient occlusion value of a triangle by calculating the average of its points. [TODO]
      1. This can be used for a cell-based clustering method in order to eliminate the noise values caused by the PODM3D's solvent accessible surface.

(Spiros) Compute Surface Statistics

1. Input: Solvent Excluded Surface
2. Output: Statistics of Surface Triangles' Edges
3. Algoritmic step:
   1. Compute min/mean/median/max edge lengths of the surface triangles [DONE]
      1. Note: Shared memory parallelism can be used [DONE]

(Angelos) Detect protein cavities

1. Input: Solvent Excluded Surface, Ambient Occlusion values of surface points, Surface Statistics
2. Output: Atomic Cavities
3. Algorithmic step:
   1. Using a modified version of the clustering method DBSCAN, that finds neighbors based on their density and ambient occlusion values, form clusters that define protein cavities. [DONE]
      1. Note: shared memory parallelism can be employed, see A Fast Parallel Clustering Algorithm for Large Spatial Databases PDBSCAN [TODO]
4. Challenges:
   1. Tune the clustering parameters eps / search-radius using the surface statistics. [DONE]
   2. Handle non-deterministic Solvent Excluded Surface of PODM3D. [TODO]
5. Ideas that might be used:
   1. A Novel Spatial Clustering Algorithm Based on Delaunay Triangulation [TODO]
   2. OPTICS: ordering points to identify the clustering structure [TODO]

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

1. Input: Solvent Excluded Surface, Atomic Cavities, Atom Types of surface points
2. Output: Labeled atomic Cavities
3. Algorithm steps:
4. Ideas that might be used:
   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]

(Spiros) Classify Atom Types of surface points

1. Input: Solvent Excluded Surface, a Polydata of the protein without Het atoms sorted by radius in descending order
2. Output: Atom Types of surface points
3. Algorithmic step:
   1. Using a KD-tree, for each surface point, find the closest atom: [DONE]
      1. If the surface point is within the radius of the closest atom, assign the atom type of the closest point to the surface point
      2. Else, assign 0, which means the surface point is part of the re-entrance surface.
      3. Note: Shared memory parallelism can be used [DONE]
4. Challenges:
   1. Find an ideal offset distance value between the surface points and their closest atoms to improve accuracy. [TODO]
      1. Idea: use the truncationError[AtomType]
5. Ideas that might be used:
   1. create a different kd-tree for each protein chain [TODO]

(Spiros/Angelos) Map Electrostatics Potentials on Surface Points

1. Input: Solvent Excluded Surface, Electrostatic Potentials Image ESP of APBS, normals of surface points, probe radius, lowerThreshold-upperThreshold
2. Output: Electrostatic values of surface points
3. Algorithmic steps:
   1. For each surface point: [DONE]
      1. Scale the point using normals, and probe radius to the solvent-accessible surface,
      2. Get the value of the closest pixel from the ESP to the scaled point:
         1. If value > upperThreshold, set electrostatic value of point = upperThreshold
         2. Else if, value < lowerThreshold, set electrostatic value of point = lowerThreshold
         3. Else get the interpolated value of the scaled points from the image:
            1. If interpolated value > upperThreshold, set electrostatic value of point = upperThreshold
            2. Else if, interpolated value < lowerThreshold, set electrostatic value of point = lowerThreshold
      3. Note: Shared memory parallelism can be used [DONE]
4. Challenges:
   1. Create reader of the OpenDx file format that is used as output by APBS [DONE]

(Angelos/Spiros) Write Labeled Solvent Surface Surface:

1. Input: Labeled Solvent Surface Surface
2. Output: solventSurfaceSurface.vtp
3. Algorithmic step:
   1. Using a vtkXMLPolyDataWriter, write the output file. [DONE]
4. Ideas that might be used:
   1. Note: a Paraview plugin may be developed for visualization purposes. [TODO]

GPCD vs MOLE

Feature	GPCD	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.
• Coronavirus envelope protein: Paper describing research on the envelope protein of coronaviruses.
• COVID-19 virion model: Paper describing a detail model of the SARS-CoV-2 virus particle.
• COVID-19 simulation: Paper describing large-scale COVID-19 simulations.

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