Structure-Based Discovery of Potential HPV E6 and EBNA1 Inhibitors: Implications for Cervical Cancer Treatment

: Cervical cancer is the fourth most diagnosed cancer and the fourth leading cause of cancer death in women globally. Its onset and progression have been attributed to high-risk human papillomavirus (HPV) types, especially 16 and 18, while the Epstein–Barr virus (EBV) is believed to also significantly contribute to cervical cancer growth. The E6 protein associated with high-risk HPV strains, such as HPV16 and HPV18, is known for its role in promoting cervical cancer and other anogenital cancers. E6 proteins contribute to the malignant transformation of infected cells by targeting and degrading tumor suppressor proteins, especially p53. On the other hand, EBV nuclear antigen 1 (EBNA1) plays a crucial role in the maintenance and replication of the EBV genome in infected cells. EBNA1 is believed to increase HPV E6 and E7 levels, as well as c-MYC, and BIRC5 cellular genes in the HeLa cell line, implying that HPV/EBV co-infection accelerates cervical cancer onset and growth. Thus, the E6 and EBNA1 antigens of HPV and EBV, respectively, are attractive targets for cervical cancer immunotherapy. This study, therefore, virtually screened for potential drug candidates with good binding affinity to all three oncoviral proteins, HPV16 E6, HPV18 E6, and EBNA1. The compounds were further subjected to ADMET profiling, biological activity predictions, molecular dynamics (MD) simulations, and molecular mechanics Poisson–Boltzmann surface area (MM/PBSA) calculations. A


Introduction
Several cancers have been linked to various viruses, including human papillomaviruses (HPVs), hepatitis viruses (HBV and HCV), Epstein-Barr virus (EBV), Kaposi's sarcoma herpesvirus (KSHV), human T cell leukemia virus (HTLV-I), and Merkel cell polyomavirus (MCPyV) [1].Virus-associated cancers account for about 10 to 15% of cancer cases globally [2,3], representing over 1.4 million cases annually [1].These two human oncoviruses, HPV and EBV, are responsible for approximately 38% of all tumor-associated cancers [3].Cervical cancer is one such cancer and is primarily associated with persistent infection by high-risk HPV types, emphasizing the pivotal role of viral oncogenesis in its development [4][5][6].Although the disease is preventable and treatable when detected early, challenges persist in ensuring widespread access to preventive measures, which include HPV vaccination and regular screening programs [7][8][9][10].Globally, approximately 604,000 women are diagnosed with cervical cancer, while over 342,000 die annually [11].The death rate is estimated to reach over 440,000 by 2030 [12].Cervical cancer ranks as the fourth most commonly diagnosed cancer and the fourth primary contributor to cancer-related fatalities among women [11,12].
HPVs typically target layered epithelial tissues and are directly linked to various cancers affecting the anogenital region, such as cervical, anal, vulvar, and vaginal carcinomas [1].Additionally, they are associated with carcinomas in other mucosal epithelial areas, particularly in the oropharyngeal region [1].HPV is diagnosed in more than 90% of cervical cancer cases [12], with the most common HPV types being 16 and 18 (approximately 70%) [13].HPV is a relatively small virus (55 nm in diameter) that possesses a non-enveloped circular double-stranded genomic DNA of approximately 8 kb in size [14][15][16].Another oncovirus, EBV, which is an oncogenic herpesvirus, is associated with a spectrum of diseases, ranging from infectious mononucleosis to various neoplastic conditions [17][18][19].These neoplastic diseases include lymphoproliferative disorders in immunocompromised individuals, Burkitt's lymphoma, Hodgkin's lymphoma, and nasopharyngeal carcinoma [20][21][22].More than 95% of healthy adults globally experience asymptomatic infection with EBV [23,24].EBV possesses a linear double-stranded DNA genome of around 175 kb, which undergoes circularization within latently infected cells [23].Within the EBV genome, there are over 80 open reading frames (ORFs), multiple non-coding RNAs, and 44 mature microRNAs (miRNA) [23,25,26].Latent membrane protein 1 (LMP1) and four EBNA proteins comprising EBNA1, EBNA2, EBNA3A, and EBNA3C are crucial EBV gene products required for B cell transformation [25].There are currently no approved EBV drugs, as the drugs that demonstrated efficacy in replicating EBV had limited clinical successes [24,27].A meta-analysis combining data from 25 articles revealed that women with EBV infection had a fourfold higher occurrence of cervical carcinoma compared to women without EBV [28].Moreover, the study indicated that EBV/HPV co-infection is associated with an elevated risk of cervical cancer [28].
Oncoproteins of EBV (LMP1, LMP2A, and EBNA1) and high-risk HPVs (E5, E6, and E7) have been documented to contribute to the onset and growth of various human cancers, including head and neck, cervical, and colorectal [29].In HPV-positive cervical, oral, and lung cancer cells, the upregulation of Foxhead box M1 (FOXM1) expression occurs through E6-mediated NKX2-1 [30,31].The induction of FOXM1 by E6 is orchestrated through the MZF1/NKX2-1 axis, and this regulatory mechanism is responsible for HPV-mediated traits such as soft agar growth, invasiveness, and stemness [30].The activation of the Wnt/β-catenin signaling pathway is implicated in mediating these effects, underscoring the intricate molecular interactions contributing to the oncogenic potential of HPV in these cancer types [30].Also, the interaction between E6 and the cellular ubiquitin ligase E6AP facilitates the degradation of p53, a tumor suppressor protein [32,33].Consequently, a crucial strategy for hindering the survival and proliferation of infected cells involves suppressing the formation of the E6-E6AP complex [34,35].
EBNA1, a nuclear antigen, is crucial for maintaining, replicating, and segregating the EBV genome and is required for preserving the episomal DNA of EBV within infected cells [27,36].It is the sole protein necessary for the replication of the latent virus [36].However, the replication and persistence of the EBV episomal genome rely significantly on EBNA1 binding to the EBV origin of plasmid replication (oriP) element [37,38].EBNA1 is believed to elevate HPV18 E6 and E7 levels along with the cellular genes c-MYC and BIRC5 in the HeLa cell line [39].This implies that co-infection of cervical cells with HPV and EBV might accelerate cervical cancer growth and progression [39].Considering the crucial roles of HPV E6 and EBNA1 in the tumorigenesis, proliferation, and progression of cervical cancer cells, it is imperative to identify drug candidates to augment the therapeutic options for combating this cancer, which is characterized by a high mortality rate.
Polypharmacology, which describes the capacity of a single drug to engage with multiple targets or biological pathways within an organism, is increasingly prominent in drug discovery research.The complexity of certain diseases, including cancers, often renders targeting a single pathway inadequate for achieving the desired therapeutic effect [40][41][42].Hence, exploring compounds with polypharmacological properties has emerged as a promising approach to addressing the multifaceted mechanisms underlying such diseases and enhancing therapeutic efficacy [42].A recent study conducted at the oncology center of the University of Gondar Comprehensive Specialized Hospital (Ethiopia) analyzed data from 124 cervical cancer patients to identify common medication-related problems (MRPs) [43].The study included all patients diagnosed with cervical cancer between 1 January 2016 and 31 December 2020.The study revealed that the most common MRPs included suboptimal dosage, requirement for supplementary drug therapy, and adverse drug reactions (ADRs), accounting for 24.4%, 22.6%, and 22% of cases, respectively [43].Furthermore, 16.1% (27 patients) experienced minor to moderate drug-drug interactions (DDIs), with 48% of these classified as moderate, necessitating close monitoring of the outcome, and 8% categorized as serious reactions, prompting the use of alternative medications in the treatment regimen [43].Interestingly, the study also found that patients with polypharmacy (≥5 drugs) were seven times more likely to experience MRPs compared to their counterparts [43].Polypharmacology presents an opportunity to mitigate drug-drug interactions and enhance patient compliance through simplified dosing regimens [44] while potentially reducing the risk of drug resistance [45,46].
Natural products provide a large and diverse source for finding these treatment options [47].Harnessing the potential of natural products may lead to the discovery of novel compounds with novel mechanisms of action against cervical cancer, contributing to the development of more effective and targeted therapeutic interventions [48,49].Additionally, the exploration of natural products aligns with the growing interest in sustainable and environmentally friendly approaches to drug discovery.Computational-based drug discovery approaches offer an efficient and cost-effective avenue to expeditiously identify potential drug candidates [50,51].Using a robust computational drug discovery pipeline, the enrichment factor for shortlisting promising drug candidates can be significantly enhanced.The integration of computational techniques, such as molecular docking, molecular dynamics simulations, biological activity prediction, and toxicity predictions enables the screening of large compound libraries, shortlisting compounds with the highest likelihood of success.
This study therefore employed computational techniques to sift through a vast traditional Chinese medicine (TCM) chemical landscape [52], seeking compounds that can potentially inhibit the E6 proteins of both HPV16 and HPV18.These compounds, displaying commendable binding affinity to both E6 proteins (docking scores ≤ −8.5 kcal/mol) were re-docked against the EBNA1 protein, unveiling potential polypharmacologic drugs that extend beyond a single viral or cancer target.The shortlisted candidates underwent further rigorous scrutiny through ADMET screening and profiling and biological activity predictions prior to molecular dynamics (MD) simulations and molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) calculations.

