Free energy perturbation calculations
SchrödingerRelease2016v2/Maestro10.6/Desmondv4.6/Glidev7.1/Phasev4.7/FEP+/LigPrepv38014/Epikv3.6
OPLS3 force field
Windows per transformation: 12
Assumed pH 7.4
Typical simulation time per window: 5 ns
Free Energy Perturbation Protocol
1) First, a crystal structure was chosen based on SHAPE screening of Schrodinger. 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 the available crystal structures and mapped them on the molecules from the spiros group and the sulfonamide group.
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.
Apart from the alignment of each D3R ligand with the native ligands, we also calculated the interactions of the spiros and sulfonamide 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.
The molecules were then visualized with the ones having the best interaction fingerprint first. From this ranking, we chose as the representative molecule the one with the pose closest to its minimum energy conformation.
A docking pose was selected for the reference ligands FXR_99, FXR_74 based on the “PosePredictionPRotocolTEmplate-ZC.txt”
2) For this docking pose, flexible ligand alignment of the spiros and sulfonamide molecules was performed using as a criterion the biggest common core among the two compound series.
3) Prime energy minimization followed for each ligand-receptor complex. Only the protein near the ligand was allowed to move (4 Angstroms). The ligand was not minimized. We chose this option because of the plasticity of the protein that we have observed when aligning the different protein structures.
4) These input poses were submitted to a FEP+/Desmond calculation. Desmond with FEP/REST provide an enhanced sampling technique for free energy perturbation, to calculate the relative binding affinities of the 26 sketched lead compounds binding to Arp2/3 complex. A total of 12 λ windows are used for both the normal FEP and FEP/REST calculations, but the FEP/REST simulation method modifies the potential energy for a localized region surrounding the binding pocket with a automated algorithm, which selects the REST region (“hot” region). The REST algorithm calculates the effective temperature of that region allowing the complex to sample the relevant structural rearrangements within the allowed sampling times. The final predictions are independent of the paths chosen to perform the transformation from one ligand to the other. As the lead compounds of each selected system are imported into the graphical interface, the FEP Mapper module of Desmond automatically generates the perturbation pathways by a variant of a mapping algorithm. Accordingly to their similarity scores the ligand pairs connect to each other by edges in a way that each edge represent one FEP/REST calculation. The minimal error in the free-energy calculation results is estimated based on the spatial conformations sampled in all simulations. As a measure of convergence of the calculations we ensured that all molecules run in closed thermodynamic cycles, the configured space of the lamda windows overlaps, and the cumulative ΔG vs time converges to one value. The systems ran on the graphics processing units (GPUs) cluster of the Greek national supercomputing facility, HPC ARIS.
5) FEP+ results gave us the final relative free energy differences ranking
6) We submitted the final snapshot of the MD trajectory for each ligand. All water molecules were removed.
7) The problem of handling charged mutations. The spiros group contained several ligands with carboxylic groups, which are deprotonated (negative charge) according to the LigPrep module of Schrodinger. However, the set contained also neutral molecules. The total charge of the system should remain conserved while performing a FEP mutation. Previous studies have shown that turning on or off a charge in a system during a FEP calculation is nontrivial when using PME for the treatment of electrostatics. [Kastenholz, M. A.; Hunenberger, P. H. J. Chem. Phys. 2006, 124 (12), 124106-124127. Kastenholz, M. A.; Hunenberger, P. H. J. Chem. Phys. 2006, 124 (22), 224501-224520]. The problem arises mainly from the fact that a charge change in PME is handled by introducing a uniform neutralizing background charge to enforce neutrality [Hummer, G., Pratt, L. R., and García, A. E. J. Phys. Chem. 1996, 100 (4), 1206-1215]. In FEP calculations this means that instead of calculating the free energy of turning off/on a charge, the calculated free energy is the sum of turning off/on the charge and the free energy of turning on/off a uniform neutralizing background charge, which might introduce an error in the overall difference in the free energy of binding.
Therefore, we chose the following strategy. Two sets of FEP mutations were performed. One contained only the set of deprotonated ligands and the other contained the same ligands in their protonated form, plus the neutral ligands. The calculations showed that the protonated and deprotonated ligands rank the same. The R^2 of two DDG values was 0.82, so the correlation of the values was very high. Therefore, we chose to keep the protected ligand results and include the neutral ligand ranking within the deprotonated set.
Contributors: Christina Athanasiou, Sofia Vasilakaki, Zoe Cournia
Biomedical Research Foundation, Academy of Athens, Greece
Physics-based pose predictions guided by native ligands
SchrödingerRelease2016v3/WaterMapv2.7/SchrödingerRelease2016v2/Maestrov10.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.
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
SP Docking
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:
1OT7
3OKI
3OMM
3OOF
3OLF
3OOK
1OSH
3FLI
3BEJ
3DCT
3DCU
3FXV
3GD2
3HC5
3HC6
1OSV
3L1B
3OKH
3P88
3RVF
3RUU
3RUT
4WVD
3OMK
3P89
4OIV
4QE8
4QE6
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.
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.
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.
Contributors: Christina Athanasiou, Sofia Vasilakaki, Zoe Cournia
Biomedical Research Foundation, Academy of Athens, Greece