Phytoconstituents of traditional Himalayan Herbs as potential inhibitors of Human Papillomavirus (HPV-18) for cervical cancer treatment: An In silico Approach

Human papillomavirus (HPV) induced cervical cancer is becoming a major cause of mortality in women. The present research aimed to identify the natural inhibitors of HPV-18 E1 protein (1R9W) from Himalayan herbs with lesser toxicity and higher potency. In this study, one hundred nineteen phytoconstituents of twenty important traditional medicinal plants of Northwest Himalayas were selected for molecular docking with the target protein 1R9W of HPV-18 E1 Molecular docking was performed by AutoDock vina software. ADME/T screening of the bioactive phytoconstituents was done by SwissADME, admetSAR, and Protox II. A couple of best protein-ligand complexes were selected for 100 ns MD simulation. Molecular docking results revealed that among all the selected phytoconstituents only thirty-five phytoconstituents showed the binding affinity similar or more than the standard anti-cancer drugs viz. imiquimod (-6.1 kJ/mol) and podofilox (-6.9 kJ/mol). Among all the selected thirty-five phytoconstituents, eriodictyol-7-glucuronide, stigmasterol, clicoemodin and thalirugidine showed the best interactions with a docking score of -9.1, -8.7, -8.4, and -8.4 kJ/mol. Based on the ADME screening, only two phytoconstituents namely stigmasterol and clicoemodin selected as the best inhibitor of HPV protein. MD simulation study also revealed that stigmasterol and clicoemodin were stable inside the binding pocket of 1R9W, Stigmasterol and clicoemodin can be used as a potential investigational drug to cure HPV infections.


Introduction
Malignant tumors have cancerous cells that grow indefinitely and proliferate locally or to the nearby body parts. They can also metastasize to distant parts via the bloodstream or the lymphatic system and affects the liver, lungs, brain, bone, and other areas. Some of the viral infections including Epstein-Barr, Human immunodeficiency virus (HIV), Human papillomavirus (HPV), Human herpesvirus 8 (HHV-8), hepatitis B and C, etc. have been well documented in increasing the risks of various cancers. Cervical cancer is identified as the fourth most common malignancy among women worldwide and even ranked after breast cancer with 2�1 million cases, colorectal cancer with 0�8 million cases, and lung cancer with 0�7 million cases worldwide [1]. Each year 0.5 million women are diagnosed with cervical cancer [2]. Viral infections are responsible for enhancing the prevalence and as well as the risks associated with approximately 15-20% of the cancers in the total human population, whereby numerous viruses play a key role in the etiology of some malignant cancers at multiple stages. HPV infection is considered to be one of the common risk factors for the development of cervical cancer, wherein about 99% of the cases are associated with cervical cancer [3,4]. HPV is a double-stranded DNA and a non-enveloped virus of the family Papovaviridae. It is 55 nm in diameter and has an icosahedral capsid comprising of 72 capsomeres, each one with at least two capsid proteins namely L1 and L2. HPV infection is widespread in humans and has been found in several animals as well, where it is unique to its host. It is reported that HPV is present in both symptomatic and asymptomatic populations and there are many more cancers that are cognate with HPV infections like anogenital, vulvar, vaginal, penile, head, and neck [5][6][7][8]. More than 120 HPV genotypes have been classified which infects human skin and among these only 13-15 types are associated with cancer involving high-risk and low-risk types depending on their oncogenic potential [9]. According to the earlier data, about 5% of the population in Europe is infected with HPV and 26% in the Sub-Saharan African countries [10], 21.4% in Eastern Europe, and 16.1% in Latin America [11]. In India, 6.3% of the population is Thalictrum foliolosum, Thymus serpyllum, Trillium grandiflorum, Viola odorata, and Zanthoxylum armatum which are traditionally used for the treatment of diabetes, gastrointestinal, hepatic, nervous, skin, and respiratory tract infections, and rheumatism [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33], were investigated for their biochemical compositions. Earlier also, many natural and plant-based compounds have been identified as promising sources for developing cytotoxic agents for the treatment and prevention of cancer in recent years [34]. Stigmasterol, an active ingredient of B. aristata has also been shown to exert anti-angiogenic and anti-cancer effects. downregulation of TNF-alpha and VEGFR-2 [35]. Therefore, the present research is focused on to study the intraction of major phytocompounds of twenty important medicinal plants of North west Himalays with target protein of HPV-18 capsid protein L1 and to study the toxicity of active phytocomplounds docking.

Ligand preparation
One hundred nineteen major phytoconstituents of twenty medicinal plants (Allium sativum, Artemisia annua, Berberis aristata, Bergenia ligulata, Cannabis sativa, Cymbopogon citrates, Dactylorhiza hatagirea, Moringa oleifera, Myristica fragrans, Oxalis corniculata, Picrorhiza kurroa, Piper nigrum, Pleurospermum brunonis, Podophyllum hexandrum, Rheum emodi, Thalictrum foliolosum, Thymus serpyllum, Trillium grandiflorum, Viola odorata and Zanthoxylum armatum) collected from the higher altitudes of Himachal Pradesh and Uttarakhand (India) were selected for the molecular docking studies. The 3-dimensional structures of the selected phytoconstituents and the two standard drugs of the HPV (podofilok and imiquimod) were downloaded from the PubChem database (www.pubchem.com) in.sdf formats. The.sdf files of the phytoconstituents were converted into PDB formats. All the selected ligands (phytoconstituents and drugs) were prepared using the open babel program from the command line on an Ubuntu terminal. The energy of all the selected phytocompounds was minimized by using chemdraw 3D 15.0. Phytoconstituents selected for the study, plant source, and their ethnomedicinal uses are enlisted in Table 1.

Protein preparation
Human papillomavirus type (HPV-18 E1) protein (PDB ID: 1R9W) [47] was selected for molecular docking with major bioactive phytoconstituents from the selected Himalayan medicinal plants as listed in Table 1 to identify its potential inhibitor. Podofilok and Imiquimod is used as drug target because of its small surface area and functional significance, the E1 dimerization interface is an obvious target for a small molecule disrupting DNA replication [47]. The 3-dimensional monomeric structure of HPV-18 E1 protein (1R9W) was downloaded from the protein databank (http://www.rscb.org/pdb). It consists of a pentamer structure, and its chain A was extracted for docking using Pymol software. Protein was prepared by using autoDock 1.5.6 program, during protein preparation all the water molecules were deleted, polar hydrogens were and kollaman charge (5.0) is added. The active site was predicted manually by grid box analysis (grid dimensions x = 40, y = 40, z = 40 Å, centered at x, y, z = 21.25, -0.773, 12.365 Å), TYR121, ASP124, VAL 126, SER129, TYR132, GLY133, GLY134, ASN135, PRO136, ASN259, ALA284, SER285, SER286, TYR288 are the key amino acids present in the active site of the protein.

Molecular docking using herbal phytoconstituents
The molecular docking of the selected ligands with the catalytic triad of the 1R9W protein was performed using AutoDock 1.5.6 program and responses were saved as a PDBQT file.

MD simulation of protein-ligand complexes
To gain deeper insights into the binding interactions of 1R9W protein with clicoemodin and stigmasterol, MD simulations were performed that accelerated the graphics processing unit (GPU) in Amber18 [54]. Amber generalized force field (GAFF) was assigned to the ligands, whilst amber ff14SB to the protein [55,56]. The atomic charges of the ligands were calculated using the restrained electrostatic potential (RESP) protocol at the HF/6-31G � level of theory 31, 32 using Gaussian 09 software [57]. The protonation states of the ionizable residues of protein structures were analyzed by pKa calculation at neutral pH using the H++ server. All the systems were prepared using tLeap program. Each system was neutralized with sodium ions and solvated using a cubic water box with TIP3P model. To eliminate the stearic clashes, each system was subjected to a total of four minimizations. First, all the sodium ions and water molecules were minimized by 2000 steps of steepest descent, followed by 3000 steps of the conjugate gradient algorithm. Consecutively, all the hydrogen and water molecules were relaxed using the same protocol. Finally, the whole system was energy-minimized for 5000 steps of steepest descent plus 5000 steps of conjugate gradients. The heat of the system was started from 0 to 300 K while running 200 ps of MD and, next, 300 ps of the density equilibration with position restraints on the protein atoms at a constant volume. Before performing the production step, all the protein-ligand systems were equilibrated with10 ns of MD without any positional restraints at a constant pressure. The temperature was maintained at 300 K by coupling to a Langevin thermostat using a collision frequency of 2 cm -1 . A cutoff of 10 Å was employed for unbonded interactions, the Particle Mesh Ewald (PME) [56] method and the SHAKE [57] algorithm were used to restrict the bond lengths involving hydrogen atoms. Finally, the MD simulations (productions) were performed using 100 ns at a temperature of 300 K. The generated trajectories were used to analyze the behavior of each complex and to access the stability of the system. The deviations of the protein and protein-ligand complex system were analyzed by calculating the root mean square deviation (RMSD), root mean square fluctuation (RMSF), a radius of gyration (RG), and solvent accessible surface area (SASA) parameters [58,59]. Furthermore, the total free energy of binding, the free energy of solvation (polar þ nonpolar solvation energies), and the potential energy (electrostatic and van der Waal's interactions) of each protein-ligand complex were calculated by molecular mechanics using Poisson-Boltzmann Surface Area (MM-PBSA) method [60]. The last 10 ns of the MD trajectory were taken for the calculation of MM-PBSA.

Molecular docking of some medicinally important phytoconstituents with target protein of HPV-18 (1R9W)
One hundred nineteen major phytoconstituents of twenty important medicinal plants of Northwest Himalayas are selected based on their reported medicinal attributes as evidenced in the earlier literature. These are widely recommended in the Ayurvedic/Unani medicines for the treatment of diabetes, gastrointestinal, hepatic, nervous, skin, and respiratory tract infections, and rheumatism (as summarized in Table 1). 3D structures of these phytoconstituents and their target proteins were retrieved from PubChem and PDB databases. Molecular docking results showed that out of one hundred nineteen phytoconstituents, only thirty-five phytoconstituents of twelve medicinal plants showed docking scores more than the standard drugs imiquimod (-6.1 kJ/mol) and podofilok (-6.9 kJ/mol) interactions with HPV-18 (1R9W). Important phytoconstituents with their respective docking score, plant sources, and interactive amino acids are summarized in Table 2 and S1 Table. From this study, we found that six phytoconstituents of Thalictrum foliolosum and Thymus serpyllum, five of Berberis aristata, four of Rheum emodi, three each from Cannabis sativa and Picrorhiza kurroa, two each from Moringa oleifera and Zanthoxylum armatum, and one each from Myristica fragrans, Oxalis corniculata, Piper nigrum, and Pleurospermum brunonis showed docking score more than the standard anti-cancerous drugs.

In silico pharmacokinetics and toxicity prediction of the bioactive phytoconstituents
The pharmacokinetics and toxicity of pharmacologically important phytoconstituents were analyzed in silico using computational facilities of SWISS ADME, admetSAR, and Protox II servers. Thirty-five active phytoconstituents of twelve medicinal plants and two standard drugs (Podofilox and Imiquimod) were selected for the pharmacokinetics and toxicity studies. Among all the thirty-five phytoconstituents, only 29 phytoconstituents follow Lipinski's rule of five, and phytoconstituents such as 3-O-b-galactopyranoside, eriodictyol-7-glucuronide, glucomoringin, pikuroside, rutin, and verbascoside do not follow Lipinski's rule of five, therefore these cannot be used as investigational drug candidates. Whereas 30 phytocompounds follows  Table 3 and S1-S3 Figs.
Based on the molecular docking drug-likeness data and toxicity data, two best phytoconstituents (stigmasterol from Berberis aristata and clicoemodin from Rheum emodi) were selected for further MD simulation analysis.

MD simulation of protein-ligand complexes
The structural stability of the protein-ligand complex was determined using RMSD analysis during a specific time (100 ns). Molecular dynamic simulation of the selected protein-ligand complexes revealed that stigmasterol is stable inside the binding pocket of 1R9W. RMSD of stigmasterol was stable from the beginning of the simulation (0 ns) and remain stable up to 100ns period, which indicates fitness of the ligand inside the binding pocket of the target protein. RMSD has a little fluctuation between 0.5-1 Å, which is acceptable (Fig 3A). Similarly, the RMSD of clicoemodin, in complex with 1R9W is stable from the beginning (0 ns) and remains stable up to 100 ns periods, and only a little fluctuation was observed which is in the acceptable range (0.5-1 Å) (Fig 3B). The RMSF value reflects the particle's deviation from the initial position of macromolecules in a protein's three-dimensional structure. The RMSF plots of stigmasterol and clicoemodin in complexation with 1R9W revealed that the interactive amino acids of stigmasterol and clicoemodin have lesser fluctuations in the alpha-helix and beta-strand of the target protein (Fig 4).
The solvent-accessible range of stigmasterol is 8500-9500 Å, whereas in the case of clicoemodin, it lies between 8000-8500 Å (Fig 5). The radius of gyration plot established the compactness of stigmasterol and clicoemodin complexes with 1R9W, which confirms the stability of the protein-ligand complexes (Fig 6).

PLOS ONE
To understand the impact of inhibitors upon complexation in terms of their binding affinities, MM-PBSA binding free energy method was utilized to calculate the binding free energies and the energy components of their complexes are shown inTable 4.
As it is evident in Table 4, the total binding free energies (ΔGbind) of clicoemodin-1R9W and stigmasterol-1R9W were -25.63 kcal/mol and -18.13 kcal/mol, respectively. Accordingly, among the two studied complexes, Clicoemodin-1R9W depicted the most favorable of ΔGbind with its lowest values of -25.63 kcal/mol. At this point, it is interesting to address the key contributions of each binding component that can impose to the total binding free energies.

Discussion
HPV has been described as the principal cause for the development of cervical cancer, and this potential is linked with its ability to produce various oncoproteins. Women infected by the human immunodeficiency virus (HIV) also have a very high susceptibility to HPV infections with the presence of multiple HPV types being typically reported [61]. HPV infections cause several benign proliferations including warts, epithelial cysts, anogenital, oro-laryngeal and -pharyngeal papillomas, intraepithelial neoplasias, keratoacanthomas, and other forms of hyperkeratoses. Brown and co-workers detected a high-risk HPV viral DNA on condylomata acuminata lesions from HIV immunosuppressed and otherwise healthy patients, always in association with a low-risk HPV type, HPV 6 or 11 [62].   In the last decades, the pharmaceutical industry has been employing computer-aided drug design (CADD) techniques to accelerate drug development, intending to reduce time, costs, and failures, mainly in the final stage [63,64]. CADD attempts to improve the efficacy of the hit discovery process through in silico methodologies by testing and screening a large number of compounds to identify a small group of candidates with desirable pharmacological properties. The virus particle binds to its specific cell surface receptors, which allows it to enter the host cell. Virions can first bind to the basement membrane before being transferred to the basal keratinocyte cell surface [65,66]. The structure-based approaches, including docking and the application of molecular dynamics simulations, use the 3D structure of the target molecule to screen potential ligands. These methods evaluate ligand recognition by the target molecule and the prediction and characterization of binding sites as well as their binding affinities [63,67].  Computational-based approaches hold a great promise for improving drug development and revolutionizing clinical research by providing appropriate interventions for women diagnosed with cervical cancer. In the current study, we found that among all the selected one

PLOS ONE
Traditional herbs as potential inhibitor of Human Papillomavirus Cervical Cancer hundred nineteen phytoconstituents of the twenty unexplored medicinal plants of northwestern Himalaya, stigmasterol from Berberis aristata, and clicoemodin from Rheum emodi can be used as a potent drug to cure the HPV. In our previous research, Kumar et al. [68] studied the interactions of traditionally used phytoconstituents (ferulic acid, ursolic acid, jaceosidin, curcumin, gingerol, indol-3-carbinol, silymarin, resveratrol, berberine, and artemisinin) against E6 oncoprotein of high-risk HPV 18 by using Autodock 4.2 tool. Sookmai et al. [69] reported the anti-papillomavirus activity of Clinacanthus nutans extract; they also suggested C. nutans as prevention of HPV. Earlier a study of Moghadamtousi and co-workers, also showed a high potential of curcumin to cure high-risk HPV 16 and 18 [70]. Similarly, various researchers revealed the anticancerous potential of stigmasterol in the treatment of skin carcinoma [71], ovarian cancer [72], and gastric cancer [73,74]. Previous reports also revealed that crude extracts or phytoconstituents obtained from medicinal plants of the Northwest Himalayas have various therapeutic attributes and show synergism of the in silico and in vitro investigations [75][76][77][78][79][80][81][82][83] and even Rolta et al. [84] reported the antiviral activity of phytoconstituents of important medicinal plants of the Northwest Himalayas against SARS-CoV-2. According to Cha et al. [85], emodin a major phytocompound of Rheum emodi inhibits the cell growth of the bladder cancer cell line. Thymus serpyllum essential oil and its two phytoconstituents thymol and carvacrol have antitumor activity against P815 mastocytoma cell line [86], HeLa (cervical carcinoma), and HepG2 (hepatocellular carcinoma) cells as reported by Hussain and coworkers, 2013 [87].
In silico docking examination by Muthukala et al. [88], has reported the inhibitory potential of quercetin and quercetin-like components towards human cervical cancer cell line proteins. Scheffner and coworkers also investigated that HPV oncoproteins i.e. E7 and E6 are activated in cervical cells during the progression of cervical cancer and also have been recognized to interact with tumor suppressor proteins viz. p53 and pRB respectively [89]. Through an in-silico approach, Kotadiya and Georrge [90] detected some putative drugs derived from the natural products against HPV, while Samant et al. [91] investigated the molecular interactions of HIV and anti-viral drugs targeting HPV-18 E6 protein. Inhibitory potential of natural products including colchicine, daphnoretin, curcumin, ellipticine, epigallocatechin-3-gallate, etc. towards HPV-16 E6 protein was also recognized by molecular docking [92]. Phytocompounds of Ficus carica is also reported as anticacinogenic drug [93]. Drug delivery systems can be developed to circumvent possible high toxicity and low availability of the compounds and, for instance, to combine this approach of L1 inhibitors with gene therapy to induce cancer cell apoptosis. Considering the costs of screening and treatment of cervical cancer and knowing that the highest incidence occurs in developing countries, a more cost-effective treatment is needed.

Conclusion
Major phytoconstituents of twenty traditionally important plants were screened against the HPV by molecular docking. Among the selected phytoconstituents stigmasterol from Berberis

PLOS ONE
aristata, and clicoemodin from Rheum emodi were found to be the best inhibitors of HPV-18 (1R9W protein). ADME/T validated the proposal of stigmasterol and clicoemodin are safer drugs for the treatment of cervical cancer. MD simulations data showed that stigmasterol and clicoemodin are stable inside the binding pocket of the target protein.
Since this study has only contemplated computational analysis, it is still necessary to validate the bioactivities of these important phytoconstituents in vitro and in vivo models before proposing their antiviral potential. We, also suggest that further research should take these results into account in both the discovery as well as the lead optimization stages to achieve the successful development of anti-HPV drugs.
Supporting information S1 Table. Binding energy of selected phytocompounds with standard drugs.