Pharmacophore based virtual screening for identification of effective inhibitors to compact HPV 16 E6 a triggered cervical cancer

doi:10.21203/rs.3.rs-1996595/v1

Download PDF

Research Article

Pharmacophore based virtual screening for identification of effective inhibitors to compact HPV 16 E6 a triggered cervical cancer

https://doi.org/10.21203/rs.3.rs-1996595/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Cervical cancer, one of the most common causes of cancer-related death among women in the world, has been linked to the presence of a particular oncoprotein that is predominantly transferred through sexual contact with an infected host. In 90% of cervical cancer deaths, a correlation has been found with the expression of the viral genome of HPV16 E6. As a result, HPV16 E6 has emerged as an optimistic therapeutic drug target for the treatment of cervical cancer. In order to develop a drug that is capable of disturbing the genome expression activity of HPV16 E6, it is imperative to identify the key chemical features of its known inhibitors. In this study, we present an investigationon identifying potential inhibitors of HPV16 E6 by utilizing pharmacophore-based virtual screening, molecular docking, ADME prediction, and molecular dynamics simulation. In the initial stage, we generated a ligand-based pharmacophore model based on the features of four known HPV16 E6 inhibitors (CA24, CA25, CA26, and CA27)via the PHASE module implanted in the Schrödinger suite. We constructedfour-point pharmacophore features, which consists of three hydrogen bond acceptors (A) and one aromatic ring (R). The common pharmacophore features further employed as a query for virtual screening against the ASINEX database via the Schrödinger suite. The pharmacophore-based virtual screening filtered out top 2000 hits, based on the fitness score. We applied the molecular docking studies for further compound filtration using Glide which provide three ligand filtering phases, namely HTVS, SP, XP. Initially, 1000 compounds were obtained from HTVS docking. Based on the glide score, they were further filtered to 500 hits by employing docking in standard precision mode. Finally, the best four hits were identified using docking in XP mode. The four compounds were then further subjected to ADME profile prediction by engaging the Qikprop module. The ADME properties of the four fell within a satisfactory range and the compounds exhibited anticipated pharmacokinetic properties. These compounds were further investigated to determine the binding stability of the protein-ligand complex at a different time scale (100 ns) by using the desmond package for a molecular dynamics simulation. These molecular dynamics simulation studies revealed that theCYS 51 and GLN 107 proteins are residues of HPV 16 E6 binding sites, and the root mean square deviation (RMSD) and root mean square fluctuations (RMSF) values for these residues were also found to be within satisfactory ranges and suggest a crucial role in enhancing the stability of the protein-ligand complex during the simulation. From these computational investigations, we concluded that the four potential compounds are appropriate for further study, and potential clinical investigation, as HPV16 E6 inhibitors.

Cervical cancer

HPV 16 E6

Pharmacophore modeling

Virtual screening

ADME

Molecular Dynamics Simulation

