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. (Angelos) Using vtkPDBReader, read the PDB file that included the protein. [DONE] [Spring 2020]
4. Challenges:
   1. (Angelos) Upgrade vtkPDBReader to be able to manage multimodel proteins [DONE] [Summer 2020]
   2. (Spiros) Fix vtkPDBReader error of setting wrong atom numbers [DONE] [Summer 2020]
   3. (Spiros) Set up vtkPDBReader to utilize The Blue Obelisk Dataset for the Covalent radius and Van-der-Waals radius [DONE] [Summer 2020]
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) Extract unique protein Chains

1. Input: A Polydata of the protein
2. Output: Unique protein chains without het-atoms and Symmetric Information
3. Algorithmic step:
   1. Implement a filter that removes the het-atoms and creates the polydata of the protein chains. [DONE] [Fall 2020]
   2. Detect the unique protein chains using Symmetry detection techniques (Kabsch algorithm) [DONE] [Spring 2021]
   3. Store the Symmetry information (rotation and translation). [DONE] [Spring 2021]
      1. Note: This way we can mesh each unique chain individually

(Spiros) Sort Atoms By Radius

1. Input: A Polydata of a unique protein chain without Het atoms
2. Output: A Polydata of a unique protein chain 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] [Summer 2020]
      1. Note: Shared memory parallelism can be employed using vtkSMPTools [DONE] [Summer 2020]

(Spiros/Αngelos) Compute the Solvent Excluded Image

1. Input: A Polydata of a unique protein chain without Het atoms sorted by radius in descending order, probe radius, spacing
2. Output: Solvent Excluded Image of a chain
3. Algorithmic steps:
   1. (Spiros) Calculate the optimal bounding box of the protein chain using CGAL's implementation [DONE] [Spring 2021]
   2. (Spiros/Angelos) Create the Solvent Accessible image by creating spheres for each atom that their radius is: Van der Waals radius + probe radius. [DONE] [Summer 2020]
      1. Note: Shared memory parallelism can be employed [DONE] [Summer 2020]
   3. (Angelos) Perform Euclidean Distance Transform on the Solvent Accessible image using the itk::SignedMaurerDistanceMapImageFilter. [DONE] [Summer 2020]
      1. Note 1: Shared memory parallelism can be employed [DONE] [Summer 2020]
      2. Note 2: GPU parallelism can be employed [TODO]
   4. (Spiros) Extract Solvent Excluded image from the EDT result by setting: [DONE] [Summer 2020]
      1. All pixel values < probe-radius - truncationError[atomType] - paddingDistance = 0
      2. All pixel values >= probe-radius - truncationError[atomType] - paddingDistance = 1
      3. Note 1: Shared memory parallelism can be employed [DONE] [Summer 2020]
      4. Note 2: GPU parallelism can be employed [TODO]
4. Challenges:
   1. (Spiros) Handle the calculation of truncationError[atomType] in a thread-safe way by sorting the atoms by radius in descending order [DONE] [Summer 2020]
      1. Truncation error is caused because of the spacing of the image.
   2. (Spiros/Angelos) Select a paddingDistance value to help PODM3D produce an accurate surface [DONE] [Summer 2020]
   3. (Spiros) Close holes that should not exist using the distance of the pixel point vs the actual atom [DONE] [Fall 2020]
   4. (Spiros) Create smaller images using the concept that is described in the EDTSurf paper [TODO]

(Spiros) Compute the Solvent Excluded Surface (SES)

1. Input: Solvent Excluded Image of a chain
2. Output: Solvent Excluded Surface of a chain
3. Algorithmic steps:
   1. Using PODM3D, produce a volumetric mesh of the image. [DONE] [Summer 2020]
      1. Note: Shared memory parallelism can be employed [DONE] [Summer 2020]
   2. Using PODM3D, extract the volumetric mesh. [DONE] [Summer 2020]
      1. Note: Shared memory parallelism can be employed [DONE] [Summer 2020]
   3. Using vtkDataSetSurfaceFilter, extract the surface of volumetric mesh. [DONE] [Summer 2020]
      1. Note 1 this filter is used because when PODM is extracting the surface, its triangles are oriented the wrong way.
      2. Note 2: Shared memory parallelism can be employed [TODO]
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. Select a default delta value. [DONE] [Fall 2020]
   2. Fix volume grading crashing. [TODO]
   3. Surface needs to be smoothed. [TODO]

(Spiros) Extract Atom Types of surface points

1. Input: A Polydata of a unique protein chain surface and a Polydata of the unique protein chain without Het atoms sorted by radius in descending order
2. Output: Atom Types of surface points of the unique protein chain surface
3. Algorithmic step:
   1. Using a KD-tree, for each surface point, find the closest atom: [DONE] [Summer 2020]
      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 employed [DONE] [Summer 2020]
4. Challenges:
   1. Find an ideal offset distance value between the surface points and their closest atoms to improve accuracy. [DONE] [Fall 2020]

