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]

Candidate Parallel Methods

1. PODM Tessellation Using Van der Waals Surface Image (Angelos)

The first method that could be used for detecting protein cavities in parallel involves using PODM to create a tessellation of an image representation of the Van der Waals surface of a protein. The Van der Waals surface is an accurate representation of the surface of a protein and represents each atom as a sphere of a radius corresponding to its Van der Waals radius.

The steps for this method are outlined below:

1. Compute the Van der Waals surface (Angelos): This polygonal surface can be computed using UCSF Chimera. The output from this is an STL file containing the polygonal surface mesh. However, because the implementation is sequential, this step can present a major bottleneck in the overall method. In order to speed the process up, a custom parallel implementation to generate this surface will be necessary in the future.
2. Generate image of the Van der Waals surface (Angelos): Because PODM requires an image as input, the polygonal Van der Waals surface needs to be converted to an image. This is done by "filling" the polygonal surface with a constant value. The segmentation module of 3D Slicer can be used for doing this, but as before, this step can present a major bottleneck because the implementation is sequential. A custom parallel implementation is once again necessary.
3. Compute Delaunay tessellation using PODM: Once we have the image of the Van der Waals surface, we use PODM to compute the Delaunay tessellation in parallel. A low enough delta should be used to ensure that the surface is as accurately captured as possible.
4. Compute Voronoi tessellation using PODM: The next step is to compute the Voronoi tessellation as the dual of the Delaunay tessellation in parallel. The Voronoi tessellation is useful for quickly identifying tunnels using Dijkstra's algorithm.
5. Identify cavities in parallel (Angelos & Spiros): With the Delaunay and Voronoi tessellations computed, the cavities now need to be identified in parallel. It is crucial that this step is done in parallel as no other method implements "proper" parallelization. This is by far the most challenging step.
6. Output detected cavities as tessellations: The last step is to output the detected cavities as a tessellated mesh. It is important to only include the identified cavities in this mesh by representing them as tetrahedra. This way, researchers can overlay the normal protein structure with the tessellations of the cavities and easily see where they are.

Figure 4: The Van der Waals surface of the COVID-19 spike glycoprotein (6vxx).

Figure 5: A tessellation based on the Van der Waals surface of the COVID-19 spike glycoprotein (6vxx). A highly-accurate surface tessellation is produced.

2. PDT Tessellation Using Point Cloud of Atoms (Spiros)

The second method involves giving PDT/PODM a point cloud of atoms and creating a tessellation from that. This would require that the Delaunay triangulation is generated and a method for computing an accurate surface for the point cloud is employed. Computing this surface is one of the most challenging steps of this method. A method that could be employed is α-shapes, but this has the disadvantage of not guaranteeing connectivity between tetrahedra. Another method probably needs to be used that can guarantee connectivity between the tetrahedra (such as the one from MOLE 2.5).

The steps for this method are outlined below:

1. Read PDB file and make point cloud of atoms (Angelos): The first step would be to read an input PDB file for a protein and create a point cloud of atoms. This involves either writing a basic PDB file reader or finding one we can use. However, writing a basic one is better because it skips the overhead of other full-featured readers as we only care about the atomic centers and potentially the type of atoms.
2. Compute Delaunay triangulation using PDT/PODM: The second step is to compute the Delaunay triangulation from the point cloud in parallel. This will produce an inaccurate surface that needs to be refined.
3. Compute accurate surface (Spiros): The Delaunay triangulation on a point cloud of atoms produces a surface that does not reflect the actual surface of the protein. As a result, an accurate surface needs to computed by refining the triangulation. This would involve removing tetrahedra until an accurate surface is found. This should be done in parallel.
   1. MOLE method:
      1. The molecular surface is approximated by iterative removal of boundary tetrahedra from the outermost layers (i.e., tetrahedra found at the interface between the molecule and the external environment).
      2. Boundary tetrahedra produced by the triangulation are removed in this step if they are sufficiently large to fit a sphere with a given probe radius
         1. Note: A tetrahedron T fits a sphere S with probe radius r if the center C of sphere S can be placed inside the tetrahedron and the distance to all vertices of T is greater than or equal to the sum of r, with the Van der Waals radius of an atom corresponding to the given vertex.
            1. Explanation: Given a boundary tetrahedron, if a Sphere S with center = Incenter C, and radius = probe radius
      3. Next, tetrahedra that are too small to fit a sphere with interior radius are removed.
         1. Explanation: Given a non-boundary tetrahedron, if the radius of its inscribed sphere is smaller than (probe?) radius, then the tetrahedron needs to be removed.
      4. Remaining tetrahedra form one or more connected components. We call the components that contain at least one tetrahedron on the molecular surface cavity diagrams.
      5. It should be noted here that the cavity diagram is a purely geometrical concept to help identify regions of space (volume) that can contain tunnels and only very loosely corresponds to the cavities shown in Figure 1B).
