Given the 3D structure of a molecule (e.g. as *.mol) and a probe radius r, I want to calculate the solvent excluded surface (SES) and the excluded volume enclosed by it. Note the difference between SAS and SES:
I'm looking for a command line program, C/C++ or Python library that can accomplish such a task.
MSMS is a rather old (1996) command line tool that computes the Solvent Excluded Surface (SES). It is wrapped in a python module in mgtool or in openstructure which may be more convenient to use.
EDTSurf[1] is a common and reasonably well regarded approach. The source code is available at https://zhanglab.ccmb.med.umich.edu/EDTSurf/
I tested the mac version and it worked with the following:
./EDTSurf_mac -i kras.pdb
which outputs kras.ply. A slight hiccup is the website describes different command line arguments - but the default is the SES. You'll see they have windows and linux versions, too.
Assuming you're comfortable with python, you could then visualize the output with Open3D:
import open3d as o3d
mesh = o3d.io.read_triangle_mesh("./kras.ply")
o3d.visualization.draw_geometries([mesh])
Now, for volume. Volume of the SES is difficult to calculate exactly, but the volume of a mesh is trivial to calculate. Open3D has a get_volume method:
mesh.get_volume()
unfortunately that froze my jupyter notebook every time - you may have better luck with a fresh install.
But the following will also work:
import numpy as np
verts = np.array(mesh.vertices)
faces = np.array(mesh.triangles)
def signed_vol_of_triangle(p1, p2, p3):
v321 = p3[0]*p2[1]*p1[2]
v231 = p2[0]*p3[1]*p1[2]
v312 = p3[0]*p1[1]*p2[2]
v132 = p1[0]*p3[1]*p2[2]
v213 = p2[0]*p1[1]*p3[2]
v123 = p1[0]*p2[1]*p3[2]
return (1 / 6)*(-v321 + v231 + v312 - v132 - v213 + v123)
def make_vol(pts):
return signed_vol_of_triangle(pts[0], pts[1], pts[2])
print(sum([make_vol(v) for v in verts[faces]]))
In short, this loops through each triangle in the mesh and calculates the signed volume of the triangular prism with the fourth point at the origin. The theory behind this is in [2], and I used the implementation in this SO answer which you might want to also consult: https://stackoverflow.com/a/1568551/3089865
While I haven't rigorously tested that volume-calculating approach for correctness, it returns 21058.959617488006 angstrom cubed, which seems reasonable to me (~21 nm^3) given KRAS fits in a ~4x3x3nm box.
[1] Xu, Dong, Hua Li, and Yang Zhang. "Protein depth calculation and the use for improving accuracy of protein fold recognition." Journal of Computational Biology 20.10 (2013): 805-816.
[2]Zhang, Cha, and Tsuhan Chen. "Efficient feature extraction for 2D/3D objects in mesh representation." Proceedings 2001 International Conference on Image Processing (Cat. No. 01CH37205). Vol. 3. IEEE, 2001.
Gepol93 is a FORTRAN library which computes three kinds of envelope surfaces :
There is also a python library called Pychemovality (which is authored by me) that calculates the solvent excluded surface area (SESA) in order to calculate ovality of molecules.
This python library uses a modified version of GEPOL93 under the hood. This library outputs the SESA and corresponding volume.
Also if you are interested in using the modified GEPOL93 Fortran script directly, then I have included compilation instructions for it on the pychemovality repo here.
