Physics-based pose predictions guided by native ligands
SchrödingerRelease2016.3/WaterMapv2.7/Schrödinger Release 2016.2/Maestro10.6/ProteinPreparationWizard/Primev4.4/Impactv71014/LigPrepv38014/Desmondv4.6/Glidev7.1/Jaguarv9.2/Macromodelv11.2/Phasev4.7/Jaguarv9.2/Epikv3.6
Using Schrödinger's Protein Preparation Wizard, we converted the raw PDB structures into all-atom, fully prepared protein models. We followed the workflow of the Protein Preparation Wizard using default parameters set by Schrodinger but also enabled Prime prediction of the missing side chains and missing loops.
The protein preparation took place using the following procedure from within the Protein Preparation Wizard, Schrödinger Release 2016-2, Maestro 10.6:
1) Automatically imported full PDB files from the PDB website
2) Automatically added missing hydrogen atoms
3) Enumerate bond orders to HET groups
4) Removed or kept co-crystallized water molecules as described below (see following section)
5) Capped protein termini with ACE and NMA residues
6) Highlighted residues with missing atoms or multiple occupancies
7) Pre-processed structures for Prime, Schrödinger's program for protein structure prediction
8) Determined the most likely ligand protonation state as well as the energy penalties associated with alternate protonation states with Epik
9) Determined optimal protonation states for histidine residues
10) Corrected potentially transposed heavy atoms in arginine, glutamine, and histidine side chains
11) Optimize the protein's hydrogen bond network by means of a systematic, cluster-based approach, which greatly decreases preparation times
12) Perform a restrained minimization that allows hydrogen atoms to be freely minimized, while allowing for sufficient heavy-atom movement to relax strained bonds, angles, and clashes using Impact
Schrödinger Release 2016-3: Schrödinger Suite 2016-3 Protein Preparation Wizard; Epik, Schrödinger, LLC, New York, NY, 2016; Impact, Schrödinger, LLC, New York, NY, 2016; Prime, Schrödinger, LLC, New York, NY, 2016.
Sastry, G.M.; Adzhigirey, M.; Day, T.; Annabhimoju, R.; Sherman, W., "Protein and ligand preparation: Parameters, protocols, and influence on virtual screening enrichments," J. Comput. Aid. Mol. Des., 2013, 27(3), 221-234
Van der Waals scaling factor for ligand= 0.8
Ligand sampling= flexible
number of poses=3
rotatables groups of Ser, Thr, Tyr were allowed to rotate
The pose prediction parameters and pose prediction methodology is described below.
1) Given the abundance of structural FXR data, we decided to utilize this information for the prediction of Downloaded 28 FXR structures that we found available in the PDB website:
We also downloaded the APO STRUCTURE from D3R
2) Based on superposition of these proteins and clustering of water molecules in the binding site, it was decided which waters to keep in the protein preparation stage which is described above.
This choice was confirmed by a WaterMap calculation on selected proteins: APO structure provided by D3R, PDB ID: 3OMM, 3RUT
3) Protein Preparation
The default parameters from Maestro Protein Preparation Module were used, except that we also added filling in missing loops and missing side chains.
Water was kept based on WaterMap calculations and alignments as explained above.
If a ligand had a carboxy group, the deprotonated form was kept.
The Maestro Protein Preparation Module was used to prepare the 29 crystal structures.
4) Ligand Preparation
LigPrep was used to prepare the ligands provided by D3R.
The LigPrep panel wass used to to set up and start ligand preparation calculations. LigPrep produced the corresponding low-energy 3D structures. LigPrep also produced multiple output structures for each input structure by generating different protonation states, stereochemistry, tautomers, and ring conformations.
Ligands were kept in their first probable protonation and tautomerization state (lowest potential energy).
5) Choice of protein structure in which each molecule was docked
After ligand and protein preparation, the next step was to dock the ligand in FXR. Given the abundance of crystal structure data, the best crystal structure for each D3R ligand had to be chosen for the docking. For this purpose, we followed different steps:
a) SHAPE screening
First, we used the SHAPE module of the PHASE tool of Schrodinger in order to identify molecules that have similar shape to a shape query, which is either a template molecule or a set of spheres of defined radius and coordinates. As a template molecule we used the native ligands from all 28 downloaded crystal structures that we mention above.
For the SHAPE screening, we used the method “Volume scoring” together with the QSAR option.
“Volume Scoring” chooses the types of atoms between which volume overlaps will be calculated. Volume overlaps are computed between atoms that have the same atom type as defined by the choice of atom typing scheme below. The QSAR option calculates the volume overlaps between atoms that have the same pharmacophore type (Acceptor, Donor, etc.) as defined for Phase QSAR models.
We then grouped the D3R ligands based on the similarity that they had with the native ligands.
b) cross-docking in all crystal structures
Glide SP was then used to dock all ligands in all 28 crystal structures. We compared the docked posed with the native ligand poses to find the D3R ligands which had the best alignment with the native ligand.
c) interaction fingerprints
Apart from the alignment of each D3R ligand with the native ligands, we also calculated the interactions of the D3R ligands with the proteins and compared them to the interaction of the native ligands with the corresponding proteins. These interactions, are termed “Structural interaction fingerprints” in Maestro and can be used on a set of ligands that are docked to a receptor to calculate ligand similarities, and cluster or filter the ligands. The fingerprints record information on proximity of the ligand to receptor residues and the types of residues with which the ligand interacts. See Deng, Z.; Chuaqui, C.; Singh, J. J. Med. Chem. 2004, 47, 337-44 and Singh, J.; Deng, Z.; Narale, G.; Chuaqui, C. Chem. Biol. Drug. Des. 2006, 67, 5-12 for a detailed description of the method and its uses.
To use this facility, choose Scripts → Cheminformatics → Interaction Fingerprints in Maestro. More information on how to use it is available in the panel and the help topic.
Based on the best match between a) similarity of the D3R ligands interaction fingerprints with the native ligands interaction fingerprints, b) the SHAPE similarity, and c) the docking alignment with the native ligand, the final structure for docking of each D3R ligand was chosen.
d) choice of crystal structure to dock benzimidazoles
For benzimidazoles we did not follow the procedure above. Instead, we categorized the D3R benzimidazoles based on their similarity with the structure of the native ligand in each benzimidazole crystal structure, which was available in the PDB.
- If a benzimidazole had an ortho substituted benzene ring, then structure 3OLF was used.
Thus, FXR_14, FXR_24, FXR_25, FXR_27, FXR_28 were docked in 3OLF
- If a benzimidazole had an non-ortho substituted benzene ring, then structure 3OOF was used.
Thus, FXR_21, FXR_29, FXR_36 were docked in 3OOF
- If a benzimidazole had an saturated ring, then structure 3OKI was used.
Thus, FXR_6, FXR_7, FXR_8, FXR_9, FXR_13, FXR_19, FXR_20, FXR_22, FXR_26, FXR_30, FXR_31, FXR_32, FXR_35 were docked in 3OKI
5) Docking procedure
The docking took place with the program Glide of Schrodinger and the SP scoring function.
Glide searches for favorable interactions between one or more ligand molecules and a receptor molecule. The combination of position and orientation of a ligand relative to the receptor, along with its conformation in flexible docking, is referred to as a ligand pose.
The shape and properties of the receptor are represented on a grid by several different sets of fields that provide progressively more accurate scoring of the ligand poses. We thus started by creating the receptor grid for each of the 28 crystal structures. The hydroxyl groups in residues such as Ser, Thr, and Tyr and the thiol group in Cys can adopt different orientations with different ligands. Glide can allow such groups to adopt different orientations when ligands are docked, to produce the most favorable interaction. For Ser and Thr, the hydroxyls can be oriented in any of the three local minima, and likewise for the thiol of Cys; for Tyr, they can be in either of the two local minima. Using the Rotatable Groups tab we chose the hydroxyl of Ser, The, and Tyr residues within a distance of 4 Angstrom from the ligand to be treated flexibly.
The ligand poses that Glide generates pass through a series of hierarchical filters that evaluate the ligand’s interaction with the receptor. The initial filters test the spatial fit of the ligand to the defined active site, and examine the complementarity of ligand-receptor interactions using a grid-based method patterned after the empirical ChemScore function.
Poses that pass these initial screens enter the final stage of the algorithm, which involves evaluation and minimization of a grid approximation to the OPLS nonbonded ligand-receptor interaction energy.
Final scoring is then carried out on the energy-minimized poses. By default, Schrödinger’s proprietary GlideScore multi-ligand scoring function is used to score the poses. In this case, we used GlideScore SP as the scoring function, and a composite Emodel score was used to rank the poses of each ligand. Emodel combines GlideScore, the nonbonded interaction energy, and the excess internal energy of the generated ligand conformation.
Whenever a common core was was present with a native ligand, a core constraint docking was performed.
Exceptions from this procedure:
-For FXR_23 a DFT/B3LYP/6-31G* relaxed coordinate scan was performed for the dihedral around amide and the isoxazole using Jaguar (Schrodinger). This calcualation was performed in order to determine whether the plane between amide bond and isoxazole should be 180o. The calculation confirmed that the minimum energy conformation has the amide bond and the isoxazole moiety in the same plane. We thus then used only planar poses resulting from docking.
- For FXR_15, FXR_16, FXR_17, we performed a DFT/B3LYP/6-31G* Jaguar simulation to calculate the minimum energy pose of these ligands.
Glide SP Docking was then performed using as starting structure the Jaguar minimized structures.
Glide XP Docking was then performed using as starting structure the Glide SP produced pose.
6) Flexible ligand alignment
Whenever a common core was was present with a native ligand, the flexible ligand alignment was performed on the docked pose. The largest common scaffold was chosen or a common SMARTS was provided.
7) After the flexible ligand alignment, prime minimization of the complex generated by the flexible ligand alignment was performed to remove steric clashes. In all cases the steric clashes were removed.
8) Binding pose Metadynamics
In some cases, a native ligand did not have a common scaffold with the docked ligand. Therefore, for this case multiple docking poses were visualized. However, there was some uncertainty as to which docking pose should be chosen. Also, induced fit docking was performed together with docking for creating initial ligand conformations. To distinguish between these different poses, we performed binding pose metadynamics to choose the best viable pose.
This procedure was followed for ligands:
FXR_1, 4, 6, 18, 22, 23, 34, 74
Contributors: Christina Athanasiou, Sofia Vasilakaki, Zoe Cournia
Biomedical Research Foundation, Academy of Athens, Greece