Acessibilidade / Reportar erro

Determining the Molecular Interactions of Natural Inhibitors with the HPV-16 E6 Oncoprotein by Docking and Dynamics Simulation Studies

Abstract

The oncoprotein E6, a pivotal player in HPV-16-induced cancer, has long been the focus of extensive research. Building upon our previous study, we identified asarinin and thiazolo[3,2-a] benzimidazole-3(2H)-one-(2-fluorobenzylideno)-7,8-dimethyl (thiazolo) as potential potent anti-HPV-16 E6 oncoprotein agents. Utilizing the UniProt-derived E6 sequence, we employed I-TASSER to model the protein's three-dimensional structure. Subsequent molecular docking via AutoDock Vina, coupled with a 1,000 ps dynamic analysis under physiological parameters, revealed that both asarinin and thiazolo have a high chance of forming stable protein-ligand complexes with HPV-16 E6, displaying distinct molecular dynamic properties. Thiazolo exhibited superior stability in simulation, evident in ligand conformation and movement graphs, while asarinin excelled in terms of contact residues. Furthermore, SASA, hydrogen bond graphs, and the DCCM graph collectively underscore the comparable potential of both drugs as robust inhibitors of the HPV-16 E6 oncoprotein.

Keywords:
Anti-viral; Cervical cancer; E6 protein; HPV16; Molecular dynamic

HIGHLIGHTS

The E6 protein is a major oncoprotein linked to the development of cancer caused by HPV-16 infections.

Asarinin and thiazolo are expected to form stable protein-ligand complexes with the HPV-16 E6 protein.

Asarinin is more stable than the other compounds in terms of contact residues, while the ligand conformation and movement graph indicate that thiazolo is more stable in simulation.

The SASA, hydrogen bond, and DCCM graphs support both compounds being equally effective against HPV-16 E6 protein.

INTRODUCTION

Cervical cancer is the second most frequent type of cancer in Indonesia and the second most common type in females between 15 and 44 [11 Nurcahyanti ADR. Cervical cancer: The case in Indonesia and natural product based therapy. J. Sci. Med. 2016; 4: 1-7., 22 Domingo E, Noviani R, Noor M, Ngelangel C, Limpaphayom K, Thuan T, et al. Epidemiology and prevention of cervical cancer in Indonesia, Malaysia, the Philippines, Thailand and Vietnam. Vaccine. 2008; 26(12): 71-9.]. It is caused when keratinocytes are infected with high-risk HPV strains, most often HPV-16 [33 Graham SV. The human papillomavirus replication cycle, and its links to cancer progression: A comprehensive review. Clin. Sci. 2017; 131: 2201-21.

4 Tan S, de Vries EG, van der Zee AG, de Jong S. Anticancer drugs aimed at E6 and E7 activity in HPV-positive cervical cancer. Curr Cancer Drug Targets. 2012; 12: 170-84.
-55 Wang X, Huang X, Zhang Y. Involvement of human papillomaviruses in cervical cancer. Front Microbiol. 2018; 9: 1-14.]. Efforts to reduce high-risk HPV infections may take one of two forms: a preventative strategy via vaccination or a curative approach through medication. State-sponsored and commercial HPV vaccinations are expected to prevent up to 90% of HPV-associated malignancies but do not eradicate persistent HPV infections or affect cancer development. As a result, discovering new HPV treatment approaches is urgently needed, especially in developing countries [66 Huh WK, Joura EA, Giuliano AR, Iversen OE, de Andrade RP, Ault KA, et al. Final efficacy, immunogenicity, and safety analyses of a nine-valent human papillomavirus vaccine in women aged 16-26 years: A randomised, double-blind trial. Lancet. 2017; 390: 2143-59., 77 El-Zein M, Richardson L, Franco EL. Cervical cancer screening of HPV vaccinated populations: Cytology, molecular testing, both or none. J. Clin. Virol. 2016; 76: 62-8.].

