The study was carried out on a Dell Workstation with a 2.26 GHz processor, 6 GB RAM, and 500 GB hard drive running in the Windows operating system. Bioinformatics software, such as AutoDock4.0 and online resources, were used in this study.

E6 protein of HPV-16 was selected as the drug target. The HPV-16 E6 (GenBank ID: AAD33252.1) protein sequence was retrieved from NCBI (http://www.ncbi.nlm.nih.gov/).

For modeling of the three-dimensional structure of the E6 protein, Phyre 2 server [20] was used. Structural refinement and energy minimization of the predicted model were carried out using the YASARA Energy Minimization Server [21]. The refined model reliability was assessed through Procheck [22], ProSA-web [23], ProQ [24], and Verify-3D [25]. The refined model was further verified by the ERRAT server [26].

Chemical structures of natural compounds (carrageenan, curcumin, EGCG, I3C, jaceosidin, WA), along with CAS (Chemical Abstracts Service) registry numbers, reported in the literature were retrieved from the PubChem database (Fig. 1) [27].

Protein-ligand docking studies were performed using the AutoDock 4.2 program [28]. It is one of the most widely used methods for protein-ligand docking. All the pre-processing steps for ligand and protein files were performed using AutoDock Tools 1.5.4 (ADT), which has been released as an extension suite to the Python Molecular Viewer [28]. ADT was used to prepare the receptor molecule (HPV-16 E6) by adding all hydrogen atoms into carbon atoms of the receptor, and Kollman charges were also assigned. For docked ligands, non-polar hydrogens were also added. Gasteiger charges were assigned and torsions degrees of freedom were allocated by ADT.

The Lamarckian genetic algorithm (LGA) was applied to model the interaction pattern between the E6 oncogenic protein and selected inhibitors. The grid maps representing the receptor proteins in the docking process were calculated using AutoGrid (part of the AutoDock package). A grid of 50, 50, and 50 points in the x, y, and z directions was centered on the p53-binding sites of E6 proteins. For all docking procedures, 10 independent genetic algorithm runs with a population size of 150 were considered for each molecule under study. A maximum number of 25 × 10⁵ energy evaluations, 27,000 maximum generations, a gene mutation rate of 0.02, and a crossover rate of 0.8 were used for the LGA. Autodock was run in order to prepare corresponding docking log file (DLG) files for further analysis.

The visualization of structure files was done using the graphical interface of the ADT tool and the PyMol molecular graphics system (http://www.pymol.org).

E6 protein of HPV type 16 has 158 amino acids in its protein sequence. Since the three-dimensional structure of E6 was not available, its structure model was predicted using the Phyre 2 server by taking chain C of the crystal structure of full-length HPV oncoprotein e62 in complex with the lxxll peptide of ubiquitin ligase e6ap (PDB ID: 4GIZ) as a structural template. The query coverage of the target-template alignment was 90% with 96% identity. The predicted structure was then subjected to the Yet Another Scientific Artificial Reality Application (YASARA) Energy Minimization Server for structural refinement. The total energy for the refined structure obtained from the YASARA Energy Minimization Server was -48,054.9 kcal/mol (score, -2.55), whereas prior to energy minimization, it was 81,463,636,897 kcal/mol (score, -4.21).

The stereochemistry of the refined model (Procheck analysis) (Fig. 2A) revealed that 99.2% of residues were situated in the most favorable region and additional allowed regions, and 0.8% (Gln106) of residues were in the generously allowed region, whereas none of the residues fell in the disallowed region of the Ramachandran plot (Fig. 2B). ProSA-web evaluation revealed a compatible Z score (Fig. 2C) value of -4.34, which is well within the range of native conformations of crystal structures [23]. The overall residue energies of the E6 3D model were largely negative except for one peak (Fig. 2D). The 3D model of HPV-16 E6 protein showed an LG score of 2.531 by Protein Quality Predictor (ProQ) tool, implying a high accuracy level of the predicted structure. A ProQ LG score >2.5 is necessary for suggesting that a model is of very good quality [24].

Similar assumptions were achieved using the ERRAT plot (Fig. 2E), where the overall quality factor was 91%. All of these outcomes recommended the reliability of the proposed model.

As all natural ligands (inhibitors) were found to be docked in various conformations and with varying binding energies, the lowest energy conformation was selected. Upon docking, the high-ranked binding energies of modeled structures of E6 oncoprotein with natural ligands were obtained (Table 1). HPV E6 degrades the p53 tumor suppressor protein. To accomplish the process of degradation, E6 proteins bind and modify the target specificity of the ubiquitin ligase E6AP. Based on the resolved structure of HPV-16 E6 by Zanier et al. [29], Cherry et al. [30] revealed in their docking analysis study that few flavonoid compounds bind in the hydrophobic pocket of the E6 and E6AP interface and mimic the leucines in the conserved α-helical motif of E6AP. Another study also showed that the biological activities of E6 are mediated by the interaction of E6 with the LXXLL motif of E6AP [31]. In our study, the E6 protein of HPV-16 showed one LXXLL motif (₁₀₃LCDLL₁₀₇). Out of six natural compounds included in our study, three (carrageenan, EGCG, and I3C) were found to interact with the LXXLL motif of HPV-16 E6 by forming a hydrogen bond with Leu106. According to Crook et al. [32], the p53-binding site is around amino acid residues 113-122 (CQKPLCPEEK) of HPV-16 E6 protein. The active site of the model was then analyzed based on the docking interaction of E6 with the p53-binding site with all 6 natural ligands.

Fig. 3 shows a docked pose of six ligand molecules with E6 protein, which showed almost similar docking results. The detailed interactions are given in Table 1. Each ligand interacts with the receptor at the p53-binding site. Among the 6 different ligands, carrageenan showed the lowest binding energy (-13.95 kcal/mol) and inhibition constant (0.059 µM) for the protein-ligand complex. Carrageenan was found to interact with the E6 residues Ala53, Phe54, Leu106, Ile108, Leu117, and Lys122 by forming hydrogen bonds. As mentioned in the introduction, carrageenan blocks HPV infection [10]; from our docking study, it might be revealed that the block of HPV infection might be due the high binding affinity of carrageenan towards the p53-binding site of E6 protein.

WA showed the second lowest binding energy of -7.58 kcal/mol and an inhibition constant of 2.77 µM. During docking with the receptor hydrogen bond, interactions were found with three residues (Ala53, Leu117, and Lys122) of E6.The other three compounds-EGCG, I3C and jaceosidin-were also found to interact with E6 protein with binding energies of -4.13, -4.06, and -4.01 kcal/mol, respectively. Leu106, Ile108, and Leu117 were the interacting residues of the receptor proteins with EGCG, whereas I3C interacted with Leu106, Pro116, and Lys122; Jaceosidin interacted with Phe54, Pro116, Leu117, and Lys122.

Curcumin was found to interact with E6 protein with less affinity (binding energy of -3.09 kcal/mol). All natural ligands were reported to inhibit HPV infection; our docking study also revealed the in silico validation of the inhibition. HPV-16 is known for causing high-risk cervical cancer. All ligands were found to interact with E6 oncoprotein of HPV-16 with significant binding energy and with amino acid residues known for p53 binding. This interaction might prevent E6 protein from interacting with host p53 protein, which may correlate with why these natural compounds are used to treat HPV infections.

Different natural products have been identified and used as promising drug sources against cancer caused by HPV. Advancement in bioinformatics and computational biology are proficient invalidating those natural drugs through in silico approaches. E6 proteins of high-risk HPV 16 inactivate p53 by inducing its degradation, a process that requires a ternary complex between E6, E6AP, and p53. Low-risk HPVs, such as HPV-6 and HPV-11, do not have the ability to degrade p53. Thus, E6 of HPV-16 is of considerable interest for designing novel inhibitors to overcome the challenges of cervical cancer. A high-quality 3D model of E6 obtained through a computational approach and docking analysis employing AutoDock 4.2 in this study provided high-throughput validation. This in silico approach may be of interest in designing new drugs from natural sources against cervical cancer.