Concatenation of molecular docking and dynamics simulation of human papillomavirus type 16 E7 oncoprotein targeted ligands: In quest of cervical cancer’s treatment

Abstract

The Human papillomaviruses type 16 E7 oncoprotein is a 98-amino-acid, 11-kilodalton acidic oncoprotein with three conserved portions. Due to its interaction with the pRb-E2F complex, CKII, CKI (mostly p21), and even HDAC1, it possesses strong transformative and carcinogenic qualities that inhibit normal differentiation and cell cycle regulation. Here, we target the E7 oncoprotein using two prior research active compounds: asarinin and thiazolo[3,2-a]benzimidazole-3(2H)-one,2-(2-fluorobenzylideno)-7,8-dimethyl (thiazolo), and valproic acid as a control. We are performing molecular docking followed by molecular dynamic analysis. By acting as competitive inhibitors in the binding site, it was hypothesized that both drugs would inhibit E7-mediated pRb degradation and E7-mediated p21 degradation, resulting in decreased cell cycle progression, immortalization, and proliferation. In addition, we expect that the direct inhibitory action of valproic acid in E7 will target the CKII-mediated phosphorylation pathway necessary for destabilizing p130 and pRb. According to the results of the dynamic simulation, stable interactions exist between every compound. Despite the instability of E7 protein, stability results indicate that both natural chemicals are preferable, with thiazolo outperforming valproic acid.

Key words
anti-cancer; cervical cancer; HPV-16 E7 protein; molecular docking; molecular dynamics simulation

INTRODUCTION

Human papillomaviruses (HPVs) are a DNA virus family that includes more than 200 different types with a strong preference for squamous epithelium (Tomaić 2016, Bzhalava et al. 2013, Graham 2017). However, they have a similar genetic make-up, various HPVs attack epithelia at different anatomic sites. The anogenital and oral mucosa infections are caused by around 30 different HPV types, which can be further classed as low-risk or high-risk based on the clinical prognosis of the lesions they cause (Münger et al. 2001, Ohlenschläger et al. 2006). While HPVs that are low-risk generate benign epithelial hyperplasia, HPVs that are high-risk create lesions with a high proclivity to proceed to malignancy. Cervical cancers are commonly linked with high-risk HPV infection in nearly all instances, and two viral proteins, E6 and E7, which are persistently expressed in the transformed phenotype, are essential for developing and maintaining the transformed phenotype, respectively (Shin et al. 2009, Münger et al. 2001). Most research on HPV inhibitors has focused on those oncoproteins due to their well-characterized activities pathways and highly differentiate high-risk from low-risk variants (Aarthy & Singh 2021).

The HPV-16 E7 oncoprotein is a 98-amino-acid-long acidic oncoprotein with a calculated molecular weight of around 11 kDa. It possesses strong transforming and carcinogenic capabilities that interfere with normal differentiation and cell cycle control (Valdovinos-Torres et al. 2008, Shin et al. 2009). The E7 protein’s three-dimensional structures have been identified, revealing an inherently disordered amino terminus and a novel metal-binding fold in their carboxyl terminus. Conserved regions 1, 2, and 3 (CR1, CR2, and CR3) are composed of two different functional domains at the N-terminus and a single functional domain at the C-terminus (Pang & Thierry 2013, Mok et al. 2018). The CR1 (amino acids 1–18) and CR2 (amino acids 19–42) segments are naturally unfolded with a high degree of flexibility and share structural and functional similarities with the adenovirus E1A protein and the SV40 large T-antigen (Tomaić 2016, Chellappan et al. 1992, García-Alai et al. 2007, Nogueira et al. 2017). The CD1 domain is required for E7 to cause S-phase progression and cellular transformation, whereas the CR2 homology domain comprises the retinoblastoma (pRB) tumor suppressor core-binding site, the LxCxE domain, a casein kinase II phosphorylation site, and its recognition site (Tomaić 2016, Bello-Rios et al. 2021, Knapp et al. 2009, Jones et al. 1997). The CR1 and CR3 have previously been implicated in destabilizing pRB-related proteins p107 and p130 by E7. The C-terminal domain contains the CR3 region, which is structurally similar to the C-terminal domain of E6 oncoprotein and contains two CxxC zinc-binding motifs separated by 29 residues required for zinc-dependent dimerization and protein stabilization (Ohlenschläger et al. 2006, Mok et al. 2018). Additionally, this domain serves as a platform for direct interaction with various tumor suppressor proteins, including p21, p27, and HDAC1. Also, the CR3 region interacts with the C-terminal region of pRb, supporting its degradation (Tomaić 2016, Todorovic et al. 2011). The full-length E7 protein has a 100-fold greater affinity for pRb than the E7 protein lacking CR1/CR2, indicating a critical function for CR3 in pRb interactions and degradations (Syafitri et al. 2021).

Numerous treatment options targeting the HPV’s main oncoproteins have been developed, including therapeutic vaccines targeting E7 protein, genome editing with antisense oligonucleotides, ribozymes, siRNA, immunotherapy with synthetic E7 DNA sequences, a proteasome inhibitor, and tumor-infiltrating T cells (Syafitri et al. 2021, Tan et al. 2012, Pal & Kundu 2020). However, they are too expensive for the world’s poorest and most vulnerable population, mainly in the middle- and low-income countries (Huh et al. 2017). Despite decades of study, no effective and conclusive therapy for chronic HPV infection has been established yet (Toots et al. 2017). As a mega biodiversity country, Indonesia offers abundant resources that may be developed medically and cost-effectively as an alternative to currently available pharmaceuticals or therapies (Syafitri et al. 2021, Salim & Munadi 2017). We examined asarinin and thiazolo from Zanthoxylum spp. and Myristica fragrans, respectively, as possible HPV-16 E7 protein inhibitors based on our prior study on HPV-16 E5 and E6 protein inhibitors (data not shown) (Asika et al. 2016, Zhang et al. 2002). There is currently no medicine that precisely targets the E7 protein, thus we utilize a histone deacetylase inhibitor (HDACi), valproic acid, as a control in this in silico study. Valproic acid is generally used to treat epilepsy, bipolar disorder, and to avoid migraines, as well as menstrual abnormalities, polycystic ovaries, and hyperandrogenism (Aronson 2015, Andoh et al. 2020, Feng et al. 2016). Some studies also suggest that it also has anticancer action, notably as an antineoplastic drug identified in 1997 by promoting apoptosis and inhibiting angiogenesis in tumor cells (Cinatl et al. 1997). It suppresses the development of HPV-positive SiHa, CasKi, and HeLa cell lines in vitro, as well as the expression of E7 mRNA. However, its structural effects on the HPV-16 E7 protein remain unknown (Tan et al. 2012, de la Cruz-Hernández et al. 2007, Dai et al. 2017, Feng et al. 2016). This study investigates the potential of asarinin and thiazolo as HPV-16 E7 protein inhibitor candidates using a few molecular dynamic parameters to determine their further potential.

