Pharmacophore Modeling and Molecular Docking Studies of potential inhibitors to E6 PBM–PDZ from Human Papilloma Virus (HPV)

High-risk human papillomaviruses (HPVs) are known to cause cervical cancer. Vaccines are now available to prevent HPV infection. However, a clinically approved drug is yet not available to treat HPV. The PDZ(PSD−95/Dlg/ZO−1)−binding motif (PBM) in the E6 protein of HPVs targets the PDZ domain (known to be associated with oncogenesis) for degradation. Therefore, it is of interest to study PBM–PDZ interaction towards its possible inhibition with a potential inhibitor. Thus, four pharmocophore models of PBM−PDZ complex were developed. In order to obtain potent small molecules for its inhibition, a commercial compound database was screened using both these pharmacophore models and molecule docking method. These efforts identified four potential compounds (1−4) towards its inhibition with the docking scores range -18.2 to -15.0.


Background:
Human papilloma viruses (HPVs) belong to the papilloma virus family with over 170 members [1]. High-risk HPVs cause cancers of vulva, vagina, penis, oropharynx, anus, and are also considered to be the main causes of cervical carcinomas that is the second major cause of female cancer-related deaths worldwide [2]. Vaccines were developed and are currently used to prevent infection of HPVs in adolescent females. However, these vaccines are only effective against defined genomic types, which have been previously designed, and there is no expectation for the effectiveness of the vaccine in previously infected patients. Therefore, the development of a molecular drug targeting HPVs is necessary.
Genetically, HPV is a double-stranded DNA virus, which consists of a genome of approximately 8000 base pairs and at least six essential early-expressed proteins (E1, E2, E4-E7) and two essential late-expressed proteins (L1 and L2) [3]. The E6 protein has been found to be expressed in almost all HPVinfected cancer cells [4], and is thought to be one of the responsible factors of viral oncogenic effects and malignant transformation. In particular, in high-risk HPVs, the E6 protein binds to the tumor suppressor p53 via E6-associated protein (E6AP), which promotes the degradation of p53 [5]. However, immortalized epithelial cells are still detected in mutants without this interaction. Another contributing factor is the interaction of the PDZ(PSD-95/Dlg/ZO-1)-binding motif (PBM) with PDZ domains in the E6 protein (Figure 1a) [6]. Therefore, E6 PBM-PDZ binding is an attractive antiviral target for the development of chemical compounds.
In the current study, we created 4 semi-empirical pharmaco phore models of the E6C (the C-terminal of the E6 protein)-PDZ interaction, and screened a commercial database of approximately 4.5 million compounds using a pharmacophorebased molecular docking method. The results of the current study will offer guidance for further investigation of lowmolecule-weight HPV inhibitors.

Methodology: General
The pharmacophore and docking studies were performed on a PC running Windows using modules of the Molecular Operating Environment (MOE) software package. Circles represent the receptor-interacting parts of threonine and valine, which are underlined as key interacting conserved residues. Pharmacophores were created at these circles. Blue circles represent H-bond donors, pink circles represent H-bond acceptors, and the green circle represents the hydrophobic region.