Behind breast cancer, colon cancer, and lung cancer, cervical cancer is one of the most common causes of cancer and the associated deaths are chiefly affecting women in developing countries [1]. On average, of the 493,000 new cases diagnosed throughout the world every year,, 233,000 of them result in fatality. The United States observed 11,070 cases of cervical cancer among women, and resulted in 4000 deaths. Human papilloma virus (HPV) is one of the precarious microorganisms mainly responsible for cervical cancer predominantly affecting women. Presently, more than 200 diverse classes of HPV have been identified and subsequently characterized into five genera: alpha (α), beta (β), gamma (γ), mu (µ), and nu (ν). Additionally, depending on the tissue that hosts the virus, HPV can be characterized into cutaneous and mucosal types [2, 3]. Cutaneous HPV types typically infect the basal epithelial cell of the hands and feet, resulting in plantar warts, common warts, and flat warts. Mucosal types infect the inner lining of tissues, such as those present in the respiratory tract, oropharyngeal region, and anogenital epithelium. Depending on the severity, the mucosal HPV infections may be designated as either high risk or low risk. High risk mucosal HPV types are frequently found in malignant tumors, such as those present in the cervical, head, and neck, while low risk types are mainly associated with benign lesions. Similarly, beta HPV is mainly responsible for skin wards, squamous cell carcinoma, and epidermodysplasia. The size of HPV genome ranges from around 6.8 kb to 8 kb in length. It possesses three core regions: The first is an early gene coding region E1, E2, E4, E5, E6 and E7, which are used for viral genome transcription, cell control, expression and transformation of HPV. The second is a late gene coding region containing major (L1) and minor (L2) capsid proteins. The third is long control region localized amongst open reading frames L1 and E6. 99.7% of the cervical cancer is directly linked to one or more oncogenic types of HPV, while low risk HPV, such as HPV6 and 11, are predominantly involved in the formation of benign genital warts and high risk HPV types, such as HPV16, 18, 31, 33, 45.They are mainly responsible for the risks associated with cervical cancer. Furthermore, 70% of all cervical cancer may be attributed to HPV16 (50%) and HPV18 (20%) [1, 3]. As a result, these forms of HPV have garnered more research interests. HPV16 E6 is an oncoprotein containing 158 amino acid residues and two zinc finger domains in which the conserved motifs are in all HPV types. The HPV16 E6 utilizes its conserved motif (LXXLL) to bind cellular proteins, namely E6AP, to form a complex of HPV16 E6/ E6AP, and then the complex binds to the central region of P⁵³ protein. Notably, these interactions can lead to degradation of the tumor suppressor p⁵³ via the ubiquitin pathway. When the viral particles of HPV enter body, they can alter the normal cell cycle and the apoptosis control, and then causes mutations which can lead to cancer. Overexpression of HPV16 E6 is often associated with squamous cellcarcinomas and other cancers including breast, vaginal, penile, vulvar, head, and neck [4, 5]. The E6 pocket protein is one of the most outstanding drug targets for the development of vaccines against HPV infection. So far, several research groups have attempted to investigate the destruction of the E6-E6AP interaction through different approaches. As evidenced by their collective works, RNA interference (RNAi), peptide aptamers, intracellular antibodies, and small molecules can lead to reinvigoration of p53 tumor suppressor activity, and eventually, the elicitation of cellular apoptosis in HPV-associated cervical cancer cells. Investigators have also attempted to find novel drugs through computer aided drug design [6, 7]. Currently, the Food and Drug Administration (FDA) has approved certain cervical cancer vaccines for prevention, control, and protection against HPV, such as Gardasil, Cervarix, and 9-valent. These products, due to their high costs, are often inaccessible for those infected with HPV and have also been linked to adverse side effects such as headaches, fevers, nausea, diarrhea, abdominal pain, fatigue, myalgia, gastrointestinal symptoms, arthralgia, and syncope [8]. There are currently no available effective inhibitors for HPV, providing an urgent need for research into methods providing protein-based inhibitors for HPV. In this work, we attempted to discover novel potent inhibitors for HPV16 E6 by using pharmacophore-based virtual screening from the ASINEX database. Initially, we constructed ligand-based pharmacophore models that were based on four known HPV16 E6 inhibitors. Then, the robust model was used to provide a query for database screening. The resulting hit compounds were selected according to fitness score. Drug-like properties were subsequently molecular docking, and evaluated for these compounds. Finally, the conformational stability of the protein-ligand complex was evaluated by Desmond.

Generation of common pharmacophore hypothesis