4. Compute Voronoi tesselation: Once the accurate surface is determined the Voronoi tessellation must be computed. This step could be done before the accurate surface is computed, but we would also have to refine the Voronoi if that is to happen.
5. Identify cavities in parallel (Angelos & Spiros): With the Delaunay and Voronoi tessellations computed, the cavities now need to be identified in parallel. It is crucial that this step is done in parallel as no other method implements "proper" parallelization. This is by far the most challenging step.
6. Output detected cavities as tessellations: The last step is to output the detected cavities as a tessellated mesh. It is important to only include the identified cavities in this mesh by representing them as tetrahedra. This way, researchers can overlay the normal protein structure with the tessellations of the cavities and easily see where they are.

Figure 6: The atomic point cloud of the COVID-19 spike glycoprotein (6vxx). There are 23,694 atoms (points) in this image.

Figure 7: An unrefined Delaunay triangulation of the COVID-19 spike glycoprotein (6vxx). Generated using the Delaunay3D filter in ParaView.

Figure 8: A Delaunay triangulation using α-shapes with α=3 of the COVID-19 spike glycoprotein (6vxx). Using α-shapes we can get a better approximation of the surface. Generated using the Delaunay3D filter in ParaView.

3. Surface-based Using Van der Waals Surface Volumetric Image

The third method involves writing a new software for doing protein cavity detection. This software would use a volumetric image (e.g. of type .mha) representing the Van der Waals surface of a protein. Uses perhaps the most faithful representation of the Van der Waals surface, which is the ground truth for a molecular surface. As long as the Van der Waal surface generation algorithm is fast, this method has the potential of being both accurate and fast as it eliminates the overhead of tessellation. Good parallelization is important for this method because there are many more voxels than there are tetrahedra in a tessellation-based method. Because this method operates on voxels, it is a candidate for future GPU parallelization. This method is also a strong candidate for deep learning, much more than the tessellation based methods, because its output can easily be used for training a deep learning model in the future.

The steps for this method are outline below:

1. Compute the Van der Waals surface: This polygonal surface can be computed using UCSF Chimera. The output from this is an STL file containing the polygonal surface mesh. However, because the implementation is sequential, this step can present a major bottleneck in the overall method. In order to speed the process up, a custom parallel implementation to generate this surface will be necessary in the future.
2. Generate image of the Van der Waals surface: Because PODM requires an image as input, the polygonal Van der Waals surface needs to be converted to an image. This is done by "filling" the polygonal surface with a constant value. The segmentation module of 3D Slicer can be used for doing this, but as before, this step can present a major bottleneck because the implementation is sequential. A custom parallel implementation is once again likely necessary in the future.
3. Identify cavities in parallel: Using the image voxels, the cavities now need to be identified in parallel (potentially using GPU parallelization in the future). It is crucial that this step is done in parallel as no other method implements "proper" parallelization. This is by far the most challenging step.
4. Output detected cavities as volumetric image: The last step is to output the detected cavities as a volumetric important. It is important to only include the identified cavities in this image. This way, researchers can overlay the normal protein structure with the volume of the cavities and easily see where they are.

Figure 9: Volumetric image of the Van der Waals surface of the COVID-19 spike glycoprotein (6vxx). There are 15.5 million voxels in this volume.

Figure 10: A slice of the volumetric image of the Van der Waals surface of the COVID-19 spike glycoprotein (6vxx). Potential cavities and voids are clearly visible.

Figure 11: An animation of a thicker slice of the volumetric image of the Van der Waals surface of the COVID-19 spike glycoprotein (6vxx). Potential cavities and voids are clearly visible.

Challenges

There are several challenges that need to be addressed if we are to have a good parallel solution for detecting protein cavities.

1. Expensive Pre-Processing

The parallel methods involving the Van der Waals surface, especially the one with the image input have expensive pre-processing. For large molecules, getting this surface can be slow because existing implementations are sequential and written in Python. We need to come up with a way quickly generate this surface and its image representation. This pre-processing does not have to be done if we operated on the point cloud using method 2 described previously, leading to a major speedup in overall performance.

Figure 12: The Van der Waals surface of the COVID-19 spike glycoprotein (6vxx).

2. Massive Parallelism (Stretch Goal)

The solution we will implement may need to support massive parallelism, ideally using MPI. The reason for this is that molecular dynamics simulations can potentially involve a massive number of atoms. If we consider the 200 million coronavirus model, then a simulation containing 5000 frames, would involve a total of 1 trillion atoms. If we are to detect the cavities in every frame in real-time, a massively parallel solution is necessary. There is no doubt about why 250,000 CPU cores were required for the COVID-19 simulation.

Figure 13: An animation produced by a molecular dynamics simulation.

Figure 14: An animation produced by the COVID-19 molecular dynamics simulation.

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.
• CAVER 3.0: Paper describing the CAVER 3.0 protein cavity detection software. This algorithm specializes in identifying tunnels.

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

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