Lead identification and optimization of novel collagenase inhibitors; pharmacophore and structure based studies

In this study, chemical feature based pharmacophore models of MMP-1, MMP-8 and MMP-13 inhibitors have been developed with the aid of HypoGen module within Catalyst program package. In MMP-1 and MMP-13, all the compounds in the training set mapped HBA and RA, while in MMP-8, the training set mapped HBA and HY. These features revealed responsibility for the high molecular bioactivity, and this is further used as a three dimensional query to screen the knowledge based designed molecules. These pharmacophore models for collagenases picked up some potent and novel inhibitors. Subsequently, docking studies were performed for the potent molecules and novel hits were suggested for further studies based on the docking score and active site interactions in MMP-1, MMP-8 and MMP-13.

Among MMPs, collagenases are intimately involved in collagen homeostasis by post-translational proteolytic degradation. They principally comprise MMP-1 (collagenase-1), MMP-8 (collagenase-2) and MMP-13 (collagenase-3) [6]. Collagenases are the only endogenous enzymes that can readily cleave the triple helical domain of fibrillar collagens I, II and III. Collagen degradation is commenced by collagenases by making a sitespecific cleavage about three-quarter of the distance from Nterminus, followed by spontaneous collagen denaturation [7]. These interstitial collagenases degrade type I, II and III collagen in cartilage; this is a committed step in the development of rheumatoid arthritis as well as osteoarthritis and is revealed by elevated levels of these collagenases [8,9].
Collagenases show interesting differences in the crystal structures, despite being highly homologous to one another. X-ray analyses of the enzyme-inhibitor complex of collagenases suggested that the S1' subunit is a selectivity pocket for collagenase inhibitors [10][11][12][13]. The S1' subsite, also called the S1'-specificity pocket, is the most prominent pocket within the catalytic domain of collagenases. Differences in the relative size and shape of the S1' pockets in MMP-1, MMP-8 and MMP-13 suggest that this pocket is a critical determinant of MMP inhibitor selectivity [1]. The quite flexible loop forms a major portion of the S1' pocket and it undergoes a conformational change on inhibitor binding [14,15]. The loop is of the same length in MMP-8 and MMP-13 and two residues are shorter in MMP-1 [16]. A comparison of the available 3D structure of MMP-1, MMP-8 and MMP-13 shows the variability of amino acid residues in the S1' loop. This variability of the amino acid residues in the S1' loop causes difference in the shape of loops [16]. The structural features of these enzymes are most decisive in determining MMP substrate specificity and thereby inhibitor specificity which is enclosed within the catalytic domain [17]. Synthetic inhibitors specifically targeting MMP-1, MMP-8 and MMP-13 are unclear. Selectivity is more vital in minimizing the detrimental effects during long term medical treatment [18]. It has been reported that side effects were observed in the clinical studies of collagenase inhibitors, because they showed broadspectrum inhibition. Therefore, specific inhibition of MMP-1, MMP-8 and MMP-13 are considered to be an attractive target in drug discovery research [19,20].
In the present study, we have generated pharmacophore models using Catalyst [21, 22] software for a diverse set of collagenase inhibitors (MMP-1, MMP-8 and MMP-13) with an aim to obtain pharmacophore model that would provide the chemical features responsible for activity. These pharmacophore features were used to screen the databases to find novel inhibitors. Further induced fit docking was performed to validate these inhibitors against MMP-1, MMP-8 and MMP-13. This in turn would be able to provide useful knowledge for developing specific new and active drug candidates targeting collagenases (MMP-1, MMP-8 and MMP-13).

Methodology:
Pharmacophore modeling using Catalyst: A set of 337 MMP-1 inhibitors with activity ranging from 0.4 nM to 100000 nM, 148 MMP-8 inhibitors with activity ranging from 0.13 nM to 78000 nM and 371 MMP-13 inhibitors with activity ranging from 0.16 nM to 100000 nM were selected from GOSTAR (gostardb.com). The molecules were divided into training and test set for the development and validation of pharmacophore models. The selection of training set is the most crucial part as it determines the quality of generated pharmacophore models. In this study, 21 of 337, 22 of 148 and 21 of 371 compounds were chosen for training set based on the diversity observed in their chemical structures and experimental activities for MMP-1, MMP-8 and MMP-13 respectively. The remaining compounds were used as test set for pharmacophore validation process. All the molecules were exported and then minimized using modified CHARMM force field in catalyst package [21, 22]. For each molecule, a maximum of 250 conformations were generated using the 'best quality' conformational search option within Catalyst's ConFirm module. It generates the conformations using 'Poling' algorithm. The molecules were then submitted to the catalyst hypothesis generation.

Model validation and Knowledge based screening
The best pharmacophore hypothesis was used initially to screen 316 MMP-1, 126 MMP-8 and 350 MMP-13 test set molecules. The same model has also been used to select potent molecules from 10,000 library molecules designed using Scaffold Hoping (Knowledge based screening). Library molecules were generated based on the knowledge of binding interaction of known ligands reported with MMP-1, MMP-8 & MMP-13 and also the common features necessary for biological activity of molecule [24][25][26]. These molecules were built using Cerius 2 software [27] and conformations for each compound were generated using best conformational analysis. These molecules were further screened for their activities using the developed pharmacophore models.
Ligprep produces a single, low energy, 3D structure for each input structure with various ring conformations, ionization states and tautomers using various criteria including molecular weight or specified numbers and types of functional groups present.

Protein preparation
Protein preparation and refinement studies were performed on MMP-1 (PDB ID: 1HFC), MMP-8 (PDB ID: 3DPE) and MMP-13 (PDB ID: 1XUC) using protein preparation module (Schrödinger suite, LLC) [28] in which the water molecules were removed, hydrogen atoms were added, bond orders were assigned and orientation of hydroxyl groups were optimized. Finally, energy minimization was carried out using default constraint of 0.3 Å RMSD and OPLS 2005 force field.

Induced fit docking
Induced fit docking method for protein structures of MMP-1, MMP-8 and MMP-13 was performed using Induced fit docking of Schrodinger package [28]. During docking process, the ligands were optimized using OPLS or MMFF force field, thus changing its conformation to find the best fit that can closely fit to the S1' pocket of MMP-1,MMP-8 and MMP-13. The binding affinity of each protein and ligand complex was reported as Glide Score [29]. All graphic images were picturised using PyMol program (www.pymol.org) [30]. Non-bonded interactions like hydrophobic was observed using LigPlot program [31] and these interactions can increase the binding affinity between target drug interfaces.

Result & Discussion:
Synthetic inhibitors taken for this study include hydroxamate, non-hydroxamate, carboxylate, phosphinate, aminocarboxylates, thiol, sulphonates, pyrrolidine, diazepine, etc., which tends to have a greater inhibition towards MMPs [18]. Hydroxamate diazepine and phenyl sulfonyl acetamide inhibitors are potent inhibitors of MMP-9 and MMP-13 both in vitro and in vivo, for osteoarthritis in rabbit model [32,33]. Since many of them have been quite potent and in some cases fairly selective against collagenases, these molecules were further taken as a basement to design the target specific inhibitors for collagenases (MMP-1, MMP-8 and MMP-13) using pharmacophore modeling.

Pharmacophore generation and validation studies using HypoGen
Ten hypotheses were generated using 21 diverse training set molecules for MMP-1 and MMP-13, and 22 molecules for MMP-8 in HypoGen. (Figure 1 a, b & c) show some of the molecules selected as the training set for MMP-1, MMP-8 and MMP-13. The best hypothesis for MMP-1 and MMP-13 consists of 1) one hydrogen bond acceptor, 2) one hydrogen bond donor and ring aromatic whereas MMP-8 consists of 1) two hydrogen bond acceptor, 2) one hydrogen bond donor and one hydrophobic.
The values of ten hypotheses such as cost, correlation (r), and root-mean-square deviations (RMSD) are statistically significant Table 1A B & 1C (see supplementary material). The pharmacophore (Hypo-1, 11 and 21 for MMP-1, MMP-8 and MMP-13 respectively) having high correlation coefficient (r), lowest total cost, and lower RMSD value was chosen to estimate the activity of test set. The best models Hypo-1,11 and 21 for MMP-1, MMP-8 and MMP-13 respectively has been given in (Figure 2) and the parameters that describe Hypo-1, 11 and 21 for MMP-1, MMP-8 and MMP-13 respectively are shown in Table 1A,   Two statistical methods were employed to rank the ten resultant hypotheses. In the first method, all the ten hypotheses were evaluated using a test set of known MMP-1, MMP-8 and MMP-13 inhibitors, which are not included in the training set. Predicted activities of the test set were calculated using all ten hypotheses and correlated with the experimental activities. Of the ten hypotheses, the best hypothesis Hypo-1 is characterized by the highest cost difference (58.88), lowest RMSD value error (0.87) and with correlation (0.89) for MMP-1. Hypo-11 is best for MMP-8, with highest cost difference (64.95), lowest RMSD value error (1.04) and correlation (0.85). Hypo-21 is the best hypothesis for MMP-13 with highest cost difference (63.90), lowest RMSD value error (1.15) and correlation (0.87). These results conclude that Hypo-1, 11 and 21 are best ranking pharmacophore for MMP-1, MMP-8 and MMP-13 respectively, among the 10 hypotheses obtained. In MMP-1 and MMP-13, all the compounds in the training set map hydrogen bond acceptor (HBA) and ring aromatic (RA), while in MMP-8, the training set map HBA and hydrophobic (HY) and these features revealed that they should be mainly responsible for the high molecular bioactivity. Thus this should be taken into account in discovering or designing novel inhibitors. The most active compounds 1, 22 and 44, has a highest fitness score of 6.80, 8.01 and 8.20 sequentially, when mapped Hypo-1, 11 and 21 to MMP-1, MMP-8 and MMP-13 respectively (Figure 3) whereas the least active compounds 21, 43 and 64 maps to a lowest value of 4.8, 4.13 and 5.73 (Figure 3). It is evident that as error, weight and configuration components are very low and not deterministic to the model, the total pharmacophore cost is also low and close to the fixed cost. Also, as total cost is less than the null cost, this model accounts for all the pharmacophore features and has a good predictability power. Analyzing the results, in MMP-1 out of 7 highly active molecules, 5 were predicted correctly as highly active, and the rest were predicted as moderately active. Among the 9 moderately active molecules, 6 molecules were predicted as moderately active, 2 were predicted as highly active and 1 was predicted as low active molecule. Out of 5 low active molecules, one was predicted as moderately active and remaining was predicted as low active. In MMP-8 out of 10 highly active molecules, 8 were predicted correctly as highly active, and the rest were predicted as moderately active. Among the 7 moderately active molecules, 2 molecules were predicted as highly active 4 were predicted as moderately active and 1 was predicted as low active molecule. Out of 5 low active molecules, 2 were predicted as moderately active and rest was predicted as low active. In MMP-13, out of 8 highly active molecules, 6 were predicted correctly as highly active, and the rest were predicted as moderately active. Among the 9 moderately active molecules, all molecules were predicted as moderately active. Out of 4 low active molecules, 1 was predicted as moderately active and rest was predicted as low active. Activities of the compounds were correctly predicted and fit values also confer a good measure of how well the pharmacophoric features of Hypo-1, Hypo-11 and Hypo-21 for MMP-1, MMP-8 and MMP-13 respectively were mapped onto the chemical features of the compounds. The best models Hypo-1,11 and 21 for MMP-1, MMP-8 and MMP-13 respectively has been given in (Figure 2) and the parameters that describe Hypo-1, 11 and 21 for MMP-1, MMP-8 and MMP-13 respectively are given in Table 1A,   Overall, a strong correlation was observed between the predicted Hypo-1, 11, 21 and the experimental activity for MMP-1, MMP-8, and MMP-13 inhibitory (IC 50 ) of the training and test compounds. However, Hypo-1, 11 and 21 models has a greater tendency to show false positives. This could be attributed to high structural similarity in the active and inactive MMP-1, MMP-8 and MMP-13 inhibitors, resulting in inability to discriminate this pattern by pharmacophore models. We have selected Hypo-1, 11 and 21 for MMP-1, MMP-8 and MMP -13 respectively as a 3D query to search a subset of knowledge based designed database of 10,000 compounds to retrieve compounds with novel structural scaffolds and desired features. The initial screening of Hypo-1, 11 and 21 yielded 3000 compounds and further cluster analysis of these hits corresponded to 220 unique cluster representatives. We further extended this study to structure-based design to limit the number of false positive hits and to further understand the binding of inhibitors to the active site of all three MMPs. (Figure 4 a, b & c) shows some of the identified and optimized potent lead molecules through virtual library screening for MMP-1, MMP-8 and MMP-13.

Docking studies for MMP-1, MMP-8 and MMP-13
For furhter validation, docking is performed for 220 hits (MMP-1, MMP-8 & MMP-13) using Induced fit docking mode of Maestro. Most of the compounds show hydrogen bonding interactions with the S1 pocket residues. Docking results shows that known (yellow green) and newly identified compounds (purple) occupy the S1' loop region with the almost similar conformations. Especially, highly active compounds are forming at least two specific or unique hydrogen bond interactions in the S1' loop and shows high glide score with collagenases (MMP-1, MMP-8 and MMP-13). The docking scores and the hydrogen bonding interaction for newly identified molecules with MMP-1, MMP-8 and MMP-13 were shown in Table 3 (see supplementary material). To investigate the ligand binding affinity at the hydrophobic S1' pocket, we further estimated the hydrophobic interactions using LIGPLOT software ( Figure 5). The hydrogen and hydrophobic interactions of highly active ligands for MMP-1, MMP-8 and MMP-13 are discussed below.

Binding mode of molecule -66 into the S1' loop of MMP-1
Docking results shows that both the ligands occupy the S1' loop region with the almost similar conformation with a glide score -7.816 and -7.089 Kcal/mol. LigPlot software was used to understand the in-depth interaction pattern between both the ligands and MMP-1. Hydrophobic interaction was identified with amino acid residues Leu 140, Phe 207, Arg 214, Ile 232, Gly 233, Leu 235, Tyr 237, Ser 239, Tyr 240, Phe 242, Ser 243, Asp 245 and Gln 247. Both the ligands form hydrogen bonding interaction with the Val 246 and specific amino acid residue Asp 245 in the S1' loop region of MMP-1 ( Figure 5).

Binding mode of molecule -69 into the S1' loop of MMP-8
Docking results shows that both the ligands occupy the S1' loop region with the almost similar conformation with a glide score -7.240 and -7.216 Kcal/mol. For the above two ligands hydrophobic interaction was seen with amino acid residues Leu 160, Leu193, Val 194, His 197, Glu 198, Ala 213, Leu 214, Pro 217, Asn 218, Tyr 219, Ala 220, Arg 222, Thr 224 and Ser 228 using LigPlot software. Phenyl group in both the ligands occupy the solvent exposed region including His 207, Ser 228 and Thr 224. Both the ligands form hydrogen bonding interaction with the Pro 217 which is unique among MMP-8 and also form pi-pi stacking with the Arg 222 in the S1' loop region of MMP-8 ( Figure 5).

Binding mode of molecule -72 into the S1' loop of MMP-13
Docking results shows that both the ligands occupy the S1' loop region with the almost similar conformation with a glide score -10.892 and -9.203 Kcal/mol. Using LigPlot software, hydrophobic interaction was seen with amino acid residues Leu 185,

Conclusion:
Statistically validated pharmacophore models were generated for collagenase inhibitors to locate the spatial arrangement of features, which are necessary for biological activity. Out of them, Hypo-1, 11 and 21 performed superior for MMP-1, MMP-8 and MMP-13 respectively and also showed excellent prediction for inhibitory activities of the test set compounds and well complemented with receptor active sites, compared to the other hypotheses. The best hypothesis for MMP-1 and MMP-13 consists of one hydrogen bond acceptor, one hydrogen bond donor and ring aromatic, whereas MMP-8 consists of two hydrogen bond acceptors, one hydrogen bond donor and one hydrophobic group. These pharmacophore models were used to retrieve the molecules from the databases which were further validated using docking studies. Finally, three structurally diverse compounds with high Glide scores and interactions with critical active site amino acids for MMP-1, MMP-8 and MMP-13 were identified. From this study, we suggest that these molecules can be used for further studies and also serve as potential leads against collagenase infections.      a + indicates that the predicted IC50 is higher than the experimental IC50; − indicates that the predicted IC50 is lower than the experimental IC50; a value of 1 indicates that the predicted IC50 is equal to the experimental IC50; b Fit value indicates how well the features in the pharmacophore overlap with the chemical features in the molecule. Fit = weight × [max (0.1−SSE)] where SSE = (D/T) 2,; D = displacement of the feature from the center of the location constraint and T = the radius of the location constraint sphere for the feature (tolerance). c Activity Scale -IC50 < 100 nM = +++ (highly active); IC50 =100 -1000 nM = ++ (moderately active ); IC50 > 1000 nM = + (low active ).