For this study, the pharmacophore models were constructed using PHASE (Schrodinger, New York, USA), which has comprehensively been used for pharmacophore-based studies which includes 3D database searching. Four well-known HPV16 E6 inhibitors (Figure 1) from recent literature were used for the generation of the pharmacophore models [9]. A Pharmacophore model is a set of chemical features that are essential for a small molecule to interact with a specific target protein and stimulate or inhibit their activity. All the 2D structure of the four known HPV16 E6 inhibitors were sketched and built using Chem Draw (Chem Draw Ultra 7.0.1(Cambridge Soft Corporation, Cambridge, MA, USA; http://www.cambridgesoft.com) [10]. After importing the constructed 2D structures into the Ligprep module that is incorporated in PHASE, hydrogens were attached onto the structures, then the 2D structures were converted into 3D, and then stereoisomers were considered. After, the most probable ionization state at user-define pH was determined. All of the structures were ionized at the neutral pH of 7. The energy minimization of the conformers were achieved using a rapid torsion angle approach provided by the MMFF’s implicit GB/SA solvent model. 1000 conformers were generated using a 100-step preprocess and a 50-step post process. The conformers were filtered through a relative energy difference of 50Kcal/mol^-1 and a maximum heavy atom deviation of 2Å, which resulted in a threshold relative of the lowest energy conformers. To obtain identical conformers, the pair distance of corresponding heavy atoms must be below 1.00 Å. They are applied only after the energy difference threshold of the two conformers are within 1 Kcal/mol^-1. The pharmacophore model was built by using the “Develop pharmacophore model” module of PHASE. PHASE offers a platform with six built-in features that include a hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic group (H), negatively ionizable (N), positively ionizable (P), and an aromatic ring (R). A four-point pharmacophore hypothesis with respect to all the conformation of the active molecules was developed, which has a set of spatial arrangement features including a minimum inter site distance of 2.0Å and a terminal pharmacophore box size of 1.0Å. The common pharmacophore hypotheses were scored by setting the root mean square deviation (RMSD) value below 1.2Å, and the vector score value to 0.5. The scoring algorithm included the contributions from the alignment of site points, vectors, selectivity volume overlap, number of ligands matched, and their relative conformational energy and activity [11]. The best pharmacophore hypothesis were selected based on survival score. Then, the robust hypothesis was subsequently applied as a query for database screening for identification of potential inhibitors.

Figure1: The 2D structure of four known HPV16 E6 inhibitors.

Pharmacophore based database screening

Pharmacophore based virtual screening is a rapid and accurate computational technique and plays a fundamental role in searching for hit candidates with desired pharmacokinetic properties. In this work, the most magnificent pharmacophore hypotheses (AAAR) were used as the 3D query for potent HPV16 E6 inhibitors from a large chemical library called ASINEX. The database screening was carried out using phase virtual screening protocol which was implemented in a Schrodinger suite. Primarily, the hit compounds were retrieved by applying fitness value, and then the retrieved compounds were further scrutinized using molecular docking [12].

Molecular Docking Studies

Molecular docking is a crucial tool for exploring the interactions between the target and small molecules. In this study, we performed the docking studies using the Glide version 5.5. Before the docking study, the protein and ligands were prepared to minimize the hindrances while docking. Initially, the 3D coordinates of the target structure, human Papilloma virus HPV16 E6 (PDB ID: 4GIZ), was retrieved from the protein databank [13]. Before importing into Maestro, bond orders were assigned, unwanted water molecules were eliminated, formal charges, and hydrogen atoms were added. All atom force field (OPLS-AA) charges were assigned using the target structure and energy minimization provided from the preparation wizard panel of Maestro [14]. The structure was then subjected to a receptor gird generation panel. The binding pocket residue is crucial for structure-based docking. So far, no co-crystal ligand was found in HPV16 E6, so we collected the following residues: H78, R10, V31, Y32, L50, L117, P116, F54, C51, R55, V62, L67, Y70, I73, R102, K122, and R131 from literature [5, 7]. The grid box was generated based on the above mentioned residues as the centroid of the grid box (20Å×20Å×20Å) and utilizing the receptor grid generation panel implanted in the Schrodinger suite [15]. Finally, 3D structures of 2000 hit compounds were obtained from pharmacophore matching, and all the ligands were incorporated into the LigPrep module for ligand preparation. In this process, we utilized the following criteria: (i) the OPLS-2005 force field was used, (ii) all possible ionization states at pH 7.0±2.0 were generated using the inonizer, (iii) the desalt option was generated, (iv) tautomers were generated for all conformers, and (v) one low energy conformer was generated per ligand. Then, the prepared ligands were subsequently taken into the structure-based virtual screening process using Glide. Glide provides three special routes for the ligand filtration: (i) Glide/HTVS is achieved for swift screening, (ii) then, Glide SP mode were used to minimize the further screening process, iii) finally, we employed Glide XP mode for the selection of probable compounds based on the glide score and the glide energy [16, 17, 18, 19].

ADME properties prediction

Most of the drug discovery assignments dropped in later clinical trials were due to the disagreeable pharmacokinetic behaviors of drug candidates. To reduce the drastic side effects and increase the success rates of drug, the prediction of the ADME properties of drug candidates is crucial for the early stages of the drug discovery project. The experimental ADME prediction is expensive, time consuming, requires more manpower, and results in high animal wastages. A minimum of 2 to 3 weeks is needed to determine ADME of the drug candidates experimentally. However, in silico it takes only a few minutes to predict the ADME of thousands of compounds. Thus, the in silico calculation was the major tool for the prediction of the ADME properties of new chemical entities. A good drug candidate possesses some ideal characteristics, including: i) it can be efficiently absorbed from the site of administration, ii) it can be distributed sufficiently in the target site, iii) it escapes extensive metabolism in the liver and other extra hepatic tissues, and iv) it is not excreted too rapidly via the renal or the hepatobiliary clearance mechanisms [20, 21]. Qikprop is one of the most popular in silico tools available for predicting the ADME properties [22]. In this study, the most promising compounds were subjected to QikProp predictions to analyze the pharmacokinetics properties.

Molecular dynamics (MD) simulation

The study of the binding stability and intermolecular interactions of the target macromolecule with promising hits after docking process is a crucial step in the in silico drug discovery and development process. Here, we employed the most significant computational packages. Desmond [23] for the determination of the binding stability of protein-ligand complexes. The orthorhombic water box was used to create a 10 Å buffer region between the protein atoms and box sides. The systems were neutralized with Na⁺ and Cl^- ions. The OPLS-2005 force field was used for energy minimization. The temperature was maintained constant at 300 K, and a 2.0 fs value was obtained in the integration step [24, 25]. We executed MD simulations for the protein-ligand complexes for 100ns. Root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), and total energy of the complexes were analyzed by using event-analysis and simulation-interaction diagrams.

Selection of robust Pharmacophore hypothesis

For the ligand-based pharmacophore modeling, a set of four known HPV16 E6 inhibitors were selected. The pharmacophore hypothesis was constructed by employing the PHASE module of the Schrödinger Suite. To find the common pharmacophore hypothesis, energy minimized structures of the four known HPV16 E6 inhibitors were imported into PHASE and the compounds that are considered as active pharmacophore sites were generated for all compounds.Eventually, 45 hypotheses were engendered as HPV16 E6 inhibitors and categorized according to the survival score. The most robust pharmacophore model sorted out was associated with the four-point pharmacophore hypotheses containing three hydrogen bond acceptors (A), one aromatic ring (R), and the best pharmacophore hypothesis was used as a query for database screening. The distance between the pharmacophoric and their sites were shown in figure 2.

Pharmacophore based virtual screening

