Prediction of binding poses to FXR using multi-targeted docking combined with molecular dynamics and enhanced sampling

Advanced molecular docking methods often aim at capturing the flexibility of the protein upon binding to the ligand. In this study, we investigate whether instead a simple rigid docking method can be applied, if combined with multiple target structures to model the backbone flexibility and molecular dynamics simulations to model the sidechain and ligand flexibility. The methods are tested for the binding of 35 ligands to FXR as part of the first stage of the Drug Design Data Resource (D3R) Grand Challenge 2 blind challenge. The results show that the multiple-target docking protocol performs surprisingly well, with correct poses found for 21 of the ligands. MD simulations started on the docked structures are remarkably stable, but show almost no tendency of refining the structure closer to the experimentally found binding pose. Reconnaissance metadynamics enhances the exploration of new binding poses, but additional collective variables involving the protein are needed to exploit the full potential of the method. Electronic supplementary material The online version of this article (doi:10.1007/s10822-017-0074-x) contains supplementary material, which is available to authorized users.

List of ligands treated by RMD in the RMD submission to the challenge, and the number of clusters considered for each ligand. More specifically, the submitted poses were obtained in the following way: First, the MD and RMD trajectories were merged and clustered together. To confirm the stability of the resulting poses, the cluster centers of the top 4-5 clusters were used as starting points for additional MD simulations, each of length 20 ns (except the MD pose itself, for which the previous MD simulation was simply adopted). Finally, each MD simulation was clustered separately and the top cluster center was used as the submitted pose in the RMD submission. The ranking of the poses was determined by a manual procedure taking into account the size of the clusters in the combined clustering, the stability of the pose in the extra MD simulation, and the ad hoc decision not to rank the MD pose as the top pose.

Ligand
Number of clusters 22 4 27 5 32 5 Table S3 Amino-acid sequences used for the various simulations (table extending into the next four pages). The secret crystal structures had greater variation in the sequence, sometimes involving insertions and deletions that could not be easily restored without a⇤ecting the geometry. We opted for a compromise: Extra terminal residues (if present) were deleted but other insertions and deletions were kept as in the crystal structures. Point mutations anywhere in the structure were put back according to the apo-protein using the graphical user interface of UCSF Chimera and the swapaa command. The Dock Prep module of UCSF Chimera was used to keep only the highest occupancy for atoms with alternative positions, and was also used to replace each truncated side-chains with a complete side-chain of the same residue type [35]. The "gaps" in the backbone (poorly resolved residues that were missing in the PDB files) were filled by adding the missing residues from the apo-template and optimizing the added sequence using the ModLoop [36] web server, which optimizes loop conformations without relying on known protein structures. The final amino acid sequence of each system is given below. In analogy with the procedure for the old crystal structures, the protonation state of the new structures were set identical to the apo structure. The average RMSD and the number of poses with RMSD below 2Å are also reported. In cases where the shading status changes, an explaining comment is also given, in which "structure X" denotes the protein structure taken from the experimental complex with ligand X. It can be seen that a hypothetical procedure, which includes only protein structures with ligands of the same type, would have given only one additional correct top-prediction (for ligand 10, for which all the structures with better scores are from complexes of the benzimidazole type), i.e. in total 20 correct poses. It can also be noted that only 9 ligands had the best score for its "own" crystal structure (i.e. the "globally first" equal to the "self first" pose); this mainly reflects the similarity between the ligands.  Table S5 Reliability of the scoring function. The table shows the best score for any docking pose that has an RMSD below 2Å towards experiment, and the di⇤erence between this score and the "globally first" score given in Table 1. The ligands are sorted according to decreasing score di⇤erence, and the table only shows the 29 ligands for which correct poses were found. A di⇤erence of zero means that there is an acceptable pose that has the best score; this is not necessarily the "globally best" pose from

MM-PBSA analysis of MD trajectories
We selected some example ligands to investigate whether MM-PBSA estimation of binding free energies provides su⇥cient statistical significance to discriminate between various binding poses and, in that case, whether it ranks the poses in agreement with their RMSD towards experiment.

Methods
The MM-PBSA calculations were performed using the stand-alone g mmpbsa tool.The binding free energy was calculated for each system using 250 snapshots from a 30 ns MD simulation with the first 5 ns discarded, i.e. with a stride of 0.1 ns. The MM-PBSA energy was computed for each snapshot as a sum of three terms: Here, E int MM is the MM interaction energy (sum of electrostatics and van der Waals energy) between the protein and the ligand. G PB is the change in Poisson-Boltzmann solvation energy upon ligand binding, with a solute dielectric constant of 1, a water dielectric of 80, a solvent probe radius of 1.4Å, the smol model for construction of the cavity, Bondii radii, and otherwise default parameters of the g mmpbsa software. Finally, G np , the nonpolar solvation energy, was calculated using the following equation: where ⇥ = 5.43 kcal mol 1 nm 2 , ⇤ = 0.92 kcal/mol, and the SASA was calculated using a probe radius of 1.4Å and default parameters of the g mmpbsa software. Note that only relative binding free energies can be evaluated in this case owing to the absence of the entropy term from the calculation.