Retrieval of target proteins
The X-ray structures of the PDZ domain and E6C fragment were obtained from Protein Data Bank (PDB, http:// www.pdb.org/) using PDB id 2I04.pdb [7]. Only monomer was used. After removing water molecules and hydrogen atoms, partial charges were added using the three-dimensional protonation module. Moreover, energy minimization was carried out using default parameters.

Creation of empirical pharmacophore models based on the structure and PDZ-E6C fragment interaction
First, MOE Ligand interaction module was used to calculate the ligand-receptor interaction. Based on previous reports showing that the X-S/T-X-V/I/L motif (Figure 1a) is critical and conserved in high-risk HPVs [8], according to the orientation of the threonine and valine residues of the E6C fragment in the cocrystal structure and protein residues, pharmacophore models with exclusion volumes were created using the Pharmacophore Query Editor.

Screening of the database based on pharmacophore models and molecular docking
A commercial database established by Namiki Shoji Co. Ltd., which comprises approximately 4.5 million compounds, was filtered using pharmacophore models. Subsequently, the MOE DOCK module was used, which contains steps for a conformation search of ligands, placement, scoring, refinement by energy minimization under a defined force field, and rescoring. Because this database is large and therefore the screening process is time-consuming, a two-step calculation was performed. First, docking was carried out without energy minimization calculation to obtain general information on whether a certain compound has the ability to bind to a specific site. In this step, the placement algorithm was set to Alpha Triangle, and the scoring function was set to London dG. Subsequently, a refined docking step with energy minimization calculation was carried out using only the top 30 poses of each molecule, under the force field MMFF94x. The same score function and other parameters were used as in the first step. The site was defined as the space of the ligand molecule (E6C) using the Site Finder module in both steps of docking. To confirm the parameters, so-called re-docking trials were carried on 2I04.pdb and 2 other similar structures (2I0L.pdb, 2I0I.pdb); 2I04.pdb showed a root mean square deviation (RMSD) as low as 0.52, suggesting sufficient repeatability (data not shown).  1 and (b) pharmacophore Model 2, the features constructed by the methyl group in the valine residues of the E6C fragment were different. In the former, one pharmacophore for the methyl group was used, and in the latter, one pharmacophore at one carbon was used. In (c) pharmacophore Model 3 and (d) pharmacophore Model 4, the white circles contain not only the features of the methyl group but also those of the carboxyl group.

Results & Discussion:
Generation of the empirical pharmacophore model Before creating pharmacophores, the interaction between PDZ and the E6C fragment within 2I04.pdb was checked using MOE Ligand interaction module (Figure 1b). We highlighted these important interactions and created pharmacophores adjacent to the interacting atoms/residues of the E6C fragment. For the valine residue, two pharmacophores at either the carbon atom of the terminal methyl group or only one large pharmacophore can be considered; therefore, a total of 6 or 7 pharmacophores were created. In addition, to exclude compounds with region(s) overlapping receptor atoms, 37 exclusion volumes were added (represented as yellow balls in Figure 2). Since it is difficult to meet all of the pharmacophores and exclusion volumes, selection and combination of these features were needed. Because valine occupied a cavity formed on the interacting surface of PDZ, whose constituent residues are conserved or substituted only by the hydrophobic leucine or isoleucine residue in high-risk HPVs, the hydrophobic interaction was thought to be important. The pharmacophore(s) shown in black are those defined as essential feature(s) in Models 1 and 2 ( Figure 2). Moreover, because of the size of the cavity, the carboxyl group of the valine residue can also be considered important. Models 3 and 4 included the pharmacophores of the carboxyl groups as essential (Figure 2).

Pharmacophore search
Using the 4 pharmacophore models, the NAMIKI database was searched, and the results are shown in Table 1 (see supplementary material). The number of features indicates the minimum number of pharmacophores that the compound's conformation matches with. The number of compounds obtained from the pharmacophore search based on Models 3 and 4 was lower than that based on Models 1 and 2. This may be simply due to the fact that the number of pharmacophores increased. In general, compared to the compounds obtained from Model 2, those obtained from Model 1 were more planar (data not shown). This is likely due to the placement of the valine pharmacophores. Next, we evaluated the similarity of the compounds obtained from Model 1 with those obtained from Model 2. SMILES files of all compounds obtained from Model 1 satisfying at least 4 features (including the essential features) and those obtained from model 2 satisfying at least 5 features (including the essential features) were compared. A total of 3147 compounds were found to be common between the two models. The same evaluation was carried out between Model 3, satisfying at least 3 features (including the essential features), and Model 4, satisfying at least 4 features (including the essential features), and 1550 compounds were found to be common. Although the structures were quite different, many of the hit compounds shared a structure of hydrophobic ring(s) oriented toward the (we refer to these compounds as Group 1), and those obtained from Model 3 (number of features = 3)/Model 4 (number of features = 4) was decreased to 4704 (referred to as Group 2). Before performing the docking studies, the database was decreased to 1/250, which substantially saved time.

Molecular docking
Molecule docking was carried out for each molecule, which provided a score according to the score function London dG. The lower the score, the more potent the predicted interaction. When the criterion was set to a score lower than -10, 707 compounds remained in Group 1 and 306 compounds remained in Group 2.
The top 4 compounds with the lowest docking scores (range -18.2 to -15.0) among compounds from Model 1 or 2 are shown in Figure 3. Compounds 1 and 2 were large, exhibiting a large number of carbon atoms and few hydrophilic groups; both compounds completely filled up the hydrophobic pocket on the surface of PDZ. Compounds 3 and 4 interact with 2 residues of PDZ, Gly463 and Phe464, by hydrophobic bonds (Figure 4). These 2 residues also interact with the valine of E6C (Figure  1b), suggesting potential to inhibit the PBM-PDZ interaction. On the other hand, although the best compound from Model 3 or Model 4 showed a dock score of -17.4 (data not shown), no other compounds showed scores as low as compounds 1-4.

Conclusion:
In order to reduce cervical carcinomas and other HPVs related diseases, there is a need to control HPV infection using a potential drug. However, to date no such drug has been approved. Therefore, it is of interest to study E6 PBM-PDZ interaction towards its inhibition using potential lead compounds. Thus, we report 4 pharmacophore models of E6C for small molecule screening. The top 4 potential compounds with features to inhibit this interaction were described. Furthermore, compared to compounds 1 and 2, compounds 3 and 4 are likely to be more active in biological assays and more suitable for further consideration, due to the fact that these compounds interact with residues Gly463 and Phe464 in PDZ. These two residues also interact with PBM in solution structure. Information on the specific mode of interaction will provide insight for the design of anti-HPV drugs.