Pharmacophore based virtual screening and docking-based virtual screening is one of the most reliable, rapid, low cost, and efficient approaches for the identification of effective inhibitors from huge chemical libraries. Initially, the well validated pharmacophore model (AAAR) of HPV 16 E6 inhibitors were used as a 3D structural query for identification of potential compound with desired properties from the ASINEX database with fitness value ranges >1. Initially, 2000 hits were obtained as a base on fitness score. Then, the selected hit compounds were further taken into molecular docking studies. In this investigation, we used Glide docking for the identification of the hits. Glide docking protocol comprises of three filtering process steps: high throughput virtual screening (HTVS), standard precision (SP), and extra precision (XP). Initially, a total of 2000 compounds were docked into the active site of HPV16 E6. 1000 compounds were obtained with HTVS, then these ligands were applied in the SP mode, resulting roughly 500 compounds. Finally four potential ligands with binding free energy lower than -8.547kcal/mol were obtained from the Glide XP mode. Then, the four compounds were further taken into the ADME properties evaluation. All of the four hits possessed desired predicted ADME properties. The chemical names of the four lead molecules are shown in Figure 3.

Docking study

The docking poses of the top four lead molecules in the active site of the HPV16E6 are shown in Figure 4 and the XP results for the four lead molecules which are listed in Table 1.

Docking simulation of Lead 1 into HPV16 E6

Upon binding, four hydrogen bonding interactions were observed between Lead 1 and HPV16 E6. The first is a hydrogen atom of the Lead 1 with the backbone oxygen atom of the hydrophobic residue of CYS51 with a bond distance of 2.12Å. The second is an oxygen atom of the Lead 1 with the backbone hydrogen atom of the CYS51 and has a bond length of 2.16Å. The third is an oxygen atom of Lead 1 with a side chain hydrogen atom of the GLN107 with a bond length of 1.79Å. The final is a hydrogen atom of Lead 1 with the backbone oxygen atom of the polar residue SER74 and has a bond length of 2.41Å. Furthermore, the following residues, LEU 50, ILE 73, ILE 128, CYS 66, VAL 53, PHE 45, and LEU67 were involved in hydrophobic interactions with Lead 1. The glide score and glide energy were calculated to be -8.896 kcal/mol and -57.416 kcal/mol respectively.

Docking simulation of Lead 2 into HPV16 E6

Upon binding, three hydrogen bonding interactions were observed between Lead 2 and HPV16 E6. The first is a strong interaction between a hydrogen atom of the Lead 2 with the backbone oxygen atom of the hydrophobic residue of CYS51 and a bond distance of 2.05 Å. The second is a strong interaction between a hydrogen atom of Lead 2 with the backbone oxygen atom of the polar residue of SER71 and a bond length of 1.81 Å. The third is an oxygen atom of Lead 2 with the backbone hydrogen atom of the CYS51 and a bond length of 1.98Å. The following residues, TYR 70, LEU 50, ILE 104, CYS 66, and VAL 31 are involved in hydrophobic interactions with Lead 2. The glide score and glide energy were calculated to be -8.856 kcal/mol and -57.690 kcal/mol respectively.

Docking simulation of Lead 3 into HPV16 E6

Upon binding, three hydrogen bonding interactions were observed between Lead 3 and HPV16 E6. The first is a strong interaction between a hydrogen atom of the Lead 3 with the backbone oxygen atom of the hydrophobic residue of CYS51 and a bond distance of 2.11 Å. The second is a strong interaction between a hydrogen atom of the Lead 3 with the backbone oxygen atom of the polar residue of SER71 and a bond length of 1.70 Å. The third is an oxygen atom of Lead 3 with the backbone hydrogen atom of the CYS51 and a bond length of 1.99 Å. The following residues, TYR 70, LEU 50, ILE 104, CYS 66, and VAL 31 are involved in hydrophobic interactions with Lead 3. The glide score and glide energy were calculated to be -8.496 kcal/mol and -65.257 kcal/mol respectively.

Docking simulation of Lead4 into HPV16 E6

Upon binding, three hydrogen bonding interactions were observed between Lead 4 and HPV16 E6. The first is a strong interaction between a hydrogen atom of the Lead 4 with the backbone oxygen atom of the hydrophobic residue of CYS51 and a bond distance of 1.95 Å. The second is a strong interaction between an oxygen atom of the Lead 4 with the backbone hydrogen atom of the hydrophobic residue of CYS51 and a bond length of 2.17 Å. The third is a strong interaction between a hydrogen atom of the Lead 4 with the backbone oxygen atom of the GLN 107 and a bond length of 2.46 Å. The last is an oxygen atom of Lead 4 with the backbone hydrogen atom of the polar residue of SER 74 and a bond length of 1.98 Å. One π-π interaction was found between His 78 with a phenyl ring in lead 4. The following residues, TYR 70, LEU 50, ILE 104, CYS 66, and VAL 31 are involved in hydrophobic interactions with Lead 4. The glide score and glide energy were calculated to be -8.462 kcal/mol and -50.920 kcal/mol respectively.

ADME Properties prediction