Error analysis
The values reported in Table S6 are the averages (µ) and standard deviations (⌅) over the N snapshots. The standard error of the mean (SE) is calculated by assuming normally distributed values while taking into account the statistical ine⇥ciency (s) of the samples, i.e. as The SE is an estimate of the statistical uncertainty and gives a lower bound to what energy di⌅erences between poses that are meaningful to discuss. More precisely, we apply Welch's unequal variances t-test to estimate the probability (p) for the mean values of two distributions being equal. Unfortunately, the systematic errors, primarily caused by the simplified treatment of solvation e⌅ects, are much more di⇥cult to estimate. In our analysis, we are not concerned with the errors in the absolute binding free energy, nor the errors in the relative binding free energy between various ligands. Instead, we are only concerned with the systematic shift between various binding poses for the same ligand. To exemplify that there still is a problem, we note that the three MD-simulated poses for ligand 30 have completely di⌅erent electrostatic contributions, ranging from 2 to 125 kcal/mol, and a similar range for the PB contributions, which are known to roughly counterbalance the electrostatic contribution, especially for charged ligands. Evidently, the electric field surrounding the ligand is very di⌅erent in the three poses, and thus one can not rely on cancellation of error between the various poses, with regards to the treatment of electrostatics and solvation. The real uncertainty is therefore much larger than the SE (which is only 1-4 kcal/mol), as the latter reflects the fluctuations around a single well-defined pose with no significant changes in the electrostatic surroundings.
Our interpretation of the results in Table S6 is that among the charged ligands (30 and 36), all comparisons except that between the "best" and experimental poses of ligand 36, probably have much larger uncertainty than indicated by the SE, possibly greater than 10 kcal/mol. For ligand 32, the experimental pose has rather di⌅erent electrostatics and thus the uncertainty in the di⌅erence towards the two other poses is larger than indicated by the SE, although probably not as large as for the charged ligands. For ligand 25, the uncertainty in the di⌅erence between the "first" and the other poses is also somewhat underestimated by the SE.

Results
A summary of the MM-PBSA results in relation to the RMSD from experiment are given in Table S7. For ligands 20, 24, 25, and 36, the ensemble started from the "best" pose is (in terms of RMSD) almost identical to that started from the experimental pose, and this is largely confirmed by the MM-PBSA results, with a maximum discrepancy of 1.4 kcal/mol (which is still somewhat larger than expected from the SEs).
For ligands 20, 24, and 25, the "first" pose has the least negative binding free energy, in agreement with its higher RMSD towards experiment, and contrasting to the prediction by the scoring function. The p-values are 0.0008, 0.0001, and 0.26, respectively, indicating that the energy difference relative to the "Exp" pose is significant in case of ligands 20 and 24, but not for ligand 25 (the apparent significancy relative to the "best" pose is not to be trusted, because the "best" and "Exp" poses are almost identical). For ligand 32, all three poses are di⌅erent, and again the Exp pose has the most negative binding free energy, in agreement with the crystal structure (p-values 0.0001 and 0.0001 for the first and best pose, respectively, relative to the Exp pose). The "first" and "best" poses are mis-ranked according to the RMSD towards experiment, but the energy di⌅erence is rather small and the real ranking might not be reflected by the RMSD anyway.
For the charged ligands, the energy ranking also agrees with the RMSD, both for the three di⌅erent poses of ligand 30, and for the two poses of ligand 36. However, although the energy di⌅erences are very large in magnitude, it should be noted that the uncertainties in these values are very di⇥cult to estimate owing to the electrostatic e⌅ects discussed above. Thus, we can not reliably classify the di⌅erences as significant.
In conclusion, this limited investigation suggests that the MM-PBSA method is useful for predicting the relative binding free energies between several poses of the same ligand, even though the actual binding free energies are much less accurate. The method is apparently more reliable for uncharged ligands, for which the systematic errors in the treatment of electrostatics tend to cancel out between the various poses. It is important to mention that standard deviations from ensembles generated by a single MD simulation are often underestimations, and thus more reliable results would have been obtained by collecting snapshots from multiple shorter MD simulations with di⌅erent initial velocities [51]. Table S6 Detailed MM-PBSA results for the 18 analyzed simulations. For each simulation, started either from the First, Best, or experimental (Exp) pose, we report the Van der Waals interaction energy (vdw), the electrostatic interaction energy (ele), the polar solvation energy (PB), the non-polar solvation energy (np), and the sum of these terms (Total). For each quantity, the mean (µ) and standard deviation ( ) over the ensemble is given, along with the standard error of the mean (SE) for the total energy. All values are in kcal/mol.    The scale on the y axis goes from 0 to 5Å. The scale on the x axis goes from 0 to 30 ns (many simulations were run until 50 ns but visual inspection ensured that no significant events occured in the 30-50 ns interval). Fig. S4 Hydrogen bond analysis for ligand 35 explaining its lack of conformational variation in the RMD simulations. A) Total number of protein-ligand hydrogen bonds during the RMD simulation. The running average is marked with orange. B) 3D representation of the ligand (cyan)-protein hydrogen bonds (red lines) generated using UCSF Chimera.

Fig. S5
Fluctuation of the ligand dihedrals in RMD simulations using only the ligand dihedrals as CVs (red), using the ligand dihedrals and side-chain dihedrals as CVs (blue), and using the ligand dihedrals and orientational CVs (cyan). The left, central, and right columns of panels show results for ligand 5, 13, and 15, respectively. Each panel shows how one of the ligand dihedrals varies with time. For ligand 5, which have only 5 rotable ligand dihedrals, the bottom three panels show the orientational CVs instead (which were only monitored for the blue and red simulations but included as CVs in the cyan simulation). The scale on the x axis always goes from 0 to 20 ns. The scale on the y axis always goes from 180 to +180 and is periodic so that the upper endpoint of the axis represents the same conformation as the lower endpoint of the axis.