Considering Rotatability of Hydroxyl Groups for the Active Site Residues of MMP-13 in Retrospective Virtual Screening Campaigns

Structural Superimposition

The 35 halo x-ray structures of MMP-13 were retrieved form PDB. For metameric structures only the first model was kept. The water molecules were removed and all of the structures were superimposed to a reference structure. The RMSD histogram for all residues was graphed. The RMSD graph was employed for conducting flexible receptor docking. All of the above tasks were done in Chimera.

Molecular Docking

The structure 3I7I was selected for the present study which according to a previous report is one of the best available MMP-13 3D structures for performing virtual screening. The Protein was prepared by Protein Preparation Wizard. In summary, the water molecules were removed, metal binding states were generated, H-bond assignment and restrained minimization were done. A grid box was centered on the center of the co-crystalized ligand and the ligands length was defined as 23 A°. Small molecules library consists of actives and decoys were prepared based on the DUD-E [10] collection for MMP-13 target. Duplicates were removed and 541 actives as well as 34521 decoys were docked in the MMP-13 active site by Glide. For most of the docking studies the HTVS (high-throughput virtual screening) setting or SP (standard precision) were employed except for the final round of screening in which the XP (extra precision) setting was employed. Due to limited computational resource, for XP docking 50 actives and 2500 decoys were randomly selected from the original DUD-E data set for MMP-13. After each retrospective virtual screening the results were assessed statistically by ROC curve analysis and enrichment factor. Glide can consider different orientations for hydroxyl group of residues Ser, Thr and Tyr. For Ser and Thr the hydroxyl group can orient in three minima while for Tyr it can adopted two orientations. We set the hydroxyl group of Tyr²⁴¹, Thr²⁴², Tyr²⁴³ and Thr²⁴⁴ rotatable with different docking precision.

Statistical Analysis

Using receiver operating characteristic (ROC) curve and enrichment factor (EF) for determining the discriminating power of a method was explained elsewhere [11]. In brief, the ROC curve and EF were applied to determine the performance of each conducted virtual screening. The increase in area under the curve (AUC) of ROC curve can be used as an indicator of improvement in discrimination between actives and decoys. EF is shown by the ability of a particular method to retrieve actives with a high rank among decoys and defined as below:

EF = (actives _sampled / actives _total) × (N _total / N _sampled)

All of the statistical test and plotting were done using R (R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/.) Including packages: enrichvs and ROCR [12].

Flexible Residues Identified by Structural Superimposition

RMSD of the residues in the active site of different experimental 3D structures for MMP-13 was shown in Fig. (1). As it was shown that the high variation in residues positon were observed in Ile²⁴⁰–Leu²⁵¹. A previous study also reported the flexibility of Thr²⁴⁴-Gly²⁵⁰ residues (in particular Thr²⁴⁴) by investigating the two available NMR structures of MMP-13 [8]. The flexibility of Thr²⁴² was also suggested previously [13]. According to the Glide available method for considering flexibility and accessibility of residues for ligands, the flexibility of hydroxyl groups of Thr²⁴², Tyr²⁴³, Thr²⁴⁴ and Tyr²⁴⁵ were investigated.

Retrospective Virtual Screening

Table 1 summarizes the results of different virtual screening performed on MMP-13. For HTVS the EF1% for the virtual screening considered all four hydroxyl groups rotatable is higher than the one obtained for rigid docking. In docking with SP mode the EF1% of both docking which considered the hydroxyl groups of either Tyr²⁴¹ or Tyr²⁴¹-Thr²⁴⁴ rotatable are higher than EF1% of rigid docking. The difference between rigid and flexible docking is more obvious in SP docking compared to HTVS one. In XP docking the best performed flexible residues of previous runs were only selected. The results again showed an improvement in EF1% by considering the flexibility for all four residues.

Table 2 shows another finding of our study. Glide score takes into account a number of parameters. Glide lipo parameter is one of them that rewards favorable hydrophobic interactions. In all of the HTVS or SP virtual screenings the Gilde lipo scoring gave higher AUC accompanying with higher EF1% compared to Glide score. In the XP mode this was not seen that it could be due to the small sample size (~2500) were used for XP docking in comparison to HTVS and SP modes (~30000). The ROC curves were shown in Fig. (2).

Inspection of the Binding Poses

By comparison of ranked ligand list it was found that some actives got better score in XP run with flexible residues compared to rigid XP run. As an example the active compound CHEMBL10335 was selected for furtherer examination. This was ranked 15 with docking score -9.861 in the rigid XP docking while ranked 3 with docking score of -12.172 in the flexible XP dock. As it is shown in Fig. (3) the hydroxyl groups of both Thr²⁴² and Thr²⁴⁴ adopted different orientations in flexible docking. However, as the residue Thr²⁴² was located near the hydrophobic part of the ligand it seems that the closer hydroxyl group to this hydrophobic part in rigid docking at least in part, was responsible for the higher energy of the docked pose for ligand CHEMBL10335. Furthermore, the hydroxyl group of the residue Tyr²⁴¹ didn’t change upon binding of ligand CHEMBL10335 and this was seen in the binding of the most of the actives that rendered it as a lesser important residue in the bonding of the ligands.