The E6 protein is a crucial oncoprotein in the HPV-16 genome historically associated with cancer formation, which makes most research on HPV-related cancer concentrate on this specific protein [88 Pal A, Kundu R. Human papillomavirus E6 and E7: The cervical cancer hallmarks and targets for therapy. Front Microbiol. 2020; 10: 1-15.

9 Yeo-Teh NSL, Ito Y, Jha S. High-risk human papillomaviral oncogenes E6 and E7 target key cellular pathways to achieve oncogenesis. Int. J. Mol. Sci. 2018; 19: 1-27.
-1010 Doorbar J, Quint W, Banks L, Bravo IG, Stoler M, Broker TR, et al. The biology and life-cycle of human papillomaviruses. Vaccine. 2012; 30: 55-70.]. Through the use of its 26S proteasome, the HPV-16 E6 protein suppresses the transactivation and degradation of p53, resulting in DNA damage and immortalization of infected cells and inhibiting innate immune responses to infected cells through the IRF3 and TYK2 signaling pathways. The E6 protein interacts with IRF3 and suppresses its transcriptional activity, primarily by reducing the synthesis of IFN-β mRNA, which is needed to regulate innate immune responses [1111 Bernard X, Robinson P, Nominé Y, Masson M, Charbonnier S, Ramirez-Ramos JR, Deryckere F, Travé G, Orfanoudakis G. Proteasomal degradation of p53 by human papillomavirus E6 oncoprotein relies on the structural integrity of p53 core domain. PLoS ONE. 2011; 6: 1-10.

12 Tummers B, van der Burg S. High-risk human papillomavirus targets crossroads in immune signaling. Viruses. 2015; 7: 2485-506.
-1313 Reiser J, Hurst J, Voges M, Krauss P, Münch P, Iftner T, et al. High-risk human papillomaviruses repress constitutive kappa interferon transcription via E6 to prevent pathogen recognition receptor and antiviral-gene expression. J. Virol. 2011; 85: 11372-80.]. To carry out these tasks, the E6 protein has a unique binding groove that interacts with E6-associated protein (E6AP) or IRF3 to create protein dimers that are highly predicted to be the active oncoprotein conformation [1414 Shah M, Anwar MA, Park S, Jafri SS, Choi S. In silico mechanistic analysis of IRF3 inactivation and high-risk HPV E6 species-dependent drug response. Sci. Rep. 2015; 5: 1-14., 1515 Zanier K, Stutz C, Kintscher S, Reinz E, Sehr P, Bulkescher J, Hoppe-Seyler K, Travé G, Hoppe-Seyler F. The E6AP binding pocket of the HPV16 E6 oncoprotein provides a docking site for a small inhibitory peptide unrelated to E6AP, indicating druggability of E6. PLoS ONE. 2014; 9: 1-12.]. Numerous studies have shown that if E6 protein does not form a dimer with E6AP, its function as an oncoprotein is significantly compromised, most notably in the p53 degradation pathway [1616 Martinez-Zapien D, Ruiz FX, Poirson J, Mitschler A, Ramirez J, Forster A, et al. Structure of the E6/E6AP/p53 complex required for HPV-mediated degradation of p53. Nature. 2016; 529: 541-5.

17 Sailer C, Offensperger F, Julier A, Kammer K-M, Walker-Gray R, Gold MG, et al. Structural dynamics of the E6AP/UBE3A-E6-p53 enzyme-substrate complex. Nat. Commun. 2018; 9: 1-12.
-1818 Mortensen F, Schneider D, Barbic T, Sladewska-Marquardt A, Kühnle S, Marx A, Scheffner M. Role of ubiquitin and the HPV E6 oncoprotein in E6AP-mediated ubiquitination. PNAS. 2015; 112: 9872-7.]. As a result, we regarded the E6 protein and its complex as attractive therapeutic targets. Several anticancer agents, such as vorinostat, have been found to have E6-suppressing activity in vitro. Vorinostat, which is orally accessible, acts as a competitive inhibitor at the E6 binding groove occupied by p53 and suppresses E6 protein synthesis in HPV-positive cells. It reduced E6-mediated p53 degradation in H1299 cells with an IC50 value of 1.79±0.50 μM [1919 Messa L, Celegato M, Bertagnin C, Mercorelli B, Nannetti G, Palù G, Loregian A. A quantitative LumiFluo assay to test inhibitory compounds blocking p53 degradation induced by human papillomavirus oncoprotein E6 in living cells. Sci. Rep. 2018; 8: 1-11.]. Vorinostat at concentrations of 1 or 5 μM impairs E6's capacity to disrupt p53 in HPV-18 raft cultures. Also, 5 or 10 μM of vorinostat in the same cultures reduced the HPV-18 DNA copy number in infected cells by 94% and 98.7%, respectively [2020 Banerjee NS, Moore DW, Broker TR, Chow LT. Vorinostat, a pan-HDAC inhibitor, abrogates productive HPV-18 DNA amplification. PNAS. 2018; 115: 11138-47.].

In our previous study, we suggested asarinin from Zanthoxylum spp and thiazolo found in Myristica fragrans as potent compounds based on molecular docking study because they are predicted to occupy the same binding groove as E6AP and IRF3, acting as direct competitive inhibitors of E6AP and IRF3 in the forming of E6 protein complexes with better stability and affinity compared to vorinostat. In this study, we conduct further analysis using a few molecular dynamic parameters to determine further the potential of asarinin and thiazolo as E6 protein complex inhibitor drug candidates.

MATERIAL AND METHODS

Data retrieval and pre-docking screening