(Spiros) Generate the Symmetric Protein chain surfaces

1. Input: A vector of the unique protein chain surfaces and the Symmetry information
2. Output: A vector of all the protein chain surfaces
3. Algorithmic step:
   1. Implement the SymmetricProteinChainsSurface filter that generates the non unique protein chain surfaces using the Symmetry information. [DONE] [Spring 2021]
      1. Note: Shared memory parallelism can be employed [DONE] [Spring 2021]

(Spiros) Assign Model and Chain arrays on SESs

1. Input: A vector of all the protein chain surfaces & model and chains information
2. Output: A vector of all the protein chain surfaces with model and chain arrays
3. Algorithmic step:
   1. Implement the AssignModelAndChain filter that assigns the model and chain arrays for each protein chain surface. [DONE] [Spring 2021]

(Spiros) Merge Protein Chain SESs

1. Input: A vector of all the protein chain surfaces with model and chain arrays
2. Output: Solvent Excluded Surface
3. Algorithmic step:
   1. Implement the MergeSurfaces filter that merges all the protein-chains surfaces into one and forms the solvent excluded surface of the protein. [DONE] [Spring 2021]
      1. Note: Shared memory parallelism can be employed [DONE] [Spring 2021]

(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. (Angelos/Spiros) 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] [Summer 2020]
      1. Note: Shared memory parallelism can be employed [DONE] [Summer 2020]
   2. (Angelos/Spiros) Using the CUDA and Optix library, calculate normals of surface points and then calculate Ambient occlusion values for all the surface points. [DONE] [Spring 2021]
      1. Note: GPU parallelism can be employed [DONE] [Spring 2021]
4. Challenges:
   1. (Spiros/Angelos) Eliminate the noise that is created due to the non-smoothed PODM3D surface by assigning for each vertex the average ambient occlusion value of its neighbors + itself. [DONE] [Fall 2020]
5. Ideas that might be used:
   1. Use the symmetry information to minimize the memory footprint for CPU and GPU [TODO]

(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] [Summer 2020]
      1. Note: Shared memory parallelism can be employed [DONE] [Summer 2020]

(Angelos/Spiros) Cavity Detection

1. Input: Solvent Excluded Surface, Ambient Occlusion values of surface points, Surface Statistics
2. Output: Protein Cavities
3. Algorithmic step:
   1. (Angelos) 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] [Summer 2020]
      1. Note: shared memory parallelism can be employed [TODO]
         1. Idea 1: A Fast Parallel Clustering Algorithm for Large Spatial Databases PDBSCAN
         2. Idea 2: Speculative Execution using PREMA
   2. (Angelos) Implement a rendering framework to visualize the clustering process using VTK. [DONE] [Spring 2021]
4. Challenges:
   1. Tune the clustering parameters:
      1. (Spiros) Search-radius = 2 * minEdgeLength of the surface. [DONE] [Spring 2021]
      2. (Angelos) Esp = 8. [DONE] [Spring 2021]
      3. (Spiros) Min cavity size = cavityFactor * delta^-2.178. [DONE] [Spring 2021]
         1. Default value of cavityFactor = 20.
   2. Limit the expansion of the cavity by using:
      1. (Spiros) A relative threshold for the difference an occlusion value from the average occlusion value of the cluster. [DONE] [Fall 2020]
      2. (Angelos) Define weights for the occlusion values, and convert the mean to a weighted mean. [DONE] [Spring 2021]
   3. Define occlusion values bounds to detect different types of cavities [TODO]
5. Ideas that might be used:
   1. OPTICS: ordering points to identify the clustering structure [TODO]
   2. Use symmetry information to quickly detect the per-model cavities [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. (Spiros/Angelos) For each surface point: [DONE] [Summer 2020]
      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] [Summer 2020]
4. Challenges:
   1. (Angelos/Spiros) Create reader of the OpenDx file format that is used as output by APBS [DONE] [Summer 2020]

(Spiros) Select Clustered cavity points

1. Input: Labeled Solvent Surface Surface
2. Output: Labeled clustered cavity points
3. Algorithmic step:
   1. Implement a SelectPoints filter. [DONE] [Spring 2021]
   2. Select only the points that are clustered at the cavity detection step. [DONE] [Spring 2021]
      1. Note: this is done to improve the speed of writing the output

(Angelos/Spiros) Write Results:

1. Input: Labeled clustered cavity points or Labeled Solvent Surface Surface
2. Output: solventSurface.vtp
3. Algorithmic step:
   1. (Angelos) Using a vtkXMLPolyDataWriter, write the output file. [DONE] [Spring 2020]
   2. (Angelos/Spiros) Create a writer, just for the cavities to minimize the IO. [TODO]
      1. Idea: See Fpocket.
4. Ideas that might be used:
   1. Note: a Paraview plugin may be developed for visualization purposes. [TODO]

Future Work

(Angelos/Spiros) Classify cavity type

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]

Results

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

Paper and Presenation

Presentation

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