Materials and Methods
This study employed a molecular docking technique to virtually screen natural products from the TCM database against three oncoproteins, EBNA1, HPV16 E6, and HPV18 E6 (Figure 1).The top compounds were subjected to safety and toxicity profiling using SwissADME and DataWarrior 5.0.Compounds that were predicted to be drug-like, non-tumorigenic, non-mutagenic, and non-irritants and had no reproductive effects were shortlisted, and their biological activities were predicted.Ten top candidates in complex with each of the proteins were further subjected to 100 ns MD simulations and MM/PBSA calculations (Figure 1).Similarity searches were then performed for the topmost compounds that were predicted to have strong binding affinity to all three oncoproteins.

Materials and Methods
This study employed a molecular docking technique to virtually screen natural products from the TCM database against three oncoproteins, EBNA1, HPV16 E6, and HPV18 E6 (Figure 1).The top compounds were subjected to safety and toxicity profiling using SwissADME and DataWarrior 5.0.Compounds that were predicted to be drug-like, nontumorigenic, non-mutagenic, and non-irritants and had no reproductive effects were shortlisted, and their biological activities were predicted.Ten top candidates in complex with each of the proteins were further subjected to 100 ns MD simulations and MM/PBSA calculations (Figure 1).Similarity searches were then performed for the topmost compounds that were predicted to have strong binding affinity to all three oncoproteins.

Figure 1.
Step-by-step procedure employed in this study to identify potential small molecule inhibitors of HPV E6 and EBNA1.Compounds from traditional Chinese medicine were screened against EBNA1 and E6 proteins of HPV16 and HPV18.Top compounds were subjected to ADMET profiling and biological activity predictions.The shortlisted compounds were further subjected to molecular dynamics simulations and MM/PBSA computations to evaluate their binding affinities.

Obtaining Protein Structures
Structures of HPV16 E6, HPV18 E6, and EBNA1 were searched and retrieved from the Research Collaboratory for Structural Bioinformatics (RCSB) protein databank (PDB) [53,54].The RCSB PDB is a repository that houses experimentally determined protein structures [53,54].In selecting the protein structures to be used in this study, factors such as resolution, the count of missing residues, and the year of structure determination were considered.The lower the resolution value of the protein structure, the higher the level of detail (or resolution) and accuracy in depicting the atomic arrangement within the protein [55].

Protein and Ligand Preparation
The 3D structure data file (sdf) format of two E6 inhibitors, gossypetin and baicalein, were obtained from PubChem [56][57][58] with corresponding compound identifiers (CIDs) Step-by-step procedure employed in this study to identify potential small molecule inhibitors of HPV E6 and EBNA1.Compounds from traditional Chinese medicine were screened against EBNA1 and E6 proteins of HPV16 and HPV18.Top compounds were subjected to ADMET profiling and biological activity predictions.The shortlisted compounds were further subjected to molecular dynamics simulations and MM/PBSA computations to evaluate their binding affinities.

Obtaining Protein Structures
Structures of HPV16 E6, HPV18 E6, and EBNA1 were searched and retrieved from the Research Collaboratory for Structural Bioinformatics (RCSB) protein databank (PDB) [53,54].The RCSB PDB is a repository that houses experimentally determined protein structures [53,54].In selecting the protein structures to be used in this study, factors such as resolution, the count of missing residues, and the year of structure determination were considered.The lower the resolution value of the protein structure, the higher the level of detail (or resolution) and accuracy in depicting the atomic arrangement within the protein [55].

Protein and Ligand Preparation
The 3D structure data file (sdf) format of two E6 inhibitors, gossypetin and baicalein, were obtained from PubChem [56][57][58] with corresponding compound identifiers (CIDs) of 5280647 and 5281605.Baicalein and gossypetin have been shown to inhibit E6 (of both HPV16 and HPV18), limit E6 binding to E6AP, and prevent E6-mediated p53 degradation [59].The 2D sdf of EBNA1 inhibitor VK-1727 [60] was also obtained from PubChem and then converted to 3D using OpenBabel [61] due to the unavailability of a 3D sdf format.The screening library consisted of 35,161 compounds from the traditional Chinese Medicine (TCM) catalogue found on ZINC15 database [62,63].The TCM compounds were pre-filtered using a molecular weight criterion of 150 to 600 g/mol, as previously performed [64].Also, duplicates were removed using OSIRIS DataWarrior (Idorsia Pharmaceuticals Ltd., version 5.0.0,Allschwil, Switzerland) [65].After filtering, 25,212 TCM compounds were obtained and used for the molecular docking study.All the ligands were then subjected to energy minimization using the Universal force field (UFF) [66,67] and were converted to protein data bank, partial charge (Q), and atom type (T) ("pdbqt") formats via PyRx 0.9 [68,69].The three protein structures (HPV16 E6, HPV18 E6, and EBNA1) were also subjected to energy minimization via GROMACS version 5.1.5using the optimized potentials for liquid simulations (OPLS) force field prior to docking.After energy minimization, each protein structure was converted to pdbqt format via PyRx using the "Make Macromolecule" option.

Molecular Docking Studies (E6 Proteins of HPV16 and HPV18)
Molecular docking, a computational structure-based approach in drug discovery, aids in the identification of potential therapeutic compounds by predicting interactions between proteins and ligands [70][71][72].A major advantage of molecular docking is its capacity to perform these tasks without the need for prior knowledge of the chemical structure of other target modulators [71].This study employed a molecular docking tool, AutoDock Vina [73] to virtually screen the TCM ligands and the two known E6 inhibitors against both E6 proteins.Compounds that demonstrated docking scores less than −8.4 kcal/mol for both proteins were shortlisted for analysis.For the HPV16 E6 protein, a grid box dimension of 51.88 × 42.19 × 43.44 Å 3 was used centering the protein at x = 39.12Å, y = 38.85Å, and z = 38.90Å.Using a center of x = 38.89Å, y = 39.18Å, and z = 39 Å, a grid box of size of 45.59 × 46.03 × 44.56 Å 3 was also set up for the HPV18 E6 protein.For all the docking runs, the exhaustiveness was set to 8.

Predictive Evaluation of ADMET Properties
The shortlisted compounds that demonstrated good binding affinities to both E6 proteins were subjected to ADMET profiling using SwissADME [74].The compounds were assessed based on Lipinski's rule of 5 [75] and Veber's rule [76].Toxicity risks, including mutagenicity, tumorigenicity, irritancy, and reproductive effects of the shortlisted compounds, were also determined using OSIRIS DataWarrior 5.0 [65].These analyses aimed to provide a more thorough understanding of the compounds' pharmacokinetic and safety profiles, contributing valuable insights into their potential as therapeutic candidates.These analyses were also performed to identify and eliminate potentially toxic candidates or compounds with low safety profiles.

Re-Docking Top Compounds against EBNA1
The shortlisted compounds with insignificant toxicities and good binding affinities to both E6 proteins were further re-docked against the EBNA1 protein using AutoDock Vina.A literature search revealed two pockets of the EBNA1 that are close and may overlap.The first is the DNA binding pocket, which comprises residues Arg469, Gly472, Lys514, Tyr518, Arg521, Pro535, and Leu536 [77,78].The other pocket, which is a deep hydrophobic pocket [79] and is also found in the DNA binding domain (DBP), includes residues Lys477, Asn480, Ile481, Gly484, Ser516, Asn519, Leu520, Leu582, Val583, Thr585, Lys586, Thr590, and Ile593 [77,80].EBNA1 inhibitor H20 was shown to bind in this pocket [77].This region is crucial for DNA binding to EBNA1 [79].The grid box with dimensions of 19.68 × 19.41 × 20.85 Å 3 was set up to cover both pockets, and the protein was centered at x = 35.33Å, y = 40.31Å, and z = 47.11Å prior to docking.

Further Re-Docking of Top Compounds against All 3 Oncoproteins Using DockThor-VS
The top TCM compounds and the known inhibitors were further re-docked against all three proteins using DockThor-VS [81,82] to provide an extra layer of confidence and reliability in the docking scores obtained via AutoDock Vina.For all 3 proteins, the grid size was set at x = 40 Å, y = 40 Å, and z = 40 Å to cover the proteins, while the discretization was 0.42.Also, "virtual screening" was selected as the algorithm precision option.Default parameters were used for number of evaluations (500,000), population size (750), initial seed (−1985), and number of runs (12).The soft docking option was also selected.For EBNA1, the grid center was set at x = 36.1225Å, y = 34.7735Å, and z = 35.9885Å.For HPV16 E6, a grid center of x = 35.7925Å, y = 36.861Å, and z = 38.2585Å was used, while x = 40.9985Å, y = 39.821Å, and z = 42.948Å was specified for HPV18 E6 docking.

Visualizing Protein-Ligand Interaction Profiles
Shortlisted compounds and the known inhibitors in complex with the various proteins were subjected to protein-ligand interaction profiling using LigPlot+ (version v.2.2.8) [83].LigPlot+ generates interaction maps of the query complex and displays the hydrogen bonds and hydrophobic interactions that exist between the protein and the ligand.

Prediction of Biological Activity of Lead Compounds
Prediction of activity spectra for substances (PASS), hosted on the Way2Drug platform, was used to determine the biological activities of the shortlisted compounds [84][85][86].PASS employs a predictive algorithm to predict diverse biological activities associated with compounds based on their chemical structures.The SMILES format of the query is submitted as the input, and the PASS algorithm evaluates the query structure by comparing it to structures in its database.PASS operates on the principle that the biological activity of a compound corresponds to its structure.The predictive accuracy of PASS averages around 95% as per leave-one-out cross-validation estimates [85].

Molecular Dynamics Simulation
The three unbound proteins (EBNA1, HPV16 E6, and HPV18 E6 proteins), as well as the protein-ligand complexes, were subjected to molecular dynamics simulations using OPLS force field via GROMACS version 5.1.5[87,88].The MD simulations were performed to analyze and compare the dynamic behavior, structural stability, and intermolecular interactions of the proteins in their unbound states and when bound to the shortlisted ligands.Trajectory analysis, root mean square deviation (RMSD), root mean square fluctuation (RMSF), the radius of gyration (Rg), and hydrogen bond analyses were performed to gain insights into the conformational changes associated with the protein-ligand interactions.These simulations aimed to provide a comprehensive understanding of the molecular mechanisms governing the binding affinity and stability of the investigated protein-ligand complexes.
For the protein-ligand complexes, prior to the MD simulations, hydrogens were added to the ligands, and their respective charges were computed using UCSF Chimera [89].The ligands were then submitted to LigParGen [90] to generate their respective topology files ("gro" and "itp" files) that are compatible with GROMACS.The 1.14*CM1A charge model was used to generate the topology files [91].Protein-ligand complexes in "gro" format were then prepared using the Sublime Text editor.The complexes were then inspected using PyMOL [92].The unbound proteins and the complexes were then centered in cubic boxes 1 nm away from the edges of the box.Each system was solvated using the TIP4P water model, and chlorine ions were then added to neutralize the charges.Each system was energy-minimized using the steepest descent algorithm through a maximum of 5000 steps or until the maximum force reached <10.0 kJ/mol, whichever occurred first.The systems were then subjected to 100 ps isothermal-isobaric ensemble (NPT) and 100 ps canonical ensemble (NVT) equilibrations before running the 100 ns MD simulations.For the MD runs, the velocity rescaling (V-rescale) thermostat and the Parrinello-Rahman barostat were used.After the simulations, Xmgrace was used to plot the RMSD, RMSF, Rg, and hydrogen bond graphs.

Molecular Mechanics Poisson-Boltzmann Surface Area Calculations
Using g_mmpbsa [93], MMPBSA decomposition was carried out to analyze the energetic contributions of the complexes.The energy terms analyzed include van der Waals (vdW), electrostatic, polar solvation, solvent-accessible surface area (SASA), and enthalpybased estimates of binding energies.The binding free energy calculations help to understand the energetics associated with the protein-ligand interactions.The per-residue energy contributions of each system were also determined to identify hotspot residues that contribute significantly to ligand binding.

Structural Similarity Search of Top Compounds
DrugBank [94,95] was employed to perform a structural similarity search of the shortlisted compounds using their SMILES and chemical structure as input.A similarity threshold of ≥0.7 was used, as this threshold is often considered indicative of substantial structural similarity [96,97].A literature search was conducted to establish the connection between the analogs and their relevance to EBV, HPV, and cervical cancer.

Molecular Docking Analysis
A total of 25,131 TCM compounds were successfully screened against both HPV16 and HPV18 E6 proteins.A total of 595 compounds had docking scores of less than −8.4 kcal/mol for both E6 proteins and were shortlisted for further analysis.

Docking against HPV16 E6
The known E6 inhibitors, baicalein and gossypetin, had docking scores of −7.9 and −6.8 kcal/mol, respectively.After virtual screening, a total of 801 natural product compounds had docking scores of −8.5 kcal/mol and lower with the HPV16 E6 and were thus shortlisted for further analysis.Compound ZINC000070454124 demonstrated the lowest docking score of −10.8 kcal/mol (Table 1), implying the strongest binding affinity with E6, followed by ZINC000103580868 and ZINC000085543515, both with a docking score of −10.4 kcal/mol.Compounds ZINC000070451100 and ZINC000103579839 both had a docking score of −10.3 kcal/mol, while ZINC000103526876 also had a good docking score of −10.2 kcal/mol.ZINC000004097766 had a docking score of −10.1, while ZINC000095912674, ZINC000085631224, and ZINC000103527856 had a docking score of −10 kcal/mol.A previous study docked 20 natural compounds with known antiviral and anticancer properties against the HPV16 E6 protein using Autodock 4.2.6 [98].The top 20 compounds had docking scores ranging from −4.33 to −8.46 kcal/mol, with ginkgetin demonstrating the strongest affinity [98].Ginkgetin had docking scores of −8.45, −8.46, and −7.22 kcal/mol when docked against the EA6P, p53, and Myc binding sites [98].These binding sites have been shown to overlap [99].Another study also docked four compounds from marine sources against HPV16 E6 using Autodock 4.2 [100].The four compounds, salicylihalamide B, salicylihalamide A, frigocyclinone, and ascorbic acid, had docking scores of −8.92, −8.76, −8.01, and −3.82 kcal/mol, respectively [100].Herein, the shortlisted compounds demonstrated better docking scores than those previously obtained, suggesting improved or better affinity (Table 1).Baicalein and gossypetin, which are known inhibitors of E6, had docking scores of −6.8 and −6.5 kcal/mol when screened against the E6 protein of HPV18.Compound ZINC000103584225 from the TCM library demonstrated the highest binding affinity to the HPV18 E6 protein with a low docking score of −11.4 kcal/mol, followed by ZINC000085593983 with −11.2 kcal/mol and ZINC000095909496 and ZINC000095911538, both with a docking score of −11.1 kcal/mol.ZINC000095913861 also had a docking score of −10.9 kcal/mol, while ZINC000103574367, ZINC000103554058, ZINC000095911347, and ZINC000103526876 had a docking score of −10.8 kcal/mol.Also, compounds ZINC000085530466, ZINC000070454124, ZINC000067902771, and ZINC000085541858 all had a docking score of −10.7 kcal/mol.In all, a total of 3063 compounds had docking scores of −8.5 kcal/mol and lower when screened against the HPV18 E6 protein.A previous study docked 12 natural products with known anti-HPV activity against the HPV18 E6 protein using Autodock 4.2, and their docking scores ranged from −2.86 to −5.85 kcal/mol [101].Another study also docked known antivirals (natural products) against the HPV18 E6 protein using Autodock 4.2, and the docking scores of the top 19 compounds ranged from −7.48 to −14.77 kcal/mol [102].Also, doxorubicin and pinoresinol monomethyl etherβ-D-glucoside (PMG), isolated from the dichloromethane extract of Jurinea macrocephala subsp.elbursensis, had docking scores of −3.08 and −4.75 when docked against the HPV18 E6 protein [103].Doxorubicin and PMG demonstrated cytotoxicity against HeLa cells with average IC 50 values of 1.33 and 10.13 µg/mL (after 24 h exposure) and 0.136 and 3.54 µg/mL (after 48 h exposure), respectively [103].However, doxorubicin and PMG are not known E6 inhibitors.The docking scores of the compounds obtained in this study imply better binding affinity to the E6, warranting further experimental testing.

Predictive Evaluation of Drug-Like Properties: ADMET Profiling
The top compounds were subjected to ADMET prediction using SwissADME (Table 2) [74].All ADME properties were predicted for all compounds except for ZINC000085571466.A total of 418 and 176 compounds passed and failed Lipinski's rule, respectively.A total of 436 and 158 compounds passed and failed Veber's rule, respectively.Veber's rule requires a drug-like or orally bioavailable compound to have a topological surface area (TPSA) less than 140 Å 2 and fewer than 10 rotatable bonds [76].Lipinski's rule, on the other hand, requires a drug-like molecule to not violate more than one of the following criteria: molecular weight should be less than 500 Daltons; number of hydrogen bond acceptor groups (usually nitrogen or oxygen atoms) should be less than or equal to 10; number of hydrogen bond donor groups (usually OH or NH groups) should be less than or equal to five; and the octanol-water partition coefficient (LogP) should be less than five [75].A total of 82 compounds failed Lipinski's rule only, 64 failed Veber's rule only, and 94 compounds failed both rules simultaneously.In all, a total of 354 compounds passed both Lipinski's and Veber's rules and were further subjected to toxicity risk prediction using OSIRIS DataWarrior 5.0 [65].A total of 29, 29, and 296 compounds were predicted to have low, high, and no mutagenic properties, respectively.For tumorigenicity, 32 compounds were predicted to be low, 32 were predicted to be high, and 290 had no tumorigenic properties.A total of 25 compounds were predicted to possess low reproductive effects, 68 had high reproductive effects, and 261 were predicted to have no reproductive effects (Figure 2).For the irritant toxicity risk prediction, 20, 63, and 271 compounds had low, high, and no irritancy, respectively (Figure 2).Altogether, 166 compounds were predicted to exhibit characteristics of non-mutagenicity, non-tumorigenicity, absence of irritant properties, and no adverse reproductive effects.These compounds were consequently chosen as reasonably safe and used for further analysis.Bar plot illustrating the count of compounds categorized as "high", "low", or "none" for each of the four toxicity risks, including mutagenicity, tumorigenicity, reproductive effect, and irritant.

Re-Docking Top Compounds against EBNA1
The top 166 TCM compounds and the known inhibitors were re-docked against the EBNA1 protein.Only compounds with docking scores ≤−7.0 kcal/mol were considered for further studies, as this AutoDock Vina threshold has been shown to effectively separate good binders from weak binders of viral proteins [104].A total of 26 compounds were obtained as potentially good binders of E6 (both HPV16 and HPV18) and EBNA1 proteins.Compound ZINC000033831887 demonstrated the highest binding affinity, with a docking score of −7.8 kcal/mol, followed by ZINC000085597263 and ZINC000085543579, both with a docking score of −7.7 kcal/mol (Table 1).Compounds ZINC000000899994, ZINC000085628525, ZINC000070454501, ZINC000085950180, and ZINC000095914132 all .Bar plot illustrating the count of compounds categorized as "high", "low", or "none" for each of the four toxicity risks, including mutagenicity, tumorigenicity, reproductive effect, and irritant.

Further Re-Docking of Top Compounds against All Three Oncoproteins Using DockThor
The top 26 TCM compounds, along with the three known inhibitors, were re-docked against all three proteins using a different docking tool, DockThor-VS [81,82].This step aimed to bolster confidence and reliability in the shortlisted compounds and docking outcomes.All three known inhibitors were successfully screened against the proteins; however, for the TCM compounds, all but ZINC000095909247, ZINC000095913339, ZINC000095914132, and ZINC000070454552 were screened successfully (Table S1).A previous study that docked drugs against various SARS-CoV-2 targets reported that DockThor-VS docking scores higher than −6.8 kcal/mol indicate binding affinities above approximately 10 µM, which were considered as low binding affinities [82].Compounds with docking scores between −6.8 and −8.2 kcal/mol were considered to have moderate binding affinities, while those with docking scores below −8.2 kcal/mol demonstrated high binding with submicromolar affinities [82].EBNA1 inhibitor, VK-1727 demonstrated high binding affinity with all three proteins with docking scores of −8.702, −8.301, and −9.363 against EBNA1, HPV16 E6, and HPV18 E6, respectively (Table S1).Further studies are required to explore the potential repurposing of VK-1727 as an effective HPV E6 inhibitor.Similarly, compounds ZINC000014588133, ZINC000085568150, ZINC000070454501, ZINC000085597263, ZINC000085568136, ZINC000085543579, and ZINC000014768164 demonstrated high binding affinities with all three proteins (Table S1).Baicalein demonstrated moderate binding affinities with all three proteins, having docking scores of −7.178, −7.463, and −7.43 kcal/mol against EBNA1, HPV16 E6, and HPV18 E6 proteins, respectively.Gossypetin, on the other hand, demonstrated moderate docking scores against EBNA1 (−6.984 kcal/mol) and HPV18 E6 (−7.28 kcal/mol) while exhibiting a low binding affinity for HPV16 E6 (−6.751 kcal/mol) [Table S1].For both HPV E6 proteins, all the TCM compounds demonstrated higher binding affinities than the two known E6 inhibitors, baicalein and gossypetin.The results obtained with DockThor-VS are consistent with those from AutoDock Vina.

Biological Activity Prediction of Selected Lead Compounds
The biological activities of the top 26 compounds obtained were predicted using the prediction of activity spectra for substances (PASS) [85,86].Biological activity predictions add depth to the previously obtained results, corroborating the potential inhibitory activities of the most promising candidates.PASS computes a compound's probability to be active (Pa) and probability to be inactive (Pi) using the activity data of structurally similar molecules [84,85,107].Only activities with Pa > 0.3 and Pa > Pi were considered for each compound.The biological activities of three compounds (ZINC000085568150, ZINC000085568136, and ZINC000103533241) could not be predicted using PASS.Most of the compounds were predicted to be antivirals.Compounds ZINC000003588336, ZINC000014768164, and ZINC000013380012 were predicted to possess antineoplastic activity specifically against cervical cancer, with Pa values of 0.311, 0.306, and 0.302 and Pi

Biological Activity Prediction of Selected Lead Compounds
The biological activities of the top 26 compounds obtained were predicted using the prediction of activity spectra for substances (PASS) [85,86].Biological activity predictions add depth to the previously obtained results, corroborating the potential inhibitory activities of the most promising candidates.PASS computes a compound's probability to be active (Pa) and probability to be inactive (Pi) using the activity data of structurally similar molecules [84,85,107].Only activities with Pa > 0.3 and Pa > Pi were considered for each compound.The biological activities of three compounds (ZINC000085568150, ZINC000085568136, and ZINC000103533241) could not be predicted using PASS.Most of the compounds were predicted to be antivirals.Compounds ZINC000003588336, ZINC000014768164, and ZINC000013380012 were predicted to possess antineoplastic activity specifically against cervical cancer, with Pa values of 0.311, 0.306, and 0.302 and Pi values of 0.017, 0.018, and 0.018, respectively.
Ten compounds, including ZINC000070454124, ZINC000085628525, ZINC000085597263, ZINC000095909247, ZINC000003588336, ZINC000070454552, ZINC000085597267, ZINC000033831887, ZINC000085950180, and ZINC000014649947, were also predicted to be antineoplastics for treating non-Hodgkin's lymphoma (NHL) (with Pa > Pi and Pa > 0.3).EBV is known to cause mononucleosis and is implicated in certain types of NHL, particularly in immunocompromised individuals [109].EBV is most commonly linked to the development of Burkitt lymphoma and post-transplant lymphoproliferative disorder (PTLD) [110].Burkitt lymphoma, in particular, has a well-established association with EBV, especially in endemic regions [22,109].In immunocompromised individuals, EBV can lead to the development of lymphomas, including diffuse large B cell lymphoma (DLBCL) [109,111].
Seventeen compounds were predicted to be apoptosis agonists, while thirteen were predicted to be Janus kinase 2 (JAK2) expression inhibitors.The activation of the JAK/STAT pathway has been proposed to induce the expression of antiapoptotic genes, potentially contributing significantly to the progression of cervical cancer tumorigenesis [115][116][117][118]. Four compounds (ZINC000085597263, ZINC000095909247, ZINC000014768164, and ZINC000085597267) were predicted to be transcription factor NF kappa B (NF-κB) inhibitors.The activation of the NF-κB signaling pathway has been implicated in cervical cancer cell proliferation, invasion, and metastasis [119,120].The HPV E6 and E7 proteins have been shown to be involved in the activation of NF-κB [119][120][121][122][123][124].

MD Simulation Analysis
Molecular dynamics (MD) simulations have been extensively used in computational drug discovery projects to investigate the conformational transitions of proteins resulting from the binding and unbinding of ligands [125].Using MD simulations, the dynamic behaviors of biomolecules at the atomic level can be explored, providing insights into the structural changes associated with ligand interactions.Herein, 100 ns were performed for the apo and holo forms of all three proteins.Average RMSD, RMSF, Rg, and number of hydrogen bonds formed were calculated for each system throughout the 100 ns MD runs (Table 4).A significant constraint is the absence of simulations conducted under the typical physiological salt concentration of ~0.15 M. Future studies could explore conducting simulations under both physiological salt concentrations and those found in cancer cells, which often tend to be higher than physiological levels [126].To evaluate convergence, a clustering analysis (gmx cluster) was employed to generate clusters from the first 80 ns of each simulation and the last 20 ns using a similarity cutoff of 0.2 nm and a time interval of 1 ns (Table S2).A small cutoff value was used to obtain well-defined clusters with structurally similar members.For each system, a total of 81 (0-80 ns) and 20 (81-100 ns) frames were generated from the first 80 ns and the last 20 ns segments, respectively.For the first 80 ns of each simulation, several distinct clusters were identified and denoted as first1, first2, first3, etc.Similarly, for the last 20 ns, clusters labeled last1, last2, last3, etc., were generated.The clusters are labeled according to decreasing cluster population size.In other words, cluster1 (first1 or last1) has the highest number of members followed by cluster2 (first2 or last2), and subsequent clusters follow suit in decreasing order.To assess convergence, the cluster representative structures obtained from the first 80 ns of the simulation were compared to those from the last 20 ns (Table S3).Aligning the C-alpha atoms of each conformation from the final 20 ns to each conformation from the first 80 ns provides a measure of how similar the structures are between the two segments of the simulation.This alignment can reveal significant structural changes or consistency between conformations from the two segments.The structural similarity was quantified using the RMSD metric (Table S3), which measures the average distance between atoms of aligned structures [127,128].By comparing the cluster populations between the first 80 ns and the last 20 ns, any significant shifts in the distribution of conformations can be observed.A simulation with a relatively converged state will exhibit stability or similarity between the populations of the two segments.However, if there are major differences in cluster populations, it may indicate ongoing dynamics or lack of convergence.Lower RMSD values indicate greater similarity between the conformations, while higher RMSD values suggest more significant structural changes.Relatively low and consistent RMSD values (0-3 Å) indicate a stable or converged simulation, whereas high RMSD values (>3 Å) suggest ongoing dynamics or a lack of convergence [129].Lower RMSD values indicate greater structural similarity between the clusters.
For the EBNA1 systems, all the RMSD values were below 1.5 Å (most were below 1 Å) (Table S3), implying no significant differences between the first 80 ns and the last 20 ns.It is also worth mentioning that for all the EBNA1 systems, the RMSD values obtained from aligning the representative structures of the largest clusters (first1 and last1) were below 1 Å except for the EBNA1-ZINC000085631200 complex, which had an RMSD of 1.148 Å (Table S3).These results suggest that all the EBNA1 systems reached a state of convergence.For all the HPV16 and HPV18 E6 systems, varying degrees of similarity were observed between clusters from the last 20 ns and those from the first 80 ns (Table S3).Most aligned clusters exhibited relatively low RMSD values, suggesting a high level of structural conservation throughout the simulation (Table S3).For the HPV16 E6 systems, the RMSD values obtained from the alignment of the representative structures from the largest clusters (first1 and last1) were below 3 Å, indicative of converged states during the simulation.A few RMSDs above 3 Å were observed when representative structures from clusters with smaller population sizes were aligned.This is not surprising as clusters with smaller populations tend to deviate (or are structurally different).For the HPV18 E6 systems, alignments of the representative structures from the largest clusters (first1 and last1) yielded RMSDs less than 2 Å, except for the unbound HPV18 E6 protein, and HPV18 E6-ZINC000085631200, HPV18 E6-ZINC000095913339, and HPV18 E6-ZINC000085543530 complexes (Table S3).Overall, the clustering analyses showed that the systems converged, suggesting the reproducibility of the simulations and that 100 ns is sufficient to study the dynamic behavior of EBNA1, HPV16 E6, and HPV18 E6 systems.

RMSD Analysis
To assess the stability of the protein-ligand complexes, the RMSDs were analyzed over the course of a 100 ns MD simulation.For the EBNA1-ligand systems, the apo-EBNA1 had an average RMSD of 0.36 nm (Figure 6a,b, and Table 4).From the plot, complexes EBNA1-VK1727, EBNA1-ZINC000070454124, EBNA1-ZINC000085568136, EBNA1-ZINC000085631200, EBNA1-ZINC000095913339, and EBNA1-ZINC000085543530 demonstrated lower RMSDs than the apo-EBNA1 (Table 4 and Figure 6a,b), implying they achieved better stability after ligand binding.For the HPV16 E6 systems, E6 complexed with baicalein, ZINC000013380012, ZINC000070454124, ZINC000085631200, ZINC000095913339, ZINC000085628525, and ZINC000085543530 had lower average RMSD values than the apo protein, implying better stability after binding to the E6 protein (Figure S4a,b).A previous MD study on HPV16 E6-salicylihalamide B complex revealed a stable RMSD with an average of 0.35 nm [100].All complexes except ZINC000014588133, ZINC000085568136, ZINC000095909247, and ZINC000085597263 demonstrated lower RMSDs than salicylihalamide B [100].For the HPV18 E6 systems, all the complexes except E6-baicalein achieved better stability (lower RMSDs) than the unbound E6 protein (Table 4 and Figure S5a,b).The average RMSD of the E6-baicalein complex was observed to be 0.46 nm, which was the same as that of the apo-HPV18 E6 protein (Table 4).
Computation 2024, 12, x FOR PEER REVIEW 20 of 40 ZINC000085631200, ZINC000095913339, ZINC000085628525, and ZINC000085543530 had lower average RMSD values than the apo protein, implying better stability after binding to the E6 protein (Figure S4a,b).A previous MD study on HPV16 E6-salicylihalamide B complex revealed a stable RMSD with an average of 0.35 nm [100].All complexes except ZINC000014588133, ZINC000085568136, ZINC000095909247, and ZINC000085597263 demonstrated lower RMSDs than salicylihalamide B [100].For the HPV18 E6 systems, all the complexes except E6-baicalein achieved better stability (lower RMSDs) than the unbound E6 protein (Table 4 and Figure S5a,b).The average RMSD of the E6-baicalein complex was observed to be 0.46 nm, which was the same as that of the apo-HPV18 E6 protein (Table 4).

Rg Analysis
The radius of gyration (Rg) measures the level of compactness in a protein's structure [130].A higher Rg value indicates reduced compactness, while a low Rg value corresponds to a high level of compactness [130].Rg is also used to assess the stability of an unbound protein or protein-ligand complexes during MD simulations, serving as an indicator of whether the structure remains stably folded or unfolds [130,131].Consistency

Rg Analysis
The radius of gyration (Rg) measures the level of compactness in a protein's structure [130].A higher Rg value indicates reduced compactness, while a low Rg value corresponds to a high level of compactness [130].Rg is also used to assess the stability of an unbound protein or protein-ligand complexes during MD simulations, serving as an indicator of whether the structure remains stably folded or unfolds [130,131].Consistency in the Rg values throughout the MD simulation suggests a stable, folded structure, while fluctuations may indicate an unfolded structure or a poorly folded structure [130].
For the MD simulations involving the apo-EBNA1 and holo-EBNA1, the average Rg of all the complexes were comparable to or lower than that of the unbound protein, with the exception of EBNA1-VK1727, EBNA1-ZINC000070454124, EBNA1-ZINC000085631200, and EBNA1-ZINC000085628525 complexes.The Rg of the unbound EBNA1 remained stable (average of 1.6 nm) until about 60 ns, at which a rise to an average of 1.65 nm was observed (Table 4 and Figure 7a,b).From the Rg plots, EBNA1 complexed with ZINC000013380012, ZINC000014588133, ZINC000085597263, ZINC000095913339, ZINC000085628525, and ZINC000085543530 demonstrated lower Rg than the unbound EBNA1 toward the end of the 100 ns simulation period.For the HPV16 E6 systems, the apo E6 had an average Rg of 1.72 nm (Table 4 and Figure S6a,b).All the protein-ligand complexes demonstrated comparable Rg to that of the apo-E6.Previous MD simulation study on HPV16 E6salicylihalamide B complex revealed a stable Rg ranging between 1.65 and 1.75 nm (average of 1.7 nm) [100], consistent with the results obtained in this study.For the HPV18 E6 systems, the unbound protein had an average Rg of 1.84 nm, while the E6-ligand complexes had lower average Rg values throughout the 100 ns simulation period (Table 4 and Figure S7a,b).This implies that upon binding, the known inhibitors and the shortlisted compounds caused the HPV18 E6 to be more compact.study.For the HPV18 E6 systems, the unbound protein had an average Rg of 1.84 nm, while the E6-ligand complexes had lower average Rg values throughout the 100 ns simulation period (Table 4 and Figure S7a,b).This implies that upon binding, the known inhibitors and the shortlisted compounds caused the HPV18 E6 to be more compact.
x FOR PEER REVIEW 22 of 40

Post-MD MM/PBSA Analysis
Aside from the dynamical structural information obtained through MD simulations, energetic information about protein-ligand interactions can also be explored using MD simulations [132].The wealth of data obtained from these simulations plays an important role in unraveling the intricate relationship between structure and function in the target, shedding light on the essence of protein-ligand interactions.Also, MD simulations contribute to enhancing the enrichment factor of virtual screening through the provision of a

Post-MD MM/PBSA Analysis
Aside from the dynamical structural information obtained through MD simulations, energetic information about protein-ligand interactions can also be explored using MD simulations [132].The wealth of data obtained from these simulations plays an important role in unraveling the intricate relationship between structure and function in the target, shedding light on the essence of protein-ligand interactions.Also, MD simulations contribute to enhancing the enrichment factor of virtual screening through the provision of a more precise ranking of hit compounds, coupled with enthalpy-based estimates of binding energy calculations [132].The utilization of the MM/PBSA method subsequent to MD simulations facilitates the quantitative assessment of binding affinity between a protein and a ligand.Herein, the contributing energy terms which include van der Waals (vdW), electrostatic, polar solvation, SASA, and the enthalpy-based estimates of binding free energies of the various protein-ligand complexes were computed using the g_mmpbsa tool [93].For each complex, the energy contributions per residue were also assessed to identify hotspot residues within the binding sites, highlighting their significance in facilitating ligand binding.
Based on the MM/PBSA analyses of the EBNA1, HPV16 E6, and HPV18 E6 protein in complex with the ligands, six compounds, namely, ZINC000013380012, ZINC000070454124, ZINC000014588133, ZINC000085568136, ZINC000095909247, and ZINC000085597263, demonstrated very strong affinity to the three proteins and are worthy of further experimental testing to ascertain their potential anti-cervical cancer activity (Table 5).

Per-Residue Energy Contributions
The per-residue energy contributions of the complexes via MM/PBSA calculations were also assessed (Figures S10-S12).The per-residue energy decomposition helps us understand the energetic aspects of molecular interactions within a complex.This approach breaks down the total binding free energy of a molecular system into contributions from individual residues.Per-residue energy decompositions provide insights into the particular amino acids or residues that play a substantial role in the binding affinity between proteins and ligands and promote the stability of the complex.Previous studies have suggested that residues that contribute energies above 5 or below −5 kJ/mol can be considered critical for ligand binding [133].In this study, for the known inhibitor, VK-1727, no EBNA1 residue contributed above 5 or below −5 kJ/mol to its binding.However, noteworthy contributions were observed for residues Lys477, Ile481, and Leu582, contributing favorably to VK-1727 binding with energies of −3.96, −4.95, and −3.81 kJ/mol, respectively, while Asp581 contributed a positive energy of 3.82 kJ/mol (Figure S10j).Based on the energy contributions in the EBNA1-VK-1727 complex, residues with contributions above 3 or below −3 kJ/mol are identified as critical for ligand binding.

Structural Similarity Search of the Potential Candidates
The literature was searched to identify the existing use of the six top compounds in the treatment of viruses, cancers, and related complications or diseases.Structural similarity searches were also performed via DrugBank to identify analogs of the top six compounds [ZINC000013380012 (broussonol B), ZINC000070454124, ZINC000014588133, ZINC000085568136, ZINC000095909247, and ZINC000085597263] (Figure 10) [94,95].Structural similarity searches help to elucidate potential biological activities by revealing connections between chemical structures and bioactive properties [134][135][136].This approach facilitates the identification of analogs, enabling the recognition of shared pharmacophores and mechanisms of action.The DrugBank-based similarity searches for the top compounds, using a similarity score threshold of 0.7, revealed analogous drugs, with the exception of ZINC000070454124 and ZINC000085597263.Compound ZINC000013380012 (broussonol B) is found in Broussonetia papyrifera, commonly known as paper mulberry and B. kazinoki (kozo) [137].The cytotoxic activity of the methanolic extracts of leaf, bark, and fruit of B. papyrifera was investigated in relation to MCF-7, HeLa, and HepG2 cell lines [138].The leaf extract exhibited significant cytotoxic activity against MCF-7 and HeLa cell lines, with IC50 values of 105 µg/mL and 110 µg/mL, confirmed by MTT assay (IC50 values: 87.5 µg/mL and 106.2 µg/mL) [138].The bark extract demonstrated notable activity against HeLa cells (IC50: 75.3 µg/mL and 88.3 µg/mL), while both leaf and bark extracts displayed moderate activity against HepG2 cells.Conversely, the fruit extract showed insignificant cytotoxicity across all tested cell lines [138].Another study also reported that polyphenols extracted from B. papyrifera stem bark induced apoptosis in HepG2 cells through mitochondria-mediated mechanisms, Compound ZINC000013380012 (broussonol B) is found in Broussonetia papyrifera, commonly known as paper mulberry and B. kazinoki (kozo) [137].The cytotoxic activity of the methanolic extracts of leaf, bark, and fruit of B. papyrifera was investigated in relation to MCF-7, HeLa, and HepG2 cell lines [138].The leaf extract exhibited significant cytotoxic activity against MCF-7 and HeLa cell lines, with IC 50 values of 105 µg/mL and 110 µg/mL, confirmed by MTT assay (IC 50 values: 87.5 µg/mL and 106.2 µg/mL) [138].The bark extract demonstrated notable activity against HeLa cells (IC 50 : 75.3 µg/mL and 88.3 µg/mL), while both leaf and bark extracts displayed moderate activity against HepG2 cells.Conversely, the fruit extract showed insignificant cytotoxicity across all tested cell lines [138].Another study also reported that polyphenols extracted from B. papyrifera stem bark induced apoptosis in HepG2 cells through mitochondria-mediated mechanisms, achieved by the deactivation of ERK and AKT signaling pathways [139].Broussonol B has also been previously shown via MTT assay to possess low cytotoxicity against A549 and HCT-8 human tumor cell lines [137].
Furthermore, a structural similarity search revealed that ZINC000013380012 is similar to troxerutin (0.717), icariin (0.707), and monoxerutin (0.703).Troxerutin has been previously reported to exhibit a growth inhibitory effect on HeLa cells, with an effective IC 50 of 320 mg/mL after 24 h of treatment [140].This inhibitory effect was associated with increased levels of apoptosis-related proteins (Bax, cleaved caspase-3, and TNF-α) and a decrease in anti-apoptotic Bcl-2 protein expression [140].The Bcl-2 inhibitor ABT-737 has been shown to be beneficial in treating post-transplant EBV-related B lymphoproliferative disorders [110].Troxerutin also induced oxidative stress by elevating malondialdehyde concentration while reducing the activity levels of glutathione peroxidase and superoxide dismutase [140].The anticancer properties of icariin have been documented in the literature [141][142][143][144][145][146].The anticancer properties of icariin against human HeLa cervical cancer cell line and normal HCvEpC cells were assessed previously [146].Icariin demonstrated dose-dependent inhibition of HeLa cell growth, with an IC 50 of 20 µM, while exhibiting minimal toxic effects on normal HCvEpC cells [146].Like troxerutin [140], the anticancer effects of icariin were attributed to the induction of apoptosis, marked by caspase 3 and 9 cleavage, the upregulation of pro-apoptotic Bax, and the downregulation of anti-apoptotic Bcl-2 [146].Additionally, icariin triggered the generation of reactive oxygen species (ROS), reduced mitochondrial membrane potential, and inhibited the mTOR/PI3K/AKT signaling cascade in HeLa cells, indicating its potent anticancer activity [146].
Another study investigated the antitumor effects of icariin in U14 tumor-bearing mice and SiHa cells [147].Icariin exhibited a dose-dependent suppression of tumor tissue growth in vivo and significantly reduced the viability of SiHa cells in vitro.In SiHa cells, icariin was found to inhibit migration, invasion capacity, and the expression of inflammatory markers such as TGF-β1, TNF-α, IL-6, IL-17A, and IL-10 [147].Additionally, icariin promoted apoptosis in cervical cancer cells by downregulating key apoptotic markers, including Ki67, survivin, Bcl-2, c-Myc, and upregulating P16, P53, and Bax levels, both in vivo and in vitro [147].The study further revealed that icariin exerted its effects through the impairment of the TLR4/MyD88/NF-κB and Wnt/β-catenin pathways [147].In another study, the role of icariin in modulating malignant behaviors by targeting the miR-875-5p/MDM4 axis in cervical cancer was investigated [148].The study revealed downregulation of miR-875-5p in cervical cancer cells, promoting malignant phenotypes and autophagy and inhibiting apoptosis [148].Conversely, the overexpression of miR-875-5p demonstrated opposing effects.MiR-875-5p was identified as a sponge for MDM4, and elevated MDM4 expression attenuated the impact of miR-875-5p mimic on autophagy and epithelial-mesenchymal transition.However, icariin was found to counteract the inhibitory effects of miR-875-5p inhibitor on autophagy and xenograft tumor growth [148].These suggest that ZINC000013380012 may serve as a promising lead molecule for the development of cervical cancer therapy, offering a potential alternative to current treatment modalities.
Compound ZINC000014588133 (melicobisquinolinone B) is found in the leaves and branches of Melicope denhamii [149] and M. pteleifolia [150].Melicobisquinolinone B has previously been shown to possess a weak anti-proliferative activity against DLD-1 human colon cancer cells [151], which suggests the need for further investigations into its potential as an anticancer agent or its optimization for enhanced efficacy.ZINC000014588133 is structurally similar to vintafolide, vinorelbine, and anhydrovinblastine, with scores of 0.704, 0.701, and 0.7, respectively.Vintafolide, antigen drug conjugate [152], has been shown to be effective in ovarian cancer by targeting the folate receptor [153][154][155][156]. Vinorelbine is used for non-small cell lung cancer treatment [157][158][159][160] and has also proven to be active against cervical cancers with moderate toxicity [161,162].Anhydrovinblastine serves as the precursor for vinblastine [163,164], an anticancer drug with a broad spectrum of activity [165,166].Moreover, ester anhydrovinblastine derivatives have been reported to demonstrate potent cytotoxicity against human non-small cell lung cancer (A549) and cervical epithelial adenocarcinoma (HeLa) cell lines [167].

Conclusions
High-risk HPV types 16 and 18 are known to cause cervical cancer in women.EBV is also believed to contribute to cervical cancer progression by elevating HPV E6 level, which in turn, suppresses p53 production.Cervical cancers are the fourth most prevalent cancers in women globally and a leading cause of death.Finding drug candidates that can inhibit cervical cancer growth and progression will be a major breakthrough.This study therefore sought to use a robust computational drug discovery pipeline to shortlist natural compounds that could target E6 proteins of HPV16 and HPV18 simultaneously.The shortlisted compounds that exhibited favorable docking scores against both HPV E6 proteins were subjected to another round of scrutiny by reassessing their binding affinity against the EBNA1 protein.This step aimed to identify potential inhibitors capable of targeting both HPV E6 and EBNA1, with the ultimate goal of providing therapeutic benefits for cervical cancer patients experiencing HPV/EBV co-infection.The good candidates were further subjected to ADMET profiling and toxicity testing.Biological activity predictions supported the antiviral and anticancer activities of the most promising candidates, setting the stage for MD simulations and MM/PBSA calculations.The biological activity predictions revealed that the compounds possess anticancer properties and may function

Figure 1 .
Figure 1.Step-by-step procedure employed in this study to identify potential small molecule inhibitors of HPV E6 and EBNA1.Compounds from traditional Chinese medicine were screened against EBNA1 and E6 proteins of HPV16 and HPV18.Top compounds were subjected to ADMET profiling and biological activity predictions.The shortlisted compounds were further subjected to molecular dynamics simulations and MM/PBSA computations to evaluate their binding affinities.
2, x FOR PEER REVIEW 11 of 40

Figure 2 .
Figure 2.Bar plot illustrating the count of compounds categorized as "high", "low", or "none" for each of the four toxicity risks, including mutagenicity, tumorigenicity, reproductive effect, and irritant.

Figure 2
Figure 2.Bar plot illustrating the count of compounds categorized as "high", "low", or "none" for each of the four toxicity risks, including mutagenicity, tumorigenicity, reproductive effect, and irritant.

Figure 3 .
Figure 3. Protein-ligand interaction profile for EBNA1 in complex with (a) VK-1727 and (b) ZINC000085568136.Ligand is represented as purple, hydrophobic bonds are represented as red arc spokes, and hydrogen bonds are represented as green dashed lines.

Figure 3 .
Figure 3. Protein-ligand interaction profile for EBNA1 in complex with (a) VK-1727 and (b) ZINC000085568136.Ligand is represented as purple, hydrophobic bonds are represented as red arc spokes, and hydrogen bonds are represented as green dashed lines.

Figure 10 .
Figure 10.Chemical structures of the three known inhibitors (VK-1727, baicalein, and gossypetin) and the six potential anti-cervical cancer compounds identified in this study.

Figure 10 .
Figure 10.Chemical structures of the three known inhibitors (VK-1727, baicalein, and gossypetin) and the six potential anti-cervical cancer compounds identified in this study.

Table 1 .
Docking scores of the known inhibitors and shortlisted compounds after virtual screening

Table 4 .
Average RMSD, Rg, RMSF, and number of hydrogen bonds computed using data from the 100 ns MD simulations of the unbound protein and protein-ligand complexes.

Table 5 .
Energy contributions from MM/PBSA calculations for the protein-ligand complexes.All energy values are in kJ/mol.