MATERIALS AND METHODS

Data retrieval and pre-docking screening

The HPV-16 E7 amino acid sequence was obtained from UniProt under the accession number P03129, and the three-dimensional structure of the protein was simulated using the I-TASSER webserver (https://zhanglab.dcmb.med.umich.edu/I-TASSER/). The 3D model was selected based on its position in C- and TM-score values rankings, among other factors. Using the Ramachandran Plot Server (https://zlab.umassmed.edu/bu/rama/), we evaluated the model’s structure, with the favored percentage coming out to 69.767%. In this study, asarinin (ID 11869417) and thiazolo (ID 1823738) were examined, as well as valproic acid (CID: 3121) as a positive control.

Molecular docking process

We used AutoDock Vina integrated within PyRx (https://pyrx.sourceforge.io/) to dock our modeled HPV-16 E7 protein (Hidayatullah et al. 2023). In this study, we look at the whole structure of the target protein, with a molecular coverage area of 36.8165×29.3955×42.2397 Å and center coordinates 51.3889×51.4450×51.5934 Å. The major docking parameters are the molecule’s affinity (measured in kcal/mol), the binding site pocket, and the interaction between the protein and the ligand as our previous study (Widiastuti et al. 2023, Putra et al. 2023).

Molecular dynamics simulation

We evaluated all of the protein-compound complexes at 37 °C, 1 atm, pH 7.4, and 0.9 % salinity, which are typical physiological conditions. For 1000 picosecond simulations, the protein and ligand complexes’ structures were generated. A macro program called Md run is used to run a molecular dynamics simulation. The simulation results were examined using the md analyze and md analyzers macro programs through yasara (Hidayatullah et al. 2022).

Target protein prediction

Three ligands chosen needed to be validated its protein target, so that the interaction between ligands and protein could be confirmed. To predict the protein target, we utilized SwissTarget webserver (http://www.swisstargetprediction.ch/). The isomeric SMILES of asarinin, thiazolo, and canonical SMILES of valproic acid were used as input (Table I).

Network analysis and GO prediction

To predict the interaction of three candidate ligands towards protein target in HPV cases, we performed networking analysis using STITCH database webserver (http://stitch.embl.de/). The webserver automatically generating the Gene Ontology including biological process, molecular function, and KEGG pathway (Heikal et al. 2021).

RESULTS

Before 3D analysis using PyMOL, we align the structure of the E7-ligand complexes before and after the MD simulation with a cycle value parameter of 5.0 and a cut-off value of 2.0 (Fig 1). The alignment results showed that the three complexes deviated around 3-4 Å. As expected, the most significant deviation was experienced by the CR1 and CR2 regions in each complex. Alignment result suggests the E7-thiazolo complex experienced the slightest deviation during the MD simulation process (RMSD: 3.910 Å), while the E7-valproic acid complex experienced the most significant structural deviation during the simulation process (RMSD: 4.942 Å). Deviation in the CR2 region is thought to play a role in causing changes in amino acid residues that interact with the three compounds tested, predominantly shifting more towards the CR3 region. These results are also reflected in the 2D visualization results, which show that less than 55% of pre-MD interactions remain intact after the simulation. Valproic acid only retains its 10% interaction before simulation, whereas asarinin and thiazolo retains some substantial amount of interaction after the simulation, around 40% for asarinin, and 54.5% for thiazolo.

Figure 1
The 3D and 2D structure visualization of HPV-16 E7 and (a) asarinin, (b) thiazolo, (c) valproic acid complex. The cyan structures in 3D visualization indicate the protein-ligand complex before MD simulation, while the green structures are after MD simulation. The colored diagram in 2D indicates protein-ligand complex after MD, while the greyscale diagram indicates protein-ligand complex before MD.

Although the three protein-ligand complexes tested contained less than 55% of the conserved residues, post-MD simulation results revealed that all tested compounds interacted with the same domain as before, with the exception of thiazolo, which interacts with additional domains and critical residues after simulation. While asarinin retains contacts with pRb binding domains (Tyr52, Ile76, Arg77), CKII phosphorylation recognition site (Glu35), and p21 binding residues (Ile54), a change in the interacting residues shows when compared to interactions in these three domains prior to the simulation. Valproic acid retains its contact with the pRb binding domain (Tyr52, Arg77) but switches its contact with the CKII phosphorylation recognition sites (Glu35, Asp36, Glu37) to direct interaction with the CKII phosphorylation site (Ser32). On the other hand, thiazolo interacts not only with the p21 binding domain (Ile54, Val55, Leu82), the pRB binding domain (Asn53, Val55), and the CKII phosphorylation recognition site (Glu35, Asp36) but also with an additional critical residue, the HDAC1 binding residue (Leu82).

We also noted that asarinin and thiazolo also remain in contact with the protein’s hydrophobic core. The potential shifting of asarinin causes the compound to interact with more hydrophobic core residues (Ile76, Val69, Ile54). After simulation, the thiazolo interacts with two hydrophobic core residues, namely Ile54 and Leu82. It retains the same amount of contact with hydrophobic core residues as before simulation, but interaction with Ile76 is now with Leu82.

Molecular dynamics simulation for 1000 ps was carried out to predict the stability, interaction pattern, and behavior of each protein-ligand complex with parameters of the physiological condition of the infected person. The energy potential graph shows that the E7-valproic acid complex has the lowest trend value, but it is insignificant compared to the other two complexes (Fig. 2). In addition, the three complexes also show a horizontal trend of graph values throughout the simulation period, around -3.13e⁵ kJ/mol, proving that the modeled complexes were energetically stable during the MD run.

Figure 2
Potential energy for HPV-16 E7 protein complex with asarinin (red line), thiazolo (green line), and valproic acid (black line) over a 1000 picosecond of simulation.

The RMSD profile (Fig. 3A) of the three control complexes shows a gradual increase in the RMSD value in the first 300 picoseconds and reaches equilibrium. The mean RMSD values of the three complexes were in the range of 3.940-4.66 Å. It is clear that the RMSD graph of the complex E7- thiazolo has a lower value than the other two complexes, as well as fluctuations the minimum during the simulation period. However, all complexes showed the most prolonged equilibration period, confirming the protein’s prolonged stability. The root means square fluctuation (RMSF) was then calculated to understand how the mutations affected protein backbone motion.

Figure 3
(a) RMSD plot, (b) RMSF plot, (c) ligand movement plot, and (d) ligand conformation for HPV-16 E7 protein complex with asarinin (red line), thiazolo (green line), and valproic acid (black line) over a 1000 picosecond of simulation.

The RMSF plot (Fig. 3b) shows that the most significant fluctuations occur in the CR1 and CR2 regions of the three protein-ligand complexes. The CR3 region was significantly more stable than CR1 and CR2 over the simulation, universally. Apart from those results, the graph of thiazolo (1.634±0.843 Å) appears to be more stable than asarinin (2.022±1.041 Å) and valproic acid (1.875±0.882 Å). The most significant fluctuations occurred in the region around the LxCxE domain, the CKII phosphorylation site, and its recognition domain, also with the N-terminal of the CR3 region. Those regions are the dominant interaction areas of the three tested compounds before the simulation, predominantly around residues number 30-38 and 52-55 (Table II). However, after simulation, all tested compounds mainly shifted to interact with the CR3 region instead, with better structural stability. These results were also confirmed by the results of the ligand movement simulation (Fig. 3c), which showed extensive ligand movement beyond 5 Å each. Thiazolo offered a stability pattern similar to the previous simulation results and has the most stable ligand movement value among the three compounds tested (5,635±1,049 Å). The most significant ligand movement was experienced by asarinin (7,523±1,339 Å), which experienced significant fluctuations in the first 100 ps of the simulation period and on the timeframes of 750 ps and 950 ps. The three tested ligands appeared to have settled at the most stable binding site after 100 ps simulation with minor fluctuations during the simulation period. The ligand conformation plot (Fig. 3d) shows asarinin has the lowest value of conformational change (1.661±0.312 Å). In contrast to the results of the two previous simulations, which were highly favored thiazolo, this compound had the highest conformational change value (1,859±0.254 Å). The three compounds tested predominantly had conformational fluctuations around 1.5 to 2 Å.

The 1000 ps SASA simulation revealed that none of the three protein-ligand complexes expanded or shrunk rapidly over the simulation period (Fig. 4a). The E7- thiazolo complex had the flattest and stable trend of SASA values with fluctuations just below 7000 2 throughout the simulation period lasts. The E7-asarinin and E7-valproic acid complexes experienced two significant changes in two different timeframes. The first fluctuation occurs at 0-300 ps, while the second fluctuation occurs after 800 ps. The E7-asarinin complex has the most drastic fluctuation under 300 ps among all of them. After 800 ps, E7-asarinin and E7-valproic acid start to divert into the different trajectory, relative to flat E7- thiazolo complex. The E7-valproic acid complex tends to shrink after 800 ps. It is reflected with its SASA value, which tends to decrease to 6500 Šat the end of the simulation. The E7-asarinin complex tends to expand, imaged with the SASA value trend to increase to around 7500 Ų at the end of the simulation period.

Figure 4
(a) SASA plot, (b) Rg plot, (c) intramolecular hydrogen bond plot, and (d) protein-solvent hydrogen bond plot for HPV-16 E7 protein complex with asarinin (red line), thiazolo (green line), and valproic acid (black line) over a 1000 picosecond of simulation.

The radius of gyration (Rg) values for the HPV-16 E7 protein’s backbone atoms were determined and plotted against simulation time (Fig. 4b). In general, the fluctuations of the three protein-ligand complexes were 13.7-14.6 Å. It is clear that E7- thiazolo complex has the lowest radius value and is stable during the simulation period (13.80±0.14 Å). Another two protein-ligand complexes fluctuate more vigorously, especially the E7-asarinin complex, which jump from 13.5 to 14.6 Å at 100-250 ps timeframe and continue to fluctuate even further to 14.76 Å at the end of the simulation period, around 10% more compared to its initial state. The E7-asarinin complex is also predicted to have an upright trajectory corresponding to its SASA plot. The E7-valproic acid, on the other hand, did not show a plateau graphic pattern showing the stability of the Rg value either and fluctuated around 13.36-14.38 Å. However, it had a lower and more stable mean value (14.02±0.23 Å) compared to the E7-asarinin complex (14.20±0.28 Å).

The intramolecular hydrogen bond plot shows fluctuating values (Fig. 4c). The highest fluctuation was observed in the E7-valproic acid complex, while the E7-asarinin complex and E7- thiazolo complex showed the highest fluctuation, but also a horizontal trajectory even though the graph is not strictly planar. Although the intramolecular hydrogen plot shows relatively significant fluctuation values, it is not apparent in the protein-solvent hydrogen plot (Fig. 4D). The three protein-ligand complexes showed a flat graphic pattern, except for the E7-valproic acid complex, which showed a planar trajectory pattern but tended to decrease even though it was less significant during the 1000 ps simulation.

The output of the DCCM plot revealed that all three complexes exhibited the same pattern of correlated movement (Fig. 5). A series of correlating motions concentrated on the CR3 region, specifically in the β₁ and β₂ structures, as well as an α-helix immediately behind them. The E7-valproic acid complex appears to have far more correlated mobility than the other two complexes. Although the movement intensities of the E7-asarinin and E7-thiazolo complexes were almost comparable, the E7-thiazolo complex exhibited a relatively lower movement intensity, indicated by lower yellow intensity on the plot.

Figure 5
Dynamic cross-correlation matrix for each protein-ligand complex. Yellow represents correlated motions, and blue represents anticorrelated motions.

In order to validate the function and cascade pathway of three candidate ligands, we performed protein target prediction (Supplementary Material - Figure S1), continued with networking analysis (Fig. S2 and S3). The protein target of asarinin shows that MCL1 (Induced myeloid leukemia cell differentiation protein Mcl-1) has the highest probability score as target (0.78). While target protein of thiazolo has similar probability score around 0.1, including ENPP1, ADORA1, ADORA2A, ADORA2B, PTPsigma, DHODH, CRHR1, ELANE, GPR39, CCNB3. Valproic target protein probability score varied between 0.07 to 0.14. Top three highest target protein score are PPARD, FFAR1, and FABP2.

Target proteins of asarinin included in various pathway, mainly in stress regulation, binding with protein and transcription factor, and playing roles in cancer pathway, renal and colorectal carcinoma, hepatitis B and HIF signaling pathway. While GO analysis of thiazolo’s target protein mainly acts in the phosphorylation activity, signaling activity, and included in the cell cycle, viral and cancer pathway (Table III and IV). The terminology did not show literal name as HIV inhibitor pathway.

The network analysis did not show any interaction between compounds and protein targets. This probably because there is no in vivo/in vitro experiment proven the protein targets of asarinin and thiazolo directly (Fig. S2 and S3).

DISCUSSION

Molecular docking and dynamics simulations indicate that asarinin, thiazolo, and valproic acid are primarily localized in the C-terminal of CR2 and CR3 regions. It is assumed that these two areas are more stable and less flexible than the CR1 and N-terminal CR2 regions (Chellappan et al. 1992, Nogueira et al. 2017, Calçada et al. 2013). The very flexible and unstructured nature of the CR1-CR2 region was further demonstrated by the comparatively high pre-and post-MD superimposed root mean square deviation alignment values (>2 Å), also further by the RMSF plot. The flexibility of the two N-terminal domains is also suggested to explain why valproic acid loses 90% of its contacts with E7 protein residues following simulation, whereas the other two compounds lose more than 50% of their initial associations with E7 protein residues. It suggests that the main reason for supporting the main important domain for E7 protein oncogenic activity is that the LxCxE motif in CR2 is available for interaction with pRb without requiring extensive structural reorganization (Ohlenschläger et al. 2006).

The three simulated compounds interacted with two key residue complexes, the pRB binding site and the CKII recognition domain, in a very similar manner. Although the major contact between E7 protein and pRB occurs in the LxCxE domain of the CR2 region, an interaction domain with pRB exists in the CR3 area as well. The interaction in the CR3 area exhibits a lesser affinity when compared to the interaction in the LxCxE domain. However, in the absence of interaction with CR3, the E7 protein’s affinity for pRB can diminish significantly. The three compounds interact primarily with the pRB binding site in the CR3 region through Tyr52, Asn53, Val55, Ile76, and Arg77. All contacts at the binding site are mediated by hydrophobic interactions, which are more robust than hydrogen bonds in the context of weak molecular force. Thus, the interaction between asarinin and thiazolo, and valproic acid with protein E7 is relatively stable and is predicted to be capable of decreasing its affinity for pRb, resulting in a slower rate of degradation of the pRb-E2F complex, which is inadequate to offset the effects of pRb-induced cell cycle arrest (Nogueira et al. 2017, Todorovic et al. 2011).

According to the simulation results, both natural compounds exhibited sustained interactions with the CKII recognition domain (aa 33–37). The domain is a five-residue expansion of the C-terminal CKII phosphorylation site. Disruption of the recognition domain disrupts the phosphorylation of pRb at Ser31 and Se32, resulting in a reduction in the affinity of E7 protein for the RbAB domain on pRb protein. This phosphorylation step is critical in the method by which the E7 protein regulates cycle cell control since phosphorylated pRb by CKII interacts with E2F transcription factors, repressing their activity and reducing DNA synthesis (Chemes et al. 2010, Basukala et al. 2019). Prior to that, the CKII required recognition of the specific location (₃₃EEEDE₃₇) in order to interact with E7 and phosphorylate pRb. Inhibition at this location impairs the protein’s phosphorylation ability (Phelps et al. 1992). However, the importance of CKII-mediated E7 phosphorylation on pRB binding and instability has been disputable (Chemes et al. 2010, Genovese et al. 2008, Firzlaff et al. 1991).

Additionally, valproic acid interacts with one of the CKII phosphorylation site’s serine residues (Ser32). The interaction, which is mediated by hydrophobic interaction and a moderate strength hydrogen bond, is thought to be relatively stable and thus would directly reduce E7’s affinity for CKII, impairing its phosphorylation capabilities even further, resulting in even less pRb being phosphorylated, impairing the cells’ immortality and replication (Jones et al. 1997, Chemes et al. 2010, Basukala et al. 2019, Genovese et al. 2008, Phillips & Vousden 1997). Additionally, interventions on this group of domains disrupt the CKII-mediated phosphorylation process necessary for destabilizing Retinoblastoma-like protein 2 (RBL-2) or p130, which is essential for promoting infected cells into the S phase via E2F-mediated transactivation (Genovese et al. 2008, Cobrinik 2005, McLaughlin-Drubin & Münger 2009).

Within the CR3 region, asarinin and thiazolo hypothesized to interfere with dimerization and stabilization mechanism of the E7 homodimer complex due to their interaction with the hydrophobic core of the CR3 region. They are critical in E7 protein dimer structure stabilization, particularly in the intermolecular contacts between the two-stranded antiparallel β-sheets on each monomer (Ohlenschläger et al. 2006, Todorovic et al. 2011, Liu et al. 2006). Meanwhile, thiazolo is also predicted to interact with a key residue of HDAC1 binding site in the region; Leu82, which overlapped with E7’s hydrophobic core residue. It is hypothesized that disrupting Leu82 and Leu83 would abolish HDAC1’s association with E7 protein and its capability to transform because altering chromatin structure occurs via deacetylation of histones, which play a critical role in HPV E7-associated transcriptional control and can extend the life of infected cells (Brehm et al. 1999, Lourenço de Freitas et al. 2021, Longworth & Laimins 2004, Helt & Galloway 2001).

Meanwhile, the p21 binding domain in CR3 overlaps with the low-affinity pRb binding region, likely since a portion of p21’s C-terminus shares a sequence with the pRb’s C-domain (Ohlenschläger et al. 2006, Liu et al. 2006). Thiazolo interacts with more residues (Ile54, Val55, Leu82) than asarinin, which may affect their inhibitory capacity, which is favored for thiazolo. Both of them are predicted to reverse p21-mediated inhibition of cyclin E–CDK2 and cyclin-A–CDK2 kinase activity and contribute to HPV E7’s ability to reestablish and sustain CDK2 activity in differentiating keratinocytes, thereby allowing for viral replication by trying to interfere with CDK2 that could significantly affect Rb-dependent E2F transcription regulation, as well as pRb phosphorylation. It is hypothesized that disrupting the E7-p21 interaction would impair the E7 protein’s capacity to stimulate DNA synthesis in differentiated cells (Shin et al. 2009, Jones et al. 1997, Parveen et al. 2016, Boshoff & Weiss 2002, Lehoux et al. 2009).

The potential energy in the three complexes was analyzed, and it was discovered that the three complexes remained stable over the simulation time. It means that the three complexes’ electrostatic and van der Waals interactions remained constant during simulation, suggesting that no anomalies occurred during the simulation process (Bavi et al. 2016, Albaugh et al. 2016). When compared to the other two complexes, the RMSD analysis suggests that the E7- thiazolo complex has the most stable backbone (Aier et al. 2016). During the simulation time, however, all three displayed a similar stabilization pattern, showing that they did not suffer considerable deformation after reaching their equilibrium state, despite their value exceeding 2 Å, owing to the instability of CR1 and CR2, as corroborated by the RMSF plot (Aier et al. 2016, Sivaramakrishnan et al. 2020, Ohlenschläger et al. 2006, Tomaić 2016). The RMSF plot reveals that all compounds’ initial binding site residues are highly flexible, contributing to extensive ligand movement and conformational changes. Each of those compounds changes conformation and location in response to the dynamics of its binding site during the simulation process (Nogueira et al. 2017, Calçada et al. 2013, Zhu et al. 2017, p.1).

Protein integrity indicators such as SASA, Rg values and hydrogen bond diagrams reveal that the integrity of all protein-ligand complexes is rather constant over the simulation period and exhibits no indications of denaturation (Zhang & Lazim 2017). These factors favor E7- thiazolo yet again. The SASA plot reveals that none of the complexes change dramatically in a short amount of time, but E7-asarinin has an upward trajectory, indicating that it is about to expand. On the other hand, E7-valproic acid displays a downward trajectory, indicating that it has reduced in size relative to its natural structure (Zhang & Lazim 2017, Kamaraj & Purohit 2013, Sneha & George Priya Doss 2016). The Rg analysis, which assesses the rigidity of the protein structure, shows that the thiazolo, and valproic acid complexes are stiffer than the asarinin complex because they have a lower Rg value than asarinin (Dash et al. 2019, Sneha & George Priya Doss 2016). Intramolecular hydrogen bond analysis revealed relatively extensive movement.

In contrast, protein-solvent hydrogen bond analysis revealed a relatively stable pattern, presumably due to critical factors such as E7’s size, CR1, and CR2 that naturally unfolded, and they moved continuously over the simulation time to rearrange themselves to achieve equilibrium (Liu et al. 2006, Ohlenschläger et al. 2006, Petukhov et al. 2004, Kuhn et al. 2010). However, they still retain their hydrophobic core, retaining the structure and staying folded (Mallamace et al. 2018, Fitzpatrick et al. 2011). Because all secondary structures are located within the CR3 region, the DCCM analysis demonstrates that all associated motions are centered inside this region. Given that both the asarinin and thiazolo complexes exhibit a similar degree of correlated movement, indicating that they are less elastic than the valproic acid complex, it is reasonable to presume that both complexes are more stable than valproic acid (Valente et al. 2020, Dash et al. 2019).

The protein targets prediction and networking analysis successfully generating some terminology related to cancer pathway and phosphorylation activity. Even though only one terminology related to viral carcinogenic found in thiazolo’s gene ontology, this could be the potential of carcinogenesis development initiated by viral infection. Each virus promotes carcinogenesis in its unique way. The majority of DNA oncogenic viruses encode oncogenes that change infected cells, with p53 and pRB being the most common targets. In addition, for certain viruses, including HBV and HPV, viral DNA integration into the human genome can play a key role in tumor growth. Because viral integration necessitates the breakage of both viral and host DNA, the rate of integration is thought to be correlated with DNA damage levels. Endogenous and exogenous causes, such as virus-induced inflammation or co-infections with other agents, environmental agents, and other factors, can cause DNA damage. Cancer usually takes years or decades to develop after an infection (Chen et al. 2014). Cervical cancer, penile cancers, and several other anogenital and head and neck cancers are linked to high-risk human papillomaviruses (HPV) 16 and 18 (some other -HPV variants are also carcinogenic) (Boshart et al. 1984, Dürst et al. 1983).

CONCLUSIONS

The HPV-16 E7 protein is one of the most well-known oncoproteins in cervical cancer cases directly involved in cell cycle progression, proliferation, and immortalization through its interaction with pRB and the pRB-related p107 and p130, CKI such as p21, and histone deacetylase. Through molecular docking and dynamic simulation, we directly target it with asarinin and thiazolo. We predict both compounds could inhibit E7-mediated pRb degradation as well as E7-mediated p21 degradation, which leads to cell cycle progression, immortalization, and proliferation. Also, we predict that the direct inhibition mechanism of valproic acid in E7 is to target the CKII-mediated phosphorylation process necessary for destabilizing p130 along with pRb. Dynamic simulation results suggest that all compounds had stable interactions. Despite the unstable nature of E7 protein, the stability results favor both natural compounds, particularly thiazolo, over valproic acid.

SUPPLEMENTARY MATERIAL

Figure S1-S3.

ACKNOWLEDGMENTS

This study was funded by PNBP Universitas Negeri Malang with contract number 5.3.619/UN32.14.1/LT/2021 and 5.3.517/UN32.14.1/LT/2021 (WEP). Authors thank Assoc. Prof. Adawiyah Suriza Shuibi, Ph.D. of Institute of Biological Sciences, Faculty of Science, Universiti Malaya, Kuala Lumpur, Malaysia for critical advices for manuscript improvement.

REFERENCES

• AARTHY M & SINGH SK. 2021. Interpretations on the interaction between protein tyrosine phosphatase and E7 oncoproteins of high and low-risk HPV: A computational perception. ACS Omega 6: 16472-16487.
• AIER I, VARADWAJ PK & RAJ U. 2016. Structural insights into conformational stability of both wild-type and mutant EZH2 receptor. Sci Rep 6: 34984.
• ALBAUGH A ET AL. 2016. Advanced potential energy surfaces for molecular simulation. J Phys Chem B 120: 9811-9832.
• ANDOH M, IKEGAYA Y & KOYAMA R. 2020. Chapter nine - Microglia in animal models of autism spectrum disorders. In: Ilieva M & Lau WK-W (Eds), Progress in Molecular Biology and Translational Science, Academic Press, p. 239-273.
• ARONSON JK. 2015. Meyler’s side effects of drugs: The international encyclopedia of adverse drug reactions and interactions, 16th ed., Amsterdam Boston Heidelberg: Elsevier.
• ASIKA AO, ADEYEMI OT, ANYASOR GN, GISARIN O & OSILESI O. 2016. GC-MS Determination of bioactive compounds and nutrient composition of Myristica fragrans Seeds. J Herbs Spices Med Plants 22: 337-347.
• BASUKALA O, MITTAL S, MASSIMI P, BESTAGNO M & BANKS L. 2019. The HPV-18 E7 CKII phospho acceptor site is required for maintaining the transformed phenotype of cervical tumour-derived cells. PLOS Pathog 15: e1007769.
• BAVI R, KUMAR R, CHOI L & LEE KW. 2016. Exploration of novel inhibitors for bruton’s tyrosine kinase by 3D QSAR modeling and molecular dynamics simulation. PLOS ONE 11: e0147190.
• BELLO-RIOS C, MONTAÑO S, GARIBAY-CERDENARES OL, ARAUJO-ARCOS LE, LEYVA-VÁZQUEZ MA & ILLADES-AGUIAR B. 2021. Modeling and molecular dynamics of the 3D structure of the HPV16 E7 protein and its variants. Int J Mol Sci 22: 1400.
• BOSHART M, GISSMANN L, IKENBERG H, KLEINHEINZ A, SCHEURLEN W & ZUR HAUSEN H. 1984. A new type of papillomavirus DNA, its presence in genital cancer biopsies and in cell lines derived from cervical cancer. EMBO J 3: 1151-1157.
• BOSHOFF C & WEISS R. 2002. Aids-related malignancies. Nat Rev Cancer 2: 373-382.
• BREHM A, NIELSEN SJ, MISKA EA, MCCANCE DJ, REID JL, BANNISTER AJ & KOUZARIDES T. 1999. The E7 oncoprotein associates with Mi2 and histone deacetylase activity to promote cell growth. EMBO J 18: 2449-2458.
• BZHALAVA D, GUAN P, FRANCESCHI S, DILLNER J & CLIFFORD G. 2013. A systematic review of the prevalence of mucosal and cutaneous human papillomavirus types. Virology 445: 224-231.
• CALÇADA EO, FELLI IC, HOŠEK T & PIERATTELLI R. 2013. The heterogeneous structural behavior of E7 from HPV16 revealed by NMR spectroscopy. ChemBioChem 14: 1876-1882.
• CHELLAPPAN S, KRAUS VB, KROGER B, MUNGER K, HOWLEY PM, PHELPS WC & NEVINS JR. 1992. Adenovirus E1A, simian virus 40 tumor antigen, and human papillomavirus E7 protein share the capacity to disrupt the interaction between transcription factor E2F and the retinoblastoma gene product. Proc Natl Acad Sci U S A 89: 4549-4553.
• CHEMES LB, SÁNCHEZ IE, SMAL C & DE PRAT-GAY G. 2010. Targeting mechanism of the retinoblastoma tumor suppressor by a prototypical viral oncoprotein: Structural modularity, intrinsic disorder and phosphorylation of human papillomavirus E7. FEBS J 277: 973-988.
• CHEN Y, WILLIAMS V, FILIPPOVA M, FILIPPOV V & DUERKSEN-HUGHES P. 2014. Viral carcinogenesis: Factors inducing DNA damage and virus integration. Cancers 6: 2155-2186.
• CINATL J, CINATL J, DRIEVER PH, KOTCHETKOV R, POUCKOVA P, KORNHUBER B & SCHWABE D. 1997. Sodium valproate inhibits in vivo growth of human neuroblastoma cells. Anticancer Drugs 8: 958-963.
• COBRINIK D. 2005. Pocket proteins and cell cycle control. Oncogene 24: 2796-2809.
• DE LA CRUZ-HERNÁNDEZ E, PÉREZ-CÁRDENAS E, CONTRERAS-PAREDES A, CANTÚ D, MOHAR A, LIZANO M & DUEÑAS-GONZÁLEZ A. 2007. The effects of DNA methylation and histone deacetylase inhibitors on human papillomavirus early gene expression in cervical cancer, an in vitro and clinical study. Virol J 4: 18.
• DAI L, CAO Y, JIANG W, ZABALETA J, LIU Z, QIAO J & QIN Z. 2017. KSHV co-infection down-regulates HPV16 E6 and E7 from cervical cancer cells. Oncotarget 8: 35792-35803.
• DASH R, ALI MC, DASH N, AZAD MAK, HOSEN SMZ, HANNAN MA & MOON IS. 2019. Structural and dynamic characterizations highlight the deleterious role of SULT1A1 R213H polymorphism in substrate binding. Int J Mol Sci 20: 6256.
• DÜRST M, GISSMANN L, IKENBERG H & ZUR HAUSEN H. 1983. A papillomavirus DNA from a cervical carcinoma and its prevalence in cancer biopsy samples from different geographic regions. Proc Natl Acad Sci U S A 80: 3812-3815.
• FENG S, YANG Y, LV J, SUN L & LIU M. 2016. Valproic acid exhibits different cell growth arrest effect in three HPV-positive/negative cervical cancer cells and possibly via inducing Notch1 cleavage and E6 downregulation. Int J Oncol 49: 422-430.
• FIRZLAFF JM, LÜSCHER B & EISENMAN RN. 1991. Negative charge at the casein kinase II phosphorylation site is important for transformation but not for Rb protein binding by the E7 protein of human papillomavirus type 16. Proc Natl Acad Sci U S A 88: 5187-5191.
• FITZPATRICK AW, KNOWLES TPJ, WAUDBY CA, VENDRUSCOLO M & DOBSON CM. 2011. Inversion of the balance between hydrophobic and hydrogen bonding interactions in protein folding and aggregation. PLOS Comput Biol 7: e1002169.
• GARCÍA-ALAI MM, ALONSO LG & DE PRAT-GAY G. 2007. The N-Terminal Module of HPV16 E7 Is an Intrinsically Disordered domain that confers conformational and recognition plasticity to the oncoprotein. Biochemistry 46: 10405-10412.
• GENOVESE NJ, BANERJEE NS, BROKER TR & CHOW LT. 2008. Casein kinase ii motif-dependent phosphorylation of human papillomavirus E7 protein promotes p130 degradation and s-phase induction in differentiated human keratinocytes. J Virol 82: 4862-4873.
• GRAHAM SV. 2017. The human papillomavirus replication cycle, and its links to cancer progression: a comprehensive review. Clin Sci 131: 2201-2221.
• HEIKAL MF ET AL. 2021. Prediction of protein-protein interaction network in malaria biomarkers and implication as therapeutic target. Malaysian J Biochem Mol Biol 24: 61-67.
• HELT A-M & GALLOWAY DA. 2001. Destabilization of the retinoblastoma tumor suppressor by human papillomavirus type 16 E7 is not sufficient to overcome cell cycle arrest in human keratinocytes. J Virol 75: 6737-6747.
• HIDAYATULLAH A, PUTRA WE, RIFA’I M, SUSTIPRIJATNO, WIDIASTUTI D, HEIKAL MF, SUSANTO H & SALMA WO. 2022. 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. Karbala Int J Mod Sci 8: 531-542.
• HIDAYATULLAH A, PUTRA WE, SUSTIPRIJATNO, WIDIASTUTI D, SALMA WO & HEIKAL MF. 2023. Molecular docking and molecular dynamics simulation-based identification of natural inhibitors against druggable human papilloma virus type 16 target. Trends Sci 20: 1-12.
• HUH WK ET AL. 2017. 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 Lond Engl 390: 2143-2159.
• JONES DL, ALANI RM & MÜNGER K. 1997. The human papillomavirus E7 oncoprotein can uncouple cellular differentiation and proliferation in human keratinocytes by abrogating p21Cip1-mediated inhibition of cdk2. Genes Dev 11: 2101-2111.
• KAMARAJ B & PUROHIT R. 2013. In silico screening and molecular dynamics simulation of disease-associated nsSNP in TYRP1 gene and its structural consequences in OCA3. BioMed Res Int 2013: e697051.
• KNAPP AA, MCMANUS PM, BOCKSTALL K & MOROIANU J. 2009. Identification of the nuclear localization and export signals of high risk HPV16 E7 oncoprotein. Virology 383: 60-68.
• KUHN B, MOHR P & STAHL M. 2010. Intramolecular hydrogen bonding in medicinal chemistry. J Med Chem 53: 2601-2611.
• LEHOUX M, D’ABRAMO CM & ARCHAMBAULT J. 2009. Molecular mechanisms of human papillomavirus-induced carcinogenesis. Public Health Genomics 12: 268-280.
• LIU X, CLEMENTS A, ZHAO K & MARMORSTEIN R. 2006. Structure of the human papillomavirus E7 oncoprotein and its mechanism for inactivation of the retinoblastoma tumor suppressor. J Biol Chem 281: 578-586.
• LONGWORTH MS & LAIMINS LA. 2004. The binding of histone deacetylases and the integrity of zinc finger-like motifs of the E7 protein are essential for the life cycle of human papillomavirus type 31. J Virol 78: 3533-3541.
• LOURENÇO DE FREITAS N, DEBERALDINI MG, GOMES D, PAVAN AR, SOUSA Â, DOS SANTOS JL & SOARES CP. 2021. Histone deacetylase inhibitors as therapeutic interventions on cervical cancer induced by human papillomavirus. Front Cell Dev Biol 8: 592868.
• MALLAMACE D, FAZIO E, MALLAMACE F & CORSARO C. 2018. The role of hydrogen bonding in the folding/unfolding process of hydrated lysozyme: A review of recent NMR and FTIR results. Int J Mol Sci 19: 3825.
• MCLAUGHLIN-DRUBIN ME & MÜNGER K. 2009. The human papillomavirus E7 oncoprotein. Virology 384: 335-344.
• MOK SC, WONG K-K, LU KH, MUNGER K & NAGYMANYOKI Z. 2018. Chapter 23 - molecular basis of gynecologic diseases. In: Coleman WB & Tsongalis GJ (Eds), Molecular Pathology (Second Edition), Academic Press, p. 507-529.
• MÜNGER K, BASILE JR, DUENSING S, EICHTEN A, GONZALEZ SL, GRACE M & ZACNY VL. 2001. Biological activities and molecular targets of the human papillomavirus E7 oncoprotein. Oncogene 20: 7888-7898.
• NOGUEIRA MO, HOŠEK T, CALÇADA EO, CASTIGLIA F, MASSIMI P, BANKS L, FELLI IC & PIERATTELLI R. 2017. Monitoring HPV-16 E7 phosphorylation events. Virology 503: 70-75.
• OHLENSCHLÄGER O ET AL. 2006. Solution structure of the partially folded high-risk human papilloma virus 45 oncoprotein E7. Oncogene 25: 5953-5959.
• PAL A & KUNDU R. 2020. Human Papillomavirus E6 and E7: The cervical cancer hallmarks and targets for therapy. Front Microbiol 10: 3116.
• PANG CL & THIERRY F. 2013. Human papillomavirus proteins as prospective therapeutic targets. Microb Pathog 58: 55-65.
• PARVEEN A, AKASH MSH, REHMAN K & KYUNN WW. 2016. Dual role of p21 in the progression of cancer and its treatment. Crit Rev Eukaryot Gene Expr 26: 49-62.
• PETUKHOV M, RYCHKOV G, FIRSOV L & SERRANO L. 2004. H-bonding in protein hydration revisited. Protein Sci Publ Protein Soc 13: 2120-2129.
• PHELPS WC, MÜNGER K, YEE CL, BARNES JA & HOWLEY PM. 1992. Structure-function analysis of the human papillomavirus type 16 E7 oncoprotein. J Virol 66: 2418-2427.
• PHILLIPS AC & VOUSDEN KH. 1997. Analysis of the interaction between human papillomavirus type 16 E7 and the TATA-binding protein, TBP. J Gen Virol 78: 905-909.
• PUTRA WE, SUSTIPRIJATNO, HIDAYATULLAH A, HEIKAL MF, WIDIASTUTI D & ISNANTO H. 2023. Analysis of three non-structural proteins, Nsp1, Nsp2, and Nsp10 of Sars-Cov-2 as pivotal target proteins for computational drug screening. J Microbiol Biotechnol Food Sci 12: 1-6.
• SALIM Z & MUNADI E. 2017. Info komoditi tanaman obat, Jakarta: Badan Pengkajian dan Pengembangan Perdagangan Kementerian Perdagangan Republik Indonesia.
• SHIN M-K, BALSITIS S, BRAKE T & LAMBERT PF. 2009. Human papillomavirus E7 oncoprotein overrides the tumor suppressor activity of p21 Cip1 in cervical carcinogenesis. Cancer Res 69: 5656-5663.
• SIVARAMAKRISHNAN M, KANDASWAMY K, NATESAN S, DEVARAJAN RD, RAMAKRISHNAN SG & KOTHANDAN R. 2020. Molecular docking and dynamics studies on plasmepsin V of malarial parasite Plasmodium vivax. Inform Med Unlocked 19: 100331.
• SNEHA P & GEORGE PRIYA DOSS C. 2016. Chapter seven - molecular dynamics: New frontier in personalized medicine. In: Donev R (Ed), Advances in Protein Chemistry and Structural Biology, Academic Press, p. 181-224.
• SYAFITRI RW, FIBRIANI A & ADITAMA R. 2021. Selection of indonesian medicinal plant active compounds as inhibitor candidates of oncoproteins E6 and E7 human papillomavirus type 16 by molecular docking. 3BIO J Biol Sci Technol Manag 3: 9-17.
• TAN S, DE VRIES GEE, VAN DER ZEE GJA & DE JONG S. 2012. Anticancer drugs aimed at E6 and E7 activity in HPV-positive cervical cancer. Curr Cancer Drug Targets 12: 170-184.
• TODOROVIC B, MASSIMI P, HUNG K, SHAW GS, BANKS L & MYMRYK JS. 2011. Systematic analysis of the amino acid residues of human papillomavirus type 16 E7 conserved region 3 involved in dimerization and transformation. J Virol 85: 10048-10057.
• TOMAIĆ V. 2016. Functional roles of E6 and E7 oncoproteins in HPV-induced malignancies at diverse anatomical sites. Cancers 8: 95.
• TOOTS M, USTAV M, MÄNNIK A, MUMM K, TÄMM K, TAMM T, USTAV E & USTAV M. 2017. Identification of several high-risk HPV inhibitors and drug targets with a novel high-throughput screening assay Lambert PF (Ed), PLOS Pathog 13: e1006168.
• VALDOVINOS-TORRES H, OROZCO-MORALES M, PEDROZA-SAAVEDRA A, PADILLA-NORIEGA L, ESQUIVEL-GUADARRAMA F & GUTIERREZ-XICOTENCATL L. 2008. Different isoforms of HPV-16 E7 protein are present in cytoplasm and nucleus. Open Virol J 2: 15-23.
• VALENTE RPP, SOUZA RC, DE MEDEIROS MUNIZ G, FERREIRA JEV, DE MIRANDA RM, LIMA AHL & VIANEZ JUNIOR JLSG. 2020. 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 10: 9625.
• WIDIASTUTI D, WARNASIH S, SINAGA SE, PUJIYAWATI E, SUPRIATNO & PUTRA WE. 2023. Identification of active compounds from Averrhoa bilimbi L. (Belimbing Wuluh) flower using LC-MS and antidiabetic activity test using in vitro and in silico approaches. Trends Sci 20; 1-9.
• ZHANG D & LAZIM R. 2017. Application of conventional molecular dynamics simulation in evaluating the stability of apomyoglobin in urea solution. Sci Rep 7: 44651.
• ZHANG S-Y, ZHOU B-J & WANG Y. 2002. Determination of L-sesamin and L-asarinin in Zanthoxylum(Roxb.) DC. by high performance liquid chromatography. 1 Jun Yi Xue Xue Bao Acad J First Med Coll PLA 22: 654-655.
• ZHU J, LV Y, HAN X, XU D & HAN W. 2017. Understanding the differences of the ligand binding/unbinding pathways between phosphorylated and non-phosphorylated ARH1 using molecular dynamics simulations. Sci Rep 7: 12439.

Publication Dates

• Publication in this collection
  17 July 2023
• Date of issue
  2023

History

• Received
  27 July 2022
• Accepted
  23 Oct 2022