The HPV-16 E6 amino acid sequence was obtained from UniProt with ID P03126. The I-TASSER webserver was used to simulate the protein's three-dimensional structure (https://zhanglab.dcmb.med. umich.edu/I-TASSER/). The 3D model was picked based on the created model's rank, with the greatest C- and TM-score values. We tested 476 bioactive compounds derived from 18 different plant species, as well as vorinostat (CID: 5311) as a positive control. All compounds were evaluated using the Lipinski’s rule of five parameters (http://www.scfbio-iitd.res.in/software/drugdesign/lipinski.jsp), then minimized and converted to the AutoDock format using the PyRx application integrated with the OpenBabel GUI [2121 Jayaram B, Singh T, Mukherjee G, Mathur A, Shekhar S, Shekhar V. Sanjeevini: A freely accessible web-server for target directed lead molecule discovery. BMC Bioinformatics. 2012; 13: 1-13., 2222 Mirzaei H, Zarbafian S, Villar E, Mottarella S, Beglov D, Vajda S, et al. Energy minimization on manifolds for docking flexible molecules. J. Chem. Theory Comput. 2015; 11: 1063-76.].

Molecular docking process

AutoDock Vina is used for docking and is integrated with PyRx (https://pyrx.sourceforge.io/) [2323 Trott O, Olson AJ. AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J. Comput. Chem. 2009; 31: 455-61.]. We analyzed the target protein's whole structure. The molecular coverage area is 44.4289 × 38.7535 × 71.0904 Å, and the center coordinates are 63.4878 × 63.4683 × 56.2079 Å. The primary docking parameters are the molecule's affinity in kcal/mol, the location of the binding site, and the protein-ligand interaction [2424 Hidayatullah A, Putra WE, Sustiprijatno, Permatasari GW, Salma WO, Widiastuti D, et al. In silico targeting DENV2's prefusion envelope protein by several natural products' bioactive compounds. CMUJ. Nat. Sci. 2021; 20(4): 1-20.

25 Putra WE, Salma WO, Rifa’i M. Anti-inflammatory activity of Sambucus plant bioactive compounds against TNF-α and TRAIL as solution to overcome inflammation associated diseases: The insight from bioinformatics study. Nat. Prod. Sci. 2019; 25(3): 215-21.

26 Putra WE, Waffareta E, Ardiana O, Januarisasi ID, Soewondo A, Rifa'i M. Dexamethasone-administrated BALB/c mouse promotes proinflammatory cytokine expression and reduces CD4+CD25+ regulatory T cells population. Biosci. Research. 2017; 14(2): 201-13.
-2727 Putra WE. In silico study demonstrates multiple bioactive compounds of Sambucus plant promote death cell signaling pathway via Fas receptor. FTST. 2018; 3(2): 682-5.].

Molecular dynamics simulation

The top two ligands with the lowest binding affinity scores and a drug control from the prior project, asarinin, thiazolo, and vorinostat, respectively, were chosen for molecular dynamic simulations against the HPV-16 E6 protein. Under typical physiological conditions (37°C, 1 atm, pH 7.4, 0.9% salinity), the protein and ligand complex structures were prepared for 1,000 picosecond simulations. The molecular dynamics simulation was done with the MD-run macro program, and the molecular dynamics data were analyzed with the MD analyze and MD analyzer macro programs as our previous protocols in molecular dynamics simulation [2828 Heikal MF, Putra WE, Sustiprijatno, Rifa'i M, Hidayatullah A, Ningsih FN, et al. In silico screening and molecular dynamics simulation of potential anti-malarial agents from Zingiberaceae as potential Plasmodium falciparum lactate dehydrogenase (PfLDH) enzyme inhibitors. Trop. Life Sci. Res. 2023. 34(2): 1-20.

29 Hidayatullah A, Putra WE, Rifa’i M, Sustiprijatno, Widiastuti D, Heikal MF, et al. Molecular docking and dynamics simulation studies to predict multiple medicinal plants’ bioactive compounds interaction and its behavior on the surface of DENV-2 E protein. KIJOMS. 2022;8(3):531-42.

30 Hidayatullah A, Putra WE, Rifa'i M, Sustiprijatno, Heikal MF, Widiastuti D, et al. Molecular docking and molecular dynamic study of multiple medicinal plants’ bioactive compounds as human papillomavirus type 16 E5 protein inhibitor. Kuwait J. Sci. 2023. 50(1A): 1-14.

31 Hidayatullah A, Putra WE, Sustiprijatno, Rifa’i M, Widiastuti D, Heikal MF, et al. Concatenation of molecular docking and dynamics simulation of human papillomavirus type 16 E7 oncoprotein targeted ligands: In quest of cervical cancer's treatment. An. Acad. Bras. Cienc. 2023; 95(Suppl. 1): 1-16.
-3232 Hidayatullah A, Putra WE, Sustiprijatno, Widiastuti D, Salma WO, Heikal MF. Molecular docking and molecular dynamics simulation-based identification of natural inhibitors against druggable human papilloma virus type 16 target. Trends Sci. 2023; 20(4): 1-12.].

RESULTS

The outputs of 3D visualization and alignment using PyMol indicated that the most significant substantial variation in the structure of the protein-ligand combination between before and after molecular dynamics simulation occurred in the E6-vorinostat complex, with a root mean square deviation of 4.541 Å (Figure 1). The protein complex formed with asarinin and thiazolo exhibited a deviation value of less than 2.500 Å between pre- and post-simulation. Compared to the complex structure of E6-asarinin before the simulation, the complex structure post-simulation has a deviation value of 2.494 Å. Meanwhile, the complex structure of E6-thiazolo after simulation had the lowest deviation value of 2.181 Å in comparison to the complex structure before simulation.

Figure 1
The 3D and 2D structure visualization of binding HPV-16 E6 protein and (A) asarinin, (B) thiazolo, and (C) vorinostat. The green structures in 3D visualization indicate the protein-ligand complex before the MD simulation, while the cyan structures are after the MD simulation. The colored diagram in 2D indicates the protein-ligand complex after MD, while the grayscale residue in that diagram indicates the protein-ligand complex before MD.

The comparison of binding site molecular dynamics before and after the simulation, as illustrated in detail by the two-dimensional visualization results, reveals that the most significant residues and interaction changes occur in the vorinostat, representing only 60% of the initial residue that continues to interact after the simulation. Additionally, some residues interact between the first and last 12 residues. Further, the composition of the interaction between E6 and vorinostat altered. Previously, Gln114, Ile111, and Ser78 exhibited a combination of hydrophobic and hydrogen bond interactions. However, Gln114 and Ser78 lost their hydrogen bonds after the simulation, whereas Tyr39 gained a hydrogen bond in lieu of its hydrogen interaction (Table 1).

Figure 2
The ligand movement plot for asarinin (red line), thiazolo (green line), and vorinostat (black line) over 1,000 picoseconds of simulation.

Figure 3
A). The ligand conformation plot, B). The solvent accessible surface area (SASA) plot, C). The intramolecular hydrogen bonds, D). The intermolecular hydrogen bonds for asarinin (red line), thiazolo (green line), and vorinostat (black line) complexes over a 1,000 picosecond of simulation.