ADME describes the destiny of drug-like compounds within living organisms, especially in the human system. The ADME properties of the four newly identified lead were assessed using QikProp. The above four lead molecules were found to have desirable drug-like properties based on Lipinski’s rule of five, which suggests a potential drug candidate have a molecular weight under 500kDa, contain less than 5 H-bond donors, have less than 10 H-bond acceptors, and a predictedlogP below 5. Then, these lead compounds were further appraised for their drug-like behavior through analysis of pharmacokinetic parameters required for absorption, distribution, metabolism, excretion and toxicity (ADMET) by using QikProp. The aqueous solubility (QPlogS) for the four lead molecules, a factor that is critical for estimation of absorption and distribution of drug within the body, ranges between -5.686 to -3.891. The cell permeability (QPPcaco2), a key factor governing drug metabolism and its access to biological membranes, ranges from~ 49 to ~1043. The predicted IC₅₀ values for the blockage of HERG K⁺channels are in the acceptable range, which is below -5. The predicted Brain/Blood barriers are under the acceptable range of ~ -2.461 to -0.644. The percentage of human oral absorption range from 62 to 100. All the pharmacokinetics parameters are fit well with the acceptable range defined for use of human. The results were listed in Table 2. According to this result, we suggest that these four compounds may be apt for human use.

Molecular dynamics simulation

The dynamic behavior of the interaction between the four Lead compounds with the HPV16 E6 was investigated using the molecular dynamics simulation package, Desmond, at different timescales under a physiological state.

MD analysis of HPV16 E6 - Lead 1

Fig. 5a shows the Root mean square fluctuations (RMSF) of the protein and Lead 1. The RMSD of the protein backbone of the HPV16 E6-Lead 1 complex showed a deviation from 2.90Å to 4.29Å during 1.26 ns to 36 ns. The protein backbone was stable from 36 ns to 77.40 ns. However, a large deviation was observed from 77.40 ns to 96.90 ns. The RMSD changed from 4.31Å to 5.78Å. The RMSD dropped back to 4.47Å at end of 100 ns. Likewise, the RMSD of the ligand was stable before 77.40 ns, but deviated significantly from 2.40 Å to 5.53 Å between 77.40 to 94.10 ns. The ligand RMSD finally dropped to 4.20 Å at the end of 100 ns. The fluctuations of each protein amino acid residue over the simulation period are displayed in Fig. 5b. The greater fluctuation indicates a decrease in the stability of the protein-ligand complex. Initially, a slight fluctuation was observed in PHE 2, ASP 56, and SER 143. The RMSF values were 6.32 Å, 2.75 Å, and 6.44Å, respectively. A large fluctuation in the side chain of SER 143 was also observed and the RMSF value was 6.39Å. The fluctuation of c-alpha showed a significant fluctuation in the GLN 3 with an RMSD of 6.25Å. The bar diagram, 5c, depicts the intermolecular interactions, including the hydrogen bonding, hydrophobic interactions, and water bridge interactions of the Lead 1 with HPV16 E6. The hydrogen bonding interaction was observed with CYS 51, SER 74, GLN 107, and has a maximum occupancy of 2.002%, 74.2%, 45.4% respectively. Hydrophobic interactions and water bridges were formed between Lead 1 with HPV16E6. Figure 5d shows that hydrogen bonding plays a crucial role in the binding stability of the HPV16 E6-Lead 1 complex. Key residues SER 74, CYS 51, GLN 107 were involved in the formation of hydrogen bonds for 63%, 99%, 100% and 88% with active site regions of HPV16 E6 during the simulation period.

MD analysis of HPV16 E6 - Lead 2

Figure 6a shows the RMSF of the protein and Lead 2. The RMSD of the protein backbone of the HPV16 E6 was steady at around 3.81 Å till the end of the 36 ns. Then, a large deviation spike was observed at 36.2 ns. The protein was stable during the remaining time and reached 2.31 Å at 99.80 ns. The RMSD of the ligand, on the other hand, showed a spike at around 60 ns. The ligand RMSD finally reached 2.75 Å at the end of 100 ns. Both protein backbones and ligands were more stable after 70 ns. The most fluctuating was found in the following amino acids: ASP 56 and SER 14., RMSF values ranged from 2.26 Å to 2.44 Å. Similarly, the following residues,GLN 3, and ASP 56, in the backbone, had more fluctuations. RMSF values were observed ranges from 2.70 Å to 2.54. Beside the residues, GLN 3, and ARG 55, in the side chain higher fluctuations and RMSF values ranged from 4.13 Å to 3.87 Å. The bar diagram of Figure 6c depicts the intermolecular interactions, which includes the hydrogen bonding, hydrophobic interactions, and water bridge interactions, of the Lead 2 with HPV16 E6. Hydrogen bonds were observed on TYR 32 and CYS 51, with a maximum occupancy of 19.8% and 35.3%, respectively. Hydrophobic interactions were observed on SER 70 with a maximum occupancy of 43.8%. Other hydrophobic interactions and water bridges were also formed between Lead 2 with HPV16E6. Figure 6d shows that CYS 52 and GLN 107 were involved in the formation of hydrogen bonds for 33% and 12%, respectively, with active site regions of HPV16 E6 during the simulation period. Other residues, such as TYR 32, PHE 69, and TYR 70 were involved in hydrophobic interactions with Lead 2. They elevated the binding stability of the protein-ligand complex during the simulation time.

MD analysis of HPV16 E6 - Lead 3

Figure 7a shows the RMSF of the protein and Lead 3. The RMSD of the protein backbone of HPV16 E6 were examined, and there were no changes in the first 5.10ns he in tprotein RMSD. The RMSD value was noted 2.66Å, but the slight fluctuation was observed from 32.20ns to 98ns and the RMSD values ranged from 3.26Å to 3.72Å, then the deviation was observed in 98.90ns and the RMSD value was noted as 4.17Å. Likewise, the ligand RMSD was observed and the RMSD value were noted in the ranges of 2.27Å to 4.17Å. The higher fluctuation of the ligand RMSD was monitored at 90.90ns and the RMSD reached 4.17Å. The most fluctuating was found in the following amino acid residue: SER 143, with an RMSF value of 5.31Å. Similarly, the backbone had more fluctuation in residue PRO 5 and SER 143, with RMSF values of 3.60Å and 5.70Å respectively. Beside the side chain of the protein themost fluctuating amino acid residues were found to be GLN 3 and SER 142, with RMSF values of 5.13Å and 5.83Å, respectively. Also, thekey residues PRO 5 and SER 143 of the c-alpha of the protein had more fluctuations were observed. The RMSF of the protein-ligand complex wasdepicted in figure 7b. The greater fluctuations indicates a decrease the stability of the protein-ligand complex. The bar diagram in Figure 7c depicts the intermolecular interactions, including the hydrogen bonding, hydrophobic interactions, and water bridge interactions, of the Lead 3 with HPV16 E6. Hydrogen bonds were observed on TYR 32 and CYS 51 with a maximum occupancy of 59.5% and 23.8% respectively. Various other inter-molecular interactions, such as hydrogen bonding, π-cation, water bridges, and hydrophobic interactions were formed between lead 3 with HPV16 E6. Figure 7d shows that CYS 51 and TYR 32 were involved in the formation of hydrogen bonds for 122% and 59% respectively, with the active site region of HPV16E6 during the simulation period. Also, one π-cation and water bridge interaction were observed on HIS 78 with a maximum occupancy of 33% and 47%, respectively. LUE 67 was predominantly involved in a hydrophobic interaction with Lead 3 and elevates the binding stability of the protein-ligand complex during the simulation time.

MD analysis of HPV16 E6-Lead 4

Figure 8a shows the RMSF of the protein and Lead 4. The RMSD of protein backbone was stable after 2.80 ns and the RMSD was noted at 2.47Å. Larger fluctuations were observed between 48 ns to 78.40 ns. The protein was more stable afterward and the RMSD reached 2.02Å. Likewise, the deviation was found in ligand RMSD between in the 48ns and 78.50ns, and the RMSD were noted ranges of 2.68Å to 2.87Å. Even though the ligand was more stable end of the simulation, the RMSD reached 2.08Å. Root mean square fluctuations (RMSF) displayed the fluctuations of each protein amino acid residue over the simulation time period (Fig. 8b).The most fluctuating amino acid residue was found to be in GLN 6, with an RMSF value of 2.15Å. Higher fluctuations indicate a decrease the interaction. Similarly, the backbone had more fluctuations in residues GLN 5 and GLN 91, with RMSF values of 2.06Å and 2.15Å, respectively. Besides the side chain of protein, the most fluctuating amino acid residues was found to be ARG 141, with an RMSF value of 6.77Å. The bar diagram Figure 8c depicts the intermolecular interactions, including the hydrogen bonding, hydrophobic interactions, and water bridge interactions, of the Lead 4 with HPV16 E6. Hydrogen bonds were observed on ASP 49and CYS 51. Among them, CYS 51 displayed hydrogen bond interactions with a maximum occupancy of 73%. Additionally, various inter-molecular interactions such as hydrogen bonding, hydrophobic interaction, as well as water bridges, formed between Lead 4 with HPV16 E6. Figure 8d shows that the hydrogen bonding, hydrophobic, water bridge was significantly responsible for the boost up of the biding stability of lead 4-HPV16E6 complex. We identified the following residues: LYS 11, ARG 10, ASP 49, and CYS 51 were involved in the formation of hydrogen bonds for 14%, 22%, 42%, and 50% with active site regions of HPV16E6 during the simulation period. Moreover, residues ARG 10, TYR 32 CYS 51, CYS 11, and ARG 102 were predominantly involved in the formation of the water bridge interactions with a maximum occupancy of 15%, 19%, 26%, 19%, 11% with Lead 4, respectively. Only one reside, LEU 67, was predominantly involved in hydrophobic interactions with Lead 4 and elevated the binding stability of the protein-ligand complex during the simulation time.

Cervical cancer is one ofmost common malignant cancer related deaths affecting women worldwide. The microorganism HPV16 E6 plays a crucial role on cervical cancer. In this study, we have identified potential HPV16 E6 inhibitors by using ligand-based pharmacophore modeling, virtual screening, ADME prediction, and molecular dynamics simulation using the Schrödinger suite. Initially, a four-point pharmacophore with three hydrogen bond acceptors (A) and one aromatic ring (R) as pharmacophore features based on the four known HPV16 E6 inhibitors. Then the best pharmacophore model was further employed as 3D search query to screen against 400,000 ASINEX database compounds. The hit compounds were concurrently subjected to filtering by checking desired features. All the 2000 hit compounds were further taken into molecular docking and ADME properties prediction for further refinement. Finally, we selected four compounds with adequate pharmacokinetic properties, and then the four hit compounds were subsequently taken into investigation on their binding stability of protein-ligand complexes. All the four complexes were stable through the whole simulation based on their RMSD. Their RMSF values and various intermolecular interactions were observed. The key residues CYS 51 and GLN 107 were the major participated residues that enhanced the stability of the complex through the entire simulation. The four new lead compounds obtained may be targeted for further experimentation leading to clinical trials.

Ethics approval and consent to participate

Not applicable. The ethical standards have been met

Consent for publication

Authors consent for publication

Availability of data and materials

From corresponding authors upon request

Competing interests

The authors declare that they have no competing interests

Funding

This compilation is a research article written by its author and required no substantial funding to stated.

Author Contributions

Conceived and designed the experiments: AM. Performed the experiments: AM. Analyzed the data: AM, GS, EF. Contributed software facility: AM. Wrote the paper: AM, H.-F.J, GS, EF. Assisted in writing the paper and referencing:H.-F. J, GS, EF.

Acknowledgements

The corresponding author would like to thank to Department of Biotechnology, Selvamm Arts and Science College (Autonomous), Namakkal for providing lab facilities as well as encouragement for this study.

Authors' information

¹Department of Biotechnology, Selvamm Arts and Science College (Autonomous), Namakkal, Tamil Nadu, India.

²Department of Chemistry, Drexel University, Philadelphia, PA 19104.

³Department of Biotechnology, Sengunthar Arts and Science College (Autonomous),Tiruchengode. Tamil Nadu,

India

Nicole S. L. Yeo-Teh, Yoshiaki Ito and Sudhakar Jha. (2018). High-Risk Human Papillomaviral Oncogenes E6 and E7 Target Key Cellular Pathways to Achieve Oncogenesis Int. J. Mol. Sci. 19(6) 1706.
VjekoslavTomaic. (2016). Functional Roles of E6 and E7 Oncoproteins in HPV-Induced Malignancies at Diverse Anatomical Sites. Cancers. 8 (10): 95.
Nicholas A. Wallace and Denise A Galloway. (2015). Novel Functions of the Human Papillomavirus E6 Oncoproteins. Annu. Rev. Virol. 2(1):403-423.
Joel Ricci-López, Abraham Vidal-Limon, MatíasZunñiga, Verónica A Jimènez, Joel B Alderete, Carlos A Brizuela, Sergio Aguila. (2019). Molecular modeling simulation studies reveal new potential inhibitors against HPV E6 protein. PLoS One.14(3): e0213028.
Srikanth Kolluru, Rosemary Momoh, Lydia Lin, Jayapal Reddy Mallareddy, John L Krstenansky . (2019). Identification of potential binding pocket on viral oncoprotein HPV16 E6: a promising anti-cancer target for small molecule drug discovery. BMC Molecular and Cell Biology. BMC Mol Cell Biol. 20(1):3s0.
Yves Nominé, Murielle Masson, Sebastian Charbonnier, Katia Zanier, TutikRistriani, François Deryckère, Annie-PauleSibler, Dominique Desplancq, Robert Andrew Atkinson, Etienne Weiss, Georges Orfanoudakis, Bruno Kieffer, Gilles Travé. (2006). Structural and functional analysis of E6 oncoprotein: insights in the molecular pathways of human papillomavirus-mediated pathogenesis. Mol Cell. 21(5):665-78.
Avinash Kumar, Ekta Rathi, Suvarna Kini. (2019). E-pharmacophore modelling, virtual screening, molecular dynamics simulations and in-silico ADME analysis for identification of potential E6 inhibitors against cervical cancer. Journal of Molecular Structure. Volume 1189, pp 299-306;
Jorge L Cervantes, Amy Hoanganh Doan. (2018). Discrepancies in the evaluation of the safety of the human papillomavirus vaccine. Mem Inst Oswaldo Cruz.113(8): e180063.
Cherry JJ, Rietz A, Malinkevich A, Liu Y, Xie M, Bartolowits M, Davisson VJ, Baleja JD, Androphy EJ. (2103). Structure based identification and characterization of flavonoids that disrupt human papillomavirus-16 E6 function. PLoS One. 8(12): e84506.
http://www.cambridgesoft.com.
Phase, version 3.6. Schrödinger, New York, NY, 2013.
Kh. Dhanachandra Singh, Muthusamy Karthikeyan, Palani Kirubakaran, SelvaramanNagamani. (2011). Pharmacophore filtering and 3D-QSAR in the discovery of new JAK2 inhibitors. Journal of Molecular Graphics and Modelling. 30: 186–197.
Zanier K, Charbonnier S, Sidi AO, McEwen AG, Ferrario MG, Poussin-Courmontagne P, Cura V, Brimer N, Babah KO, Ansari T, Muller I, Stote RH, Cavarelli J, Vande Pol S, Travé G. (2013). Structural basis for hijacking of cellular LxxLL motifs by papillomavirus E6 oncoproteins. Science. 2013 339(6120):694-8.
14. Protein preparation Wizard Maestro (2009) New York: Schrodinger LLC.
Maestro, version 9.1 (2009) Schrodinger, LLC, New York.
Lucía Pérez-Regidor, Malik Zarioh, Laura Ortega,and Sonsoles Martín-Santamaría. (2016). Virtual Screening Approaches towards the Discovery of Toll-Like Receptor Modulators. Int J Mol Sci. 17(9): 1508.
Ligprep, Version 2.3, Schrodinger, LLC, New York, NY,
Glide, Version 5.5 (2009) Schrodinger, LLC, New York, NY.
Dik-Lung Ma, Daniel Shiu-HinChan and Chung-Hang Leung. (2011). Molecular docking for virtual screening of natural product databases. Chem. Sci 2: 1656-1665.
Fidele NK. (2013). An insilico evaluation of the ADMET profile of the Streptome DB database. Springer plus. 2:353.
Saeed A. (2107). In silico ADME-Tox modeling: progress and prospects. Expert Opin Drug MetabToxicol. 13(11): 1147-1158.
QikProp, Version 2.3, LLC, New York, NY, 2010.
Ganesh Kumar V, Kranthi Raj K. Leela MC, Jayakumar SB, Venkateswara Rao T. (2105). High-throughput virtual screening with e-pharmacophore and molecular simulations study in the designing of pancreatic lipase inhibitors. Drug Des DevelTher. 9: 4397–4412.
Ramin ES. Ayhan U, Muhammet B. Mine Y. Mert M. Serdar D. (2017). Virtual Screening of Small Molecules Databases for Discovery of Novel PARP-1 Inhibitors: Combination of in silico and in vitro Studies. J Biomol Struct Dyn 35(9):1899-1915.
Subramaniyan V. Palani M. Srinivasan P. Ram A. Sanjeev kumar S. (2018). Novel ligand-based docking; molecular dynamic simulations; and absorption, distribution, metabolism, and excretion approach to analyzing potential acetylcholinesterase inhibitors for Alzheimer's disease. J Pharm Anal8(6): 413–420.

Table1. Glide extra-precision (XP) results for the four lead molecules, by use of Schrodinger 10.2.

S. No

Compound

Glide score (Kcal/mol)

Glide Energy

(Kcal/mol)

No. of hydrogen bonds

Interaction residues

Distance (Å)

Lead1

-8.896

-57.416

CYS 51 (2)

GLN 107

SER 74

2.12Å

2.16Å

1.79Å

2.41Å

Lead2

-8.856

-57.690

SER 71

CYS 51 (2)

1.81

2.05

1.98

Lead3

-8.496

-65.257

SER 71

CYS 51 (2)

1.70

2.11

1.99

Lead4

-8.462

-50.920

CYS 51 (2)

GLN 107

SER 74

1.95

2.17

2.46

1.98

Glide score (Kcal/mol), Glide energy (Kcal/mol), No. of hydrogen bond interaction, Interacting residues, and Distance between the protein and ligand (Å).

Table 2:Assessment of drug-like properties of the lead molecules as verified by QikProp (Schrodinger 10.2).

Compound	MW	HBD	HBA	QPLogPo/w	Caco-2	QPLogHERG	QPLogS	QPLogB/B	% Oral absorption
Lead 1	357.37	2.00	11.70	0.173	99	-5.372	-3.891	-1.748	64
Lead 2	429.43	4.00	11.60	0.820	49	-6.105	-4.640	-2.461	62
Lead 3	339.35	2.00	9.70	0.685	133	-5.384	-3.977	-1.455	69
Lead 4	470.51	1.00	8.00	4.513	1043	-6.624	-5.686	-0.644	100

Molecular weight (<500 Da); Hydrogen bond donors (<5); Hydrogen bond acceptors (<10); Predicted octanol/water partition co-efficient log p (acceptable range: 22.0 to 6.5); Predicted aqueous solubility; S in mol/L (acceptable range: 26.5 to 0.5); Predicted IC50 value for blockage of HERG K+ channels (acceptable range: below 25.0); Predicted Caco-2 cell permeability in nm/s (acceptable range: 25 is poor and 500 is great); QP log BB for brain/blood (-3.0 / 1.2); Predicted percentage of human oral absorption (25 % is poor).

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Pharmacophore based virtual screening for identification of effective inhibitors to compact HPV 16 E6 a triggered cervical cancer

Status:

Version 1

Abstract

Figures

Introduction

Materials And Methods

Results And Discussion

Conclusion

Declarations

References

Tables

Additional Declarations

Status:

Version 1