What is the correct method of identifying the interacting residues for specific molecular docking?

Question

I'm trying to do a specific molecular docking using Autodock Vina in PyRx. So, one way to do that is to mark all of the interacting residues and make sure the grid box encompasses all of that interacting residues.

So, my protein already has a co-crystalized ligand attached to it. However, I'm not sure what method to use.

After searching online, I found out that I should identify the interacting residues of that ligand. Then, I have opened this protein-ligand complex in 2 applications and 1 web server.

Here are the information of the protein-ligand complex:

• The protein PDB ID is: 5TL8 (Naegleria fowleri CYP51-posaconazole complex)
• The ligand is Posaconazole
• The cofactor and other small molecules: HEME, 1,2-ETHANEDIOL, calcium ion

I tried opening the complex in Discovery Studio Ligplus

and Protein Plus webserver (this server has more interacting residues). All of them have different residues.

Which one should I follow? Should we included the HEME or do we just look at the Posaconazole itself?

Answer 1

Not an answer either. But a series of caveats.

First, HEM parameterisation. For molecular thermodynamics tools you will need to generate a topology file for this bad-boi. The iron has 2 donor bonds and 4 single bonds, one of which with the coordinating cysteine and the other with an oxygen, which will attach the substrate, in this case a steroid —let's ignore the radical state of the activated bound oxygen, but do keep in mind that the partial charges and atom type auto-assigned will be wrong. A cheat would be constraining said molecular oxygen and using a porphyrin ring with no iron: not reality, but would make the experiment work, is as physically questionable as any iron-bound topology you'll get and the oxygen is the key player here.

Second, cytochromes act on greasy substrates, there are few if any polar interactions generally, so constraining to polymer residues may not be that useful. Specificity is controlled by the gating and entranceway (think of trucks reversing into a loading bay), which is annoying the gating loops are floppy loops that are generally not crystallised, so tethered MD from a bound to unbound state (undocking) is sometimes done. There's a lot of literature on P450 and P411 cytochromes, but often it's not as clear cut and simple as say a hydrolase.

However, if your substrate you want placed is similar to one present you may be interested in a tool of mine Fragmenstein, which can be used to convert a compound to a similar one (stripping the haem and adding it later: I had a Rosetta params file for HEM somewhere but not sure where it is nowadays)

Answer 2

not an answer, just more inputs, but https://proteins.plus/5tl8#poseview

gives me this:

references:

ProteinsPlus: interactive analysis of protein–ligand binding interfaces

Molecular complexes at a glance: automated generation of two-dimensional complex diagrams

Drawing the PDB: Protein−Ligand Complexes in Two Dimensions

from this last paper :

Hydrogen bonds are implemented mainly following the measures published by Desiraju and Steiner in 2001.14 The optimal distance between two atoms connected by a hydrogen bond is set to 1.9 Å with a tolerance of 0.5 Å. Additional to this measure, the acceptor−hydrogen−donor angle must not fall below 120°. Hydrogen atoms that are bound to a noncarbon atom are treated as hydrogen donor candidates. Potential acceptor atoms are either nitrogen, oxygen, or sulfur atoms provided that they are uncharged or negatively charged and that their surface is accessible.

Metal interactions are calculated between metal atoms embedded in the receptor and metal acceptor atoms, which are identical to the hydrogen bond acceptor atoms. Their geometry is based on the calculated coordination geometry of the metal.15 Each coordination point that is not occupied by a receptor atom is checked for close ligand atoms, and the maximal distance deviation is set to 0.8 Å. If no geometry can be calculated, a sphere with a radius of 2 Å is placed around the metal. In this case, all atoms lying on the sphere, again with a tolerance of 0.8 Å, are regarded as potential interaction partners.

In contrast to the formerly mentioned interactions, hydrophobic contacts are estimated based on the distance between two hydrophobic atoms only. They are visualized not by a dashed interaction line but by drawing the label of the contacting residue and a spline segment denoting the hydrophobic part of the ligand. Because many atoms typically form a hydrophobic subpocket, this representation reflects the interaction geometry better. A prerequisite for a hydrophobic contact is that at least three hydrophobic ligand atoms lie in the range of the currently examined receptor residue. Hydrophobic atoms are in this context carbon atoms with accessible surface and halogens. The maximum distance is set to the sum of the van der Waals radii of atoms in question and a tolerance of 0.8 Å.

EDITED :

Adding more on a side note I came accross this paper, I am completely naive about the topic but the title :

“Observing Noncovalent Interactions in Experimental Electron Density for Macromolecular Systems: A Novel Perspective for Protein-Ligand Interaction Research.” Journal of chemical information and modeling vol. 62,7 (2022): 1734-1743. doi:10.1021/acs.jcim.1c01406

is kind of intriguing, I am not an expert and cannot vauch for its strictness but to me sounds like a good why to validate an interaction like seeing it even if in abstract terms:

... Non-covalent interactions (NCI) can be observed by pinpointing the ED saddle point, i.e. (3,−1) critical points and further quantified by measuring ED deviation from a homogeneous electron distribution using density and its first derivative (s = [1/(2(3π2)1/3)]|∇ρ|/ρ4/3)....

Please feel free to add to this, its not my field.

https://www.biorxiv.org/content/10.1101/2022.01.24.468575v2.full

Answer 3

I think I've found the answer to my own question.

So, I tried using this PRANK web server for predicting the binding site of my protein. The results is suprisingly very similar to the one from the Protein Plus web server. So, I think it's better to use the one from Protein Plus.

You can see the results on the link below.

References:

1. https://prankweb.cz/viewer?id=5tl8&database=v3-conservation-hmm&created=2022-01-27T07:50:59&lastChange=2022-01-27T07:50:59%20%20&structureName=structure.pdb&predictionName=5TL8&predictedStructure=undefined
2. https://proteins.plus/5tl8