There was no change in the number of interacting residues in thiazolo, but post-molecular dynamic data revealed that five of eight residues were preserved throughout the molecular dynamic simulation process. Tyr77 also interacts with thiazolo through a moderate hydrogen bond (2.95 Å). Asarinin was designated the most consistent compound in terms of contact residues and type of interaction created since it exhibited no change in contact residues or type of interaction between before and after molecular dynamics simulation.

The simulation for 1,000 ps indicates that vorinostat's ligand movement was considerably more significant than the other ligands (Figure 2). Vorinostat exhibited a mean movement RMSD of 5.394±1.705 Å with a minimum and maximum of 0.687 and 7.739 Å, respectively, which is more than twice the trajectory observed for thiazolo (2.692±0.638 Å). The only ligand with a maximum movement of less than 4 Å (3.999 Å) is thiazolo. Meanwhile, asarinin, which had the highest affinity value in the previous research, had a ligand movement about 34% lower than vorinostat, with a mean value of 3.539±0.916 Å. Thiazolo is stable after 100 ps, asarinin's ligand movement is stable after 400 ps, and vorinostat is steady after 650 ps of simulation.

Table 1
Binding site residues comparison before and after molecular dynamics simulation

The ligand conformational RMSD values were within 3 Å for all tested ligands, including control (Figure 3A). However, vorinostat is recorded to have the most extensive conformational changes (mean: 2.231±0.409 Å), and vorinostat conformation appears stable after 200 ps. Asarinin and thiazolo consistently have maximum conformational changes below 2 Å. Asarinin was recorded to have a lower average conformational change (mean: 1.129±0.308 Å) than thiazolo (mean: 1.439±0.258 Å), but the RMSD chart shows that the conformation graph of thiazolo relative is more stable than asarinin. Solvent accessible surface area (SASA) simulations (Figure 3B) showed that the entire protein-ligand complex tested had fluctuating SASA values in the range between 10,500 and just below 11,500 Å2. The whole protein-ligand complex's SASA value tends to stabilize after passing 50 ps.

During the 1,000 ps simulation period, there was a decrease in intramolecular hydrogen bonds in the initial 50 ps of the simulation period by 38% (Figure 3C). Intramolecular hydrogen bonds fluctuated but remained stable between 100 to 120 until the simulation period ended. A similar result was also shown in the protein-solvent hydrogen bonds simulation (Figure 3D), which showed no significant fluctuations during the simulation period. The number of protein-solvent hydrogen bonds stabilized in around 350 bonds starting from 50 ps until the end of the simulation.

The DCCM visualizes (Figure 4) coordinated movements between protein residues by emphasizing highly correlated motions between specific residues in yellow and strongly anticorrelated motions in blue [3333 Valente RPDP, Souza RC, de Medeiros Muniz G, Ferreira JEV, de Miranda RM, E Lima AHL, et al. Using accelerated molecular dynamics simulation to elucidate the effects of the T198F mutation on the molecular flexibility of the West Nile virus envelope protein. Sci. Rep. 2020; 10: 1-6.]. The motion plots of all protein-ligand complexes are identical; only the degree of significance differs. Correlated mobility is observed in the majority of E6N and E6C α-helix and β-sheet structures. The E6-vorinostat complex appears to have a significantly higher degree of correlation and anticorrelation than the other two complexes. All α-helixes structures on E6N exhibit considerable correlation, whereas all α-helixes structures on E6C exhibit significant anticorrelation. The C-terminal IRF3 binding pocket and the C-terminal E6AP binding pocket also exhibited strongly correlated movement compared to the major alpha-helix and N-terminal of the two binding pockets. Similar movement patterns were seen in the E6-asarinin and E6-thiazolo matrices, but their intensities were much lower in both. In comparison to asarinin, thiazolo seems to have a lower movement intensity.

Figure 4
Dynamic cross-correlation matrix for each protein-ligand complex

DISCUSSION

The results of the three protein complexes' three-dimensional structure alignment revealed that all three complexes had structural deviations, most notably the E6-vorinostat complex, which had the most substantial deviation. However, visual inspection of the 3D visualization data revealed that all three compounds tested had moved and reoriented inside the cavity and had not escaped, with the vorinostat exhibiting the most significant deviation. The simulation results indicate that the vorinostat undergoes the most substantial motion and conformational changes. Meanwhile, thiazolo seemed to be more stable than asarinin throughout the simulation process, owing to the compound's ability to sustain a ligand movement RMSD of less than 3 Å. While asarinin may exhibit less conformational changes than thiazolo complex throughout the simulation period, thiazolo exhibits more stable conformation. Thus, during the simulation process, even though asarinin had the lowest conformational change, thiazolo was more stable when interacting with the E6 protein [3434 Aier I, Varadwaj PK, Raj U. Structural insights into conformational stability of both wild-type and mutant EZH2 receptor. Sci. Rep. 2016; 6: 1-10., 3535 Chopdar KS, Dash GC, Mohapatra PK, Nayak B, Raval MK. Monte-Carlo method-based QSAR model to discover phytochemical urease inhibitors using SMILES and GRAPH descriptors. J. Biomol. Struct. Dyn. 2021; 0: 1-10.].

Despite the deviation, asarinin retained residual contact with the E6AP/IRF3 binding pocket during the simulation, whereas thiazolo had some contact residues changed and, with the addition of hydrogen bonds with Tyr77, the contact residues remained unchanged. Vorinostat lost hydrogen bonds with Gln114, Ile111, and Ser78, and hydrophobic contact substitution with hydrogen bonds with Tyr39 occurred. On the other hand, it generated additional hydrophobic interactions with residues in the E6AP/IRF3 binding pocket, such as Arg136, Val60, Leu74, and Val38. These findings corroborated the simulation results, which indicated that, despite the deviation, the three compounds remained inside the E6AP/IRF3 binding region [3636 Shah M, Wadood A, Rahman Z, Husnain T. Interaction and inhibition of Dengue envelope glycoprotein with mammalian receptor DC-Sign, An in-silico approach. PLoS ONE. 2013; 8: 1-10.

37 Vande Pol SB, Klingelhutz AJ. Papillomavirus E6 oncoproteins. Virology. 2013; 445: 115-37.

38 Rietz A, Petrov D, Bartolowits M, DeSmet M, Davisson V, Androphy E. Molecular probing of the HPV-16 E6 protein alpha helix binding groove with small molecule inhibitors. PloS ONE. 2016; 11: 1-20.
-3939 Ricci-López J, Vidal-Limon A, Zunñiga M, Jimènez VA, Alderete JB, Brizuela CA, et al. Molecular modeling simulation studies reveal new potential inhibitors against HPV E6 protein. PLoS ONE. 2019; 14: 1-22.].

The SASA analysis revealed that the values fluctuated during the simulation procedure, indicating the three protein complexes were expanding [4040 Schmidt C, Macpherson JA, Lau AM, Tan KW, Fraternali F, Politis A. Surface accessibility and dynamics of macromolecular assemblies probed by covalent labeling mass spectrometry and integrative modeling. Anal Chem. 2017; 89: 1459-68.

41 Marsh JA, Teichmann SA. Relative solvent accessible surface area predicts protein conformational changes upon binding. Structure. 2011; 19: 859-867.
-4242 Candotti M, Pérez A, Ferrer-Costa C, Rueda M, Meyer T, Gelpí JL, Orozco M. Exploring early stages of the chemical unfolding of proteins at the proteome scale. PLoS Comput. Biol. 2013; 9: 1-11.]. However, the SASA value, typically constant after 50 ps, suggests that all proteins' structures stay folded and do not collapse due to damage to their hydrophobic cores [4343 Tarek M, Tobias DJ. The role of protein-solvent hydrogen bond dynamics in the structural relaxation of a protein in glycerol versus water. Eur. Biophys. J. 2008; 37: 701-709., 4444 Zhang D, Lazim R. Application of conventional molecular dynamics simulation in evaluating the stability of apomyoglobin in urea solution. Sci. Rep. 2017; 7: 1-12.]. Hydrogen bonds are required for the formation and maintenance of protein structures. In general, hydrogen bonds are classified into two types: intramolecular and protein-solvent. The intramolecular hydrogen bond's primary function is to preserve the native protein's stability and structure under folded circumstances, while the protein-solvent hydrogen bond's primary function is to promote the protein's dynamics [4343 Tarek M, Tobias DJ. The role of protein-solvent hydrogen bond dynamics in the structural relaxation of a protein in glycerol versus water. Eur. Biophys. J. 2008; 37: 701-709., 4545 Petukhov M, Rychkov G, Firsov L, Serrano L. H-bonding in protein hydration revisited. Protein Sci. 2004; 13: 2120-2129., 4646 Kuhn B, Mohr P, Stahl M. Intramolecular hydrogen bonding in medicinal chemistry. J Med Chem. 2010; 53: 2601-2611.]. The intramolecular hydrogen bond and protein-solvent parameters fluctuated for each protein complex due to the continuous movement of molecular atoms throughout the simulation time [4747 Ishak SNH, Aris SNAM, Halim KBA, Ali MSM, Leow TC, Kamarudin NHA, Masomian M, Rahman RNZRA. Molecular dynamic simulation of space and earth-grown crystal structures of thermostable T1 lipase Geobacillus zalihae revealed a better structure. Molecules. 2017; 22: 1-13.

48 Chikalov I, Yao P, Moshkov M, Latombe J-C. Learning probabilistic models of hydrogen bond stability from molecular dynamics simulation trajectories. BMC Bioinformatics. 2011; 12: 1-6.
-4949 England JL, Haran G. Role of solvation effects in protein denaturation: From thermodynamics to single molecules and back. Annu. Rev. Phys. Chem. 2011; 62: 257-77.]. The graph of the two kinds of reasonably stable hydrogen bonds indicates that the structure of such complex proteins is relatively stable and does not denature spontaneously throughout the simulation [4040 Schmidt C, Macpherson JA, Lau AM, Tan KW, Fraternali F, Politis A. Surface accessibility and dynamics of macromolecular assemblies probed by covalent labeling mass spectrometry and integrative modeling. Anal Chem. 2017; 89: 1459-68.]. Also, although the three complexes showed nearly identical movement matrix patterns, asarinin, and thiazolo complexes were significantly more rigid and stable than the control [3333 Valente RPDP, Souza RC, de Medeiros Muniz G, Ferreira JEV, de Miranda RM, E Lima AHL, et al. Using accelerated molecular dynamics simulation to elucidate the effects of the T198F mutation on the molecular flexibility of the West Nile virus envelope protein. Sci. Rep. 2020; 10: 1-6., 5050 Dash R, Ali MC, Dash N, Azad MAK, Hosen SMZ, Hannan MA, Moon IS. Structural and dynamic characterizations highlight the deleterious role of SULT1A1 R213H polymorphism in substrate binding. Int. J. Mol. Sci. 2019; 20: 1-22.]. They create stable protein-ligand complexes with the HPV-16 E6 protein, according to the simulation. Since this interaction occurs precisely in the E6AP/IRF3 binding sites, the two compounds are believed to function as concurrent competitive inhibitors of E6AP and IRF3. Disrupting the formation of the E6-E6AP complex is expected to reduce the degree of p53 degradation by inhibiting E6 protein conformational changes required for binding to p53 [1717 Sailer C, Offensperger F, Julier A, Kammer K-M, Walker-Gray R, Gold MG, et al. Structural dynamics of the E6AP/UBE3A-E6-p53 enzyme-substrate complex. Nat. Commun. 2018; 9: 1-12., 2828 Heikal MF, Putra WE, Sustiprijatno, Rifa'i M, Hidayatullah A, Ningsih FN, et al. In silico screening and molecular dynamics simulation of potential anti-malarial agents from Zingiberaceae as potential Plasmodium falciparum lactate dehydrogenase (PfLDH) enzyme inhibitors. Trop. Life Sci. Res. 2023. 34(2): 1-20.]. By inhibiting the development of the E6-IRF3 complex, those potential compounds are believed to restore signaling pathways and activate the IRF3 pathway, thus improving the overall innate immune response's performance against HPV-16 infection [1212 Tummers B, van der Burg S. High-risk human papillomavirus targets crossroads in immune signaling. Viruses. 2015; 7: 2485-506., 3636 Shah M, Wadood A, Rahman Z, Husnain T. Interaction and inhibition of Dengue envelope glycoprotein with mammalian receptor DC-Sign, An in-silico approach. PLoS ONE. 2013; 8: 1-10., 5151 Howie HL, Katzenellenbogen RA, Galloway DA. Papillomavirus E6 proteins. Virology. 2009; 384: 324-34.].

CONCLUSION

The E6 protein is a major oncoprotein linked to the development of cancer caused by HPV-16 infections. According to different molecular dynamic characteristics, asarinin and thiazolo are expected to form stable protein-ligand complexes with the HPV-16 E6 protein. While asarinin is more stable than the other compounds in terms of contact residues, the ligand conformation and movement graph indicate that thiazolo is more stable in simulation. Additionally, the SASA, hydrogen bond, and DCCM graphs support both compounds being equally effective against HPV-16 E6 protein. This research may pave the path for the future use of this compound as an E6 protein inhibitor.

Acknowledgments

We thank the Computational Laboratory, Department of Biology, Brawijaya University for providing the research instruments.

REFERENCES

  • 1
    Nurcahyanti ADR. Cervical cancer: The case in Indonesia and natural product based therapy. J. Sci. Med. 2016; 4: 1-7.
  • 2
    Domingo E, Noviani R, Noor M, Ngelangel C, Limpaphayom K, Thuan T, et al. Epidemiology and prevention of cervical cancer in Indonesia, Malaysia, the Philippines, Thailand and Vietnam. Vaccine. 2008; 26(12): 71-9.
  • 3
    Graham SV. The human papillomavirus replication cycle, and its links to cancer progression: A comprehensive review. Clin. Sci. 2017; 131: 2201-21.
  • 4
    Tan S, de Vries EG, van der Zee AG, de Jong S. Anticancer drugs aimed at E6 and E7 activity in HPV-positive cervical cancer. Curr Cancer Drug Targets. 2012; 12: 170-84.
  • 5
    Wang X, Huang X, Zhang Y. Involvement of human papillomaviruses in cervical cancer. Front Microbiol. 2018; 9: 1-14.
  • 6
    Huh WK, Joura EA, Giuliano AR, Iversen OE, de Andrade RP, Ault KA, et al. Final efficacy, immunogenicity, and safety analyses of a nine-valent human papillomavirus vaccine in women aged 16-26 years: A randomised, double-blind trial. Lancet. 2017; 390: 2143-59.
  • 7
    El-Zein M, Richardson L, Franco EL. Cervical cancer screening of HPV vaccinated populations: Cytology, molecular testing, both or none. J. Clin. Virol. 2016; 76: 62-8.
  • 8
    Pal A, Kundu R. Human papillomavirus E6 and E7: The cervical cancer hallmarks and targets for therapy. Front Microbiol. 2020; 10: 1-15.
  • 9
    Yeo-Teh NSL, Ito Y, Jha S. High-risk human papillomaviral oncogenes E6 and E7 target key cellular pathways to achieve oncogenesis. Int. J. Mol. Sci. 2018; 19: 1-27.
  • 10
    Doorbar J, Quint W, Banks L, Bravo IG, Stoler M, Broker TR, et al. The biology and life-cycle of human papillomaviruses. Vaccine. 2012; 30: 55-70.
  • 11
    Bernard X, Robinson P, Nominé Y, Masson M, Charbonnier S, Ramirez-Ramos JR, Deryckere F, Travé G, Orfanoudakis G. Proteasomal degradation of p53 by human papillomavirus E6 oncoprotein relies on the structural integrity of p53 core domain. PLoS ONE. 2011; 6: 1-10.
  • 12
    Tummers B, van der Burg S. High-risk human papillomavirus targets crossroads in immune signaling. Viruses. 2015; 7: 2485-506.
  • 13
    Reiser J, Hurst J, Voges M, Krauss P, Münch P, Iftner T, et al. High-risk human papillomaviruses repress constitutive kappa interferon transcription via E6 to prevent pathogen recognition receptor and antiviral-gene expression. J. Virol. 2011; 85: 11372-80.
  • 14
    Shah M, Anwar MA, Park S, Jafri SS, Choi S. In silico mechanistic analysis of IRF3 inactivation and high-risk HPV E6 species-dependent drug response. Sci. Rep. 2015; 5: 1-14.
  • 15
    Zanier K, Stutz C, Kintscher S, Reinz E, Sehr P, Bulkescher J, Hoppe-Seyler K, Travé G, Hoppe-Seyler F. The E6AP binding pocket of the HPV16 E6 oncoprotein provides a docking site for a small inhibitory peptide unrelated to E6AP, indicating druggability of E6. PLoS ONE. 2014; 9: 1-12.
  • 16
    Martinez-Zapien D, Ruiz FX, Poirson J, Mitschler A, Ramirez J, Forster A, et al. Structure of the E6/E6AP/p53 complex required for HPV-mediated degradation of p53. Nature. 2016; 529: 541-5.
  • 17
    Sailer C, Offensperger F, Julier A, Kammer K-M, Walker-Gray R, Gold MG, et al. Structural dynamics of the E6AP/UBE3A-E6-p53 enzyme-substrate complex. Nat. Commun. 2018; 9: 1-12.
  • 18
    Mortensen F, Schneider D, Barbic T, Sladewska-Marquardt A, Kühnle S, Marx A, Scheffner M. Role of ubiquitin and the HPV E6 oncoprotein in E6AP-mediated ubiquitination. PNAS. 2015; 112: 9872-7.
  • 19
    Messa L, Celegato M, Bertagnin C, Mercorelli B, Nannetti G, Palù G, Loregian A. A quantitative LumiFluo assay to test inhibitory compounds blocking p53 degradation induced by human papillomavirus oncoprotein E6 in living cells. Sci. Rep. 2018; 8: 1-11.
  • 20
    Banerjee NS, Moore DW, Broker TR, Chow LT. Vorinostat, a pan-HDAC inhibitor, abrogates productive HPV-18 DNA amplification. PNAS. 2018; 115: 11138-47.
  • 21
    Jayaram B, Singh T, Mukherjee G, Mathur A, Shekhar S, Shekhar V. Sanjeevini: A freely accessible web-server for target directed lead molecule discovery. BMC Bioinformatics. 2012; 13: 1-13.
  • 22
    Mirzaei H, Zarbafian S, Villar E, Mottarella S, Beglov D, Vajda S, et al. Energy minimization on manifolds for docking flexible molecules. J. Chem. Theory Comput. 2015; 11: 1063-76.
  • 23
    Trott O, Olson AJ. AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J. Comput. Chem. 2009; 31: 455-61.
  • 24
    Hidayatullah A, Putra WE, Sustiprijatno, Permatasari GW, Salma WO, Widiastuti D, et al. In silico targeting DENV2's prefusion envelope protein by several natural products' bioactive compounds. CMUJ. Nat. Sci. 2021; 20(4): 1-20.
  • 25
    Putra WE, Salma WO, Rifa’i M. Anti-inflammatory activity of Sambucus plant bioactive compounds against TNF-α and TRAIL as solution to overcome inflammation associated diseases: The insight from bioinformatics study. Nat. Prod. Sci. 2019; 25(3): 215-21.
  • 26
    Putra WE, Waffareta E, Ardiana O, Januarisasi ID, Soewondo A, Rifa'i M. Dexamethasone-administrated BALB/c mouse promotes proinflammatory cytokine expression and reduces CD4+CD25+ regulatory T cells population. Biosci. Research. 2017; 14(2): 201-13.
  • 27
    Putra WE. In silico study demonstrates multiple bioactive compounds of Sambucus plant promote death cell signaling pathway via Fas receptor. FTST. 2018; 3(2): 682-5.
  • 28
    Heikal MF, Putra WE, Sustiprijatno, Rifa'i M, Hidayatullah A, Ningsih FN, et al. In silico screening and molecular dynamics simulation of potential anti-malarial agents from Zingiberaceae as potential Plasmodium falciparum lactate dehydrogenase (PfLDH) enzyme inhibitors. Trop. Life Sci. Res. 2023. 34(2): 1-20.
  • 29
    Hidayatullah A, Putra WE, Rifa’i M, Sustiprijatno, Widiastuti D, Heikal MF, et al. Molecular docking and dynamics simulation studies to predict multiple medicinal plants’ bioactive compounds interaction and its behavior on the surface of DENV-2 E protein. KIJOMS. 2022;8(3):531-42.
  • 30
    Hidayatullah A, Putra WE, Rifa'i M, Sustiprijatno, Heikal MF, Widiastuti D, et al. Molecular docking and molecular dynamic study of multiple medicinal plants’ bioactive compounds as human papillomavirus type 16 E5 protein inhibitor. Kuwait J. Sci. 2023. 50(1A): 1-14.
  • 31
    Hidayatullah A, Putra WE, Sustiprijatno, Rifa’i M, Widiastuti D, Heikal MF, et al. Concatenation of molecular docking and dynamics simulation of human papillomavirus type 16 E7 oncoprotein targeted ligands: In quest of cervical cancer's treatment. An. Acad. Bras. Cienc. 2023; 95(Suppl. 1): 1-16.
  • 32
    Hidayatullah A, Putra WE, Sustiprijatno, Widiastuti D, Salma WO, Heikal MF. Molecular docking and molecular dynamics simulation-based identification of natural inhibitors against druggable human papilloma virus type 16 target. Trends Sci. 2023; 20(4): 1-12.
  • 33
    Valente RPDP, Souza RC, de Medeiros Muniz G, Ferreira JEV, de Miranda RM, E Lima AHL, et al. Using accelerated molecular dynamics simulation to elucidate the effects of the T198F mutation on the molecular flexibility of the West Nile virus envelope protein. Sci. Rep. 2020; 10: 1-6.
  • 34
    Aier I, Varadwaj PK, Raj U. Structural insights into conformational stability of both wild-type and mutant EZH2 receptor. Sci. Rep. 2016; 6: 1-10.
  • 35
    Chopdar KS, Dash GC, Mohapatra PK, Nayak B, Raval MK. Monte-Carlo method-based QSAR model to discover phytochemical urease inhibitors using SMILES and GRAPH descriptors. J. Biomol. Struct. Dyn. 2021; 0: 1-10.
  • 36
    Shah M, Wadood A, Rahman Z, Husnain T. Interaction and inhibition of Dengue envelope glycoprotein with mammalian receptor DC-Sign, An in-silico approach. PLoS ONE. 2013; 8: 1-10.
  • 37
    Vande Pol SB, Klingelhutz AJ. Papillomavirus E6 oncoproteins. Virology. 2013; 445: 115-37.
  • 38
    Rietz A, Petrov D, Bartolowits M, DeSmet M, Davisson V, Androphy E. Molecular probing of the HPV-16 E6 protein alpha helix binding groove with small molecule inhibitors. PloS ONE. 2016; 11: 1-20.
  • 39
    Ricci-López J, Vidal-Limon A, Zunñiga M, Jimènez VA, Alderete JB, Brizuela CA, et al. Molecular modeling simulation studies reveal new potential inhibitors against HPV E6 protein. PLoS ONE. 2019; 14: 1-22.
  • 40
    Schmidt C, Macpherson JA, Lau AM, Tan KW, Fraternali F, Politis A. Surface accessibility and dynamics of macromolecular assemblies probed by covalent labeling mass spectrometry and integrative modeling. Anal Chem. 2017; 89: 1459-68.
  • 41
    Marsh JA, Teichmann SA. Relative solvent accessible surface area predicts protein conformational changes upon binding. Structure. 2011; 19: 859-867.
  • 42
    Candotti M, Pérez A, Ferrer-Costa C, Rueda M, Meyer T, Gelpí JL, Orozco M. Exploring early stages of the chemical unfolding of proteins at the proteome scale. PLoS Comput. Biol. 2013; 9: 1-11.
  • 43
    Tarek M, Tobias DJ. The role of protein-solvent hydrogen bond dynamics in the structural relaxation of a protein in glycerol versus water. Eur. Biophys. J. 2008; 37: 701-709.
  • 44
    Zhang D, Lazim R. Application of conventional molecular dynamics simulation in evaluating the stability of apomyoglobin in urea solution. Sci. Rep. 2017; 7: 1-12.
  • 45
    Petukhov M, Rychkov G, Firsov L, Serrano L. H-bonding in protein hydration revisited. Protein Sci. 2004; 13: 2120-2129.
  • 46
    Kuhn B, Mohr P, Stahl M. Intramolecular hydrogen bonding in medicinal chemistry. J Med Chem. 2010; 53: 2601-2611.
  • 47
    Ishak SNH, Aris SNAM, Halim KBA, Ali MSM, Leow TC, Kamarudin NHA, Masomian M, Rahman RNZRA. Molecular dynamic simulation of space and earth-grown crystal structures of thermostable T1 lipase Geobacillus zalihae revealed a better structure. Molecules. 2017; 22: 1-13.
  • 48
    Chikalov I, Yao P, Moshkov M, Latombe J-C. Learning probabilistic models of hydrogen bond stability from molecular dynamics simulation trajectories. BMC Bioinformatics. 2011; 12: 1-6.
  • 49
    England JL, Haran G. Role of solvation effects in protein denaturation: From thermodynamics to single molecules and back. Annu. Rev. Phys. Chem. 2011; 62: 257-77.
  • 50
    Dash R, Ali MC, Dash N, Azad MAK, Hosen SMZ, Hannan MA, Moon IS. Structural and dynamic characterizations highlight the deleterious role of SULT1A1 R213H polymorphism in substrate binding. Int. J. Mol. Sci. 2019; 20: 1-22.
  • 51
    Howie HL, Katzenellenbogen RA, Galloway DA. Papillomavirus E6 proteins. Virology. 2009; 384: 324-34.
  • Funding:

    This research received no external funding

Edited by

Editor-in-Chief:

Paulo Vitor Farago

Associate Editor:

Andressa Novatski

Publication Dates

  • Publication in this collection
    31 May 2024
  • Date of issue
    2024

History

  • Received
    15 Nov 2022
  • Accepted
    27 Feb 2024
Instituto de Tecnologia do Paraná - Tecpar Rua Prof. Algacyr Munhoz Mader, 3775 - CIC, 81350-010 Curitiba PR Brazil, Tel.: +55 41 3316-3052/3054, Fax: +55 41 3346-2872 - Curitiba - PR - Brazil
E-mail: babt@tecpar.br