In Silico Analysis of the Apoptotic and HPV Inhibitory Roles of Some Selected Phytochemicals Detected from the Rhizomes of Greater Cardamom.

Occurrence of cervical cancer, caused due to persistent human papilloma virus (HPV) infection, is common in women of developing countries. As the conventional treatments are expensive and associated with severe side effects, there is a need to find safer alternatives, which is affordable and less toxic to the healthy human cells. Present study aimed to evaluate the anti-HPV and apoptotic potential of four compounds from the greater cardamom (Amomum subulatum Roxb. var. Golsey), namely rhein, phytosphingosine, n-hexadecenoic acid and coronarin E. Their anti-HPV and apoptotic potential were studied against viral E6, E7 and few anti-apoptotic proteins of host cell (BCL2, XIAP, LIVIN) by in silico docking technique. Phytochemicals from the plant extract were analysed and identified by LC/MS and GC/MS. Involvement of the target proteins in various biological pathways was determined through KEGG. Structural optimization of the three-dimensional structures of the ligands (four phytochemicals and control drug) was done by Avogadro1.1. Receptor protein models were built using ProMod3 and other advanced tools. Pharmacophore modelling of the selected phytochemicals was performed in ZINCPharmer. Swiss ADME studies were undertaken to determine drug likeness. The ligands and proteins were digitally docked in DockThor docking program. Protein flexibility-molecular dynamic simulation helped to study protein-ligand stability in real time. Finally, the correlation of evaluated molecules was studied by the use of principal component analysis (PCA) based on the docking scores. All the ligands were found to possess apoptotic and anti-cancer activities and did not violate Lipinsky criteria. n-Hexadecanoic acid and its analogues showed maximum efficacy against the target proteins. All the protein-ligand interactions were found to be stable. The uncommon phytochemicals identified from rhizomes of greater cardamom have anti-cancer, apoptotic and HPV inhibitory potentials as analysed by docking and other in silico studies, which can be utilized in drug development after proper experimental validation.

Introduction

Cervical cancer, the fourth most common cancer in women, is a cause of concern in countries with low human development index (HDI). According to Globocan 2020, there were 604,127 new cases and 341,831 deaths across the globe, while 85% of these deaths were from the developing world [1,2] and more than one third cases were from India and China. Several factors, mainly, weakened immune system, poor nutritional status, multiparity, having multiple sexual partners, lack of proper hygiene practices, smoking, uncontrolled use of oral contraceptives, poor diagnosis and genetic predisposition are the common risk factors for developing cervical cancer. Human papilloma virus (HPV) is the primary causal organism for this cancer [3]. Conventional treatment involves chemotherapy, surgery and radiotherapy, which is not only associated with severe side effects but also unaffordable for the economically weaker section. Considering this backdrop, it is of paramount importance for the rapid development of an inexpensive and safe therapeutic regime to treat this virus-induced cancer.

HPV is a non-enveloped DNA virus belonging to the family Papillomaviridae with a diameter of 50–55 nm. There are as many as 300 genotypes of this virus; however, 200 types are pathogenic to humans. Among all these genotypes, HPV16 is considered the high-risk category. With 7.9 Kb genomic size, HPV16 codes for many viral oncoproteins, including E6 and E7. E6 and E7 are majorly responsible for integration and viral replication within the host cell, inducing several cancer markers such as angiogenesis, proliferation, metastasis, hyper telomerase activity, etc. [3], and at the same time, the anti-apoptotic proteins such as BCL2, XIAP, LIVIN [4,5], etc. upregulate in the HPV-infected cancer cells can be a targeted therapy to induce cell survival.

To date, therapeutic options are scarce for the treatment of HPV-induced cervical cancer. Vaccines are the most promising strategy for the prevention of HPV persistence. However, high expense and the presence of multiple variants of the virus currently restricted the universal applications of such vaccines, especially in the developing nations [6,7]. Recently, Gao et al. (2020) [8] introduced poly(amide-amine)-poly(β-amino ester) hyperbranched copolymer/linear poly(β-amino ester)-CRISPR/Cas9 therapeutic plasmids polyplex to treat papilloma virus-infected cervical cancer. Zamulaeva et al. (2021) [9] concluded that radiation might not be very successful for such carcinogenesis due to multiple factors. In a review article, Moga et al. (2021) [10] showed that several bioactive molecules from the seaweeds, such as glucans, lutein, carotene, stypoldione, etc., could be the potential anti-HPV anti-cancer agents.

Recent advancements in computer-aided techniques gave acceleration in the momentum of the drug discovery screening strategy. Inclusion of in silico docking experiments for screening of candidate drugs before wet-lab validation can expedite the process and reduce the cost significantly. Bioinformatics witnesses many applications in the drug discovery process [11,12]. Natural products are being considered safe and effective therapeutics when compared with their synthetic counterparts. However, the purification of natural compounds, especially those are sourced from botanical origin, is expensive and time-consuming. Therefore, screening a few potential bioactive compounds among many molecules can give a significant advantage in the drug discovery process and integration of targeted molecular docking technique may provide a substantial edge over it. A few activities of natural products against anti-HPV proteins were reported in the literature. Prasasty et al. (2017) [13] studied five natural products against HPV proteins by using the molecular docking technique. Recently, Kotadiya et al. (2020) [14] screened 17,944 herbal molecules to identify anti HPV E6 and E7 agents by applying a similar strategy.

Certain phytochemical classes such as labdane diterpene, anthraquinone, long-chain fatty acid and amino alcohols are reported to have significant activities. Hsieh et al. (2020) [15] showed that labdane diterpenes could induce apoptosis in 5-fluorouracil resistance human oral cancer cells. Furthermore, the role of anthraquinone moiety against a range of carcinomas, such as breast carcinoma, colon carcinoma, prostate carcinoma and neuroblastoma was also reported in the literature [16]. The role of long-chain fatty acid and amino alcohol as anti-cancer therapeutic options [17] were also reported. Though, significant literature is available on the anti-cancer properties of these phytochemical classes, their therapeutic potential against HPV-induced cervical cancer is yet to be established. The botanical family Zingiberaceae is considered the powerhouse of phytochemicals and reported to contain such phytochemical classes in abundance [18]. Amomum subulatum belongs to Zingiberaceae and has substantial pharmacological properties such as antioxidant, anti-inflammatory, anti-viral and anti-cancer activities [19,20].

However, the scientific data on the effect of major phytochemicals of large cardamom on HPV-induced cervical cancer is limited to the best of our knowledge.

In this context, our study aimed to evaluate anti-HPV potential of four identified uncommon compounds from the greater cardamom (Amomum subulatum), namely rhein, phytosphingosine, n-hexadecenoic acid and coronarin E. Furthermore, anti-HPV potential of these compounds was studied against the target oncoproteins HPVE6, HPVE7 and anti-apoptotic proteins BCL2, XIAP and LIVIN by using in silico docking technique. Pharmacophore screening was performed to find a few novel candidate chemical entities against the target proteins, using the phytochemicals as templates. Finally, the correlation of evaluated molecules was studied by using the statistical tool principal component analysis (PCA) based on the docking scores.

Material and Methods

Materials

Plant Material

Rhizomes of the plant Amomum subulatum var. Golsey was obtained from Kabi research farms of Indian Cardamom Research Institute in Pangthang, Sikkim, and the identified plant specimen was submitted to the Calcutta University Herbarium (CUH) Kolkata, India (accession number 20087).

Chemicals

Ethanol (LC–MS grade, SRL, India), hexane (GC–MS grade, SRL, India) and deionized double distilled water (Milli Q).

Methods

Preparation of Plant Extract

The rhizomes were cut from the main stem, washed in tap water to clean off the mud and soil matter and cut into small pieces and allowed to air dry until complete moisture loss. The dried pieces were ground to powder in an electronic grinder. Solvent extraction was carried out in 95% ethanol with the ground sample in a 1:10 ratio, including three intermittent solvent changes by filtration through non-absorbent cotton for over 10 days. The Amomum subulatum ethanolic extracts were then set to desired concentrations by air drying and re-dissolving the remnant extracted matter in the solvent.

LC/MS (Liquid Chromatography/Mass Spectrometry) Analysis

LC/MS analysis was carried out at the Sophisticated Analytical Instrument Facility (SAIF), Indian Institute of Technology, Mumbai. Ethanolic extract of A. subulatum (1 mg/mL) was subjected to liquid chromatography/mass spectrometry (LC/MS) analysis (Agilent Technologies) coupled with Q-TOF (quadrupole time of flight) mass spectrometer (model — G6550A) with a dual AJS ESI ion source, a HiP sampler with binary pump (model — G422013) with a 30-min run time at a flow rate of 0.3 mL/min and pressure limit 1200 bar. The injection volume was set as 5 μL and the chromatography column was Hypersil gold (3 μm particle size, 100 × 2.1 mm dimension). The flow gradient was maintained as 95% water and 5% acetonitrile for 0–18 min, 100% acetonitrile for 18–26 min and 95% water and 5% acetonitrile for 26–30 min. Phytochemicals were identified by Agilent LC/MS Phytochemical library.

GC/MS (Gas Chromatography/Mass Spectrometry) Analysis

A total of 2 g of dried extract was dissolved in 10 mL of hexane at 37 °C for 3 days. Residual water was removed by adding a pinch of sodium sulphate (Merck), followed by low-speed centrifugation and then finally filtered by 0.45-micron nylon membrane for gas chromatography/mass spectrometry (GC–MS) analysis. The analysis was performed in a GC interface with a mass spectrometer using an SH-Rxi™-5Sil (Shimadzu) column. The oven temperatures ranged from 70 to 260 °C. A total of 2 μL of the filtered extract was injected. A 49-min run was programmed at 5 °C/min ramping with initial hold time for 1 min and final hold time for 10 min. Helium was used as the carrier gas at a flow rate of 1 mL/min. The injection port was set at 250 °C. MS operating parameters include electron impact ionization at 70 eV with a mass range of 50–600 amu. NIST mass spectral library 2017 (Version 1.0) (Version 1.0) [21].

Pathway Analysis

KEGG or Kyoto Encyclopaedia of Genes and Genomes (https://www.genome.jp/kegg/) is an open-source database containing high-level functional information of the biological systems such as the cell, the organism and the ecosystem [22,23]. Diverse pathways based on a large molecular and high throughput experimental data are available in this database. For this study, we identified individual pathways involving the target proteins (HPVE6, HPVE7, XIAP, LIVIN and BCL2) of HPV-induced carcinogenesis and based on these pathways, an integrated pathway was constructed.

Preparation of Ligands

Based on the GC/MS and LC/MS results, four major phytochemicals, coronarin E (CID 9,971,144), n-hexadecanoic acid (CID 985), phytosphingosine (CID 122,121) and rhein (CID 10,168), and three dimensional (3D) structures of the selected phytochemicals and control kaempferol (CID 5,280,863) were downloaded from the PubChem database. PubChem is an open-source database (https://pubchem.ncbi.nlm.nih.gov/) of the National Institute of Health (NIH) and contains data on chemical structures, identifiers, chemical and physical properties, biological activities, patents, health, safety, toxicity data, etc. [24]. Furthermore, structural optimization was performed by Avogadro 1.1 software by applying a universal force field (UFF) algorithm. It utilizes frequencies, single-point energies and geometric optimization data [25].

Target Proteins (Receptors)

Three dimensional (3D) structures of the proteins, HPV (human papillomavirus) oncoprotein E6 (PDB id 4GIZ, 2.55 Å, X ray diffraction), BCL-2 (PDB id 2XA0, 2.70 Å, X ray diffraction) and BIR-3 domain of XIAP proteins (PDB id 1F9X, solution NMR), were downloaded from the RCSB Protein Data Bank (https://www.rcsb.org/).

Due to the unavailability of the 3D structures of the proteins, human papillomavirus oncoprotein E7 and LIVIN in the RCSB PDB, we modelled these proteins by using SWISS-MODEL server (https://swissmodel.expasy.org/) in two steps. In the first step, protein sequences of E7 (accession no. AAD33253.1) and LIVIN (accession no. AAQ89195.1) were downloaded from NCBI protein resources (https://www.ncbi.nlm.nih.gov/protein). In the second step, the sequences were uploaded to the SWISS-MODEL server for further modelling. SWISS MODEL relies on structural template(s) identification, target sequence and template structure(s) alignment, model-building and quality evaluation to build the model from the template library. It uses modelling engine ProMod3 and advanced tools such as Monte Carlo, graph-based TreePack algorithm, CHARMM22/CMAP force field and SCWRL4. The quality evaluation of the modelled structures was done by Ramachandran plot and QMEAN [26]. Structural rendering was performed by the open-source software UCSF Chimera 1.17 software [27].

Pharmacophore Screening Phytochemical Analogues

The pharmacophore is defined as a three-dimensional molecular interaction of ligand and receptor based on the chemical features such as hydrogen bond (H-bond), ionic charges, lipophilic and aromatic contact and hydrophobic groups. These chemical features or descriptors of the ligand are grossly responsible for a biological and functional response towards the receptor or the target proteins. Here, identified phytochemicals were used as pharmacophore candidates and chemicals features were deduced by uploading the ZINCPharmar server (http://zincpharmer.csb.pitt.edu/pharmer.html) [28], by using ‘add feature’ function. ZINCPharmar is an open access server used to screen small molecule entities against submitted candidate pharmacophore models. ZINC database contains 22,724,825 purchasable molecules. Among the output results, the top three small molecules were selected for each of the phytochemical inputs based on the RMSD values and three-dimensional structures of the ZINC ligands were downloaded from the database.

Evaluation of Drug Likeness

The canonical smile formats of the selected four phytochemicals and top three ZINC outputs were downloaded from the PubChem and uploaded in the open source SwissADME server for evaluation of drug-like properties of the phytocompounds. The drug-like properties of molecules were evaluated by multiple parameters such as lipophilicity (iLogP) water solubility (ESol), Lipinski rule of 5 and topological polar surface area (TPSA). Furthermore, lipophilicity (WLOGP) and TPSA were plotted in the Boiled-Egg graphical interface. This graph further explains the bioavailability of the drug through PGP (p-glycoprotein) responses. PGP is a membrane-bound transporter that induces efflux of the substrates (PGP +), thereby reducing the intracellular concentrations and inhibiting the drug’s bioavailability [29].

Molecular Docking

After preparation of ligands and proteins, molecular docking was performed by the web server DockThor (https://dockthor.lncc.br/v2/) [[[[. DockThor uses MMFFLigand and PdbThorBox inhouse tools for it docking algorithm along with the MMFF94S53 force field [[[[. For docking set up, grid parameters are shown in Table 1. Two anti-cancer drugs, namely kaempferol and ciglitazone were used as natural and synthetic controls, respectively.

Table 1

Grid parameters

S/N	Proteins	Centre grid parameters
1	XIAP	− 41.9855 × 139.931 × 15.3355
2	LIVIN	112.459 × 137.2445 × 174.7715
3	BCL2	33.446 ×  − 14.0395 ×  − 14.7705
4	E6	2 × 60 × 29
5	E7	16.954 × 22.4675 × 40.3375

Principal Component Analysis (PCA) and Molecular Pharmacophore Alignment

PCA is a multivariate statistical that is used to understand the dimensionality of a large set of data. It transforms the measured variables into unrelated ones, i.e. principal components (PCs). Each of the PCs is orthogonal, i.e. unrelated, and PC1 covers maximum variability of the dataset followed by subsequent PCs (PC2, PC3, etc.). PCA aids in clustering the data points based on their inherent variations [33,34]. In this study, the docking experiment’s binding energies (kcal mol−1) were used as inputs of PCA and thereby, it clusters the ligand molecules based on their similar properties. Furthermore, based on the best clusters identified, molecular pharmacophore alignment was performed by an open-source PharmaGist web server (https://bioinfo3d.cs.tau.ac.il/PharmaGist/php.php) [35] to determine the common chemical features of compounds. PharmaGist relies on the structure of the ligand rather than the receptor protein structure.

Protein Flexibility-Molecular Dynamic (MD) Simulation

The protein flexibility of the top-ranked ligand–protein complexes was evaluated by CABS-flex 2.0 server (http://biocomp.chem.uw.edu.pl/CABSflex2) and presented with RMSF (root mean square fluctuation). CABS-flex offers fast protein flexibility simulation and generates protein dynamic simulation at highly reduced system requirements. Flexibility simulation obtained from this server was reported to highly correlate with the NMR results [36,37]. CABS-flex provides high resolution (10-ns) protein near-native protein dynamics simulation and hence is very effective for evaluation of protein–ligand stability on real-time basis. Simulation in CABS-flex was set with default parameters, with 50 cycles.

Results and Discussion

Identification of Phytochemicals

In this study, four compounds, namely coronarin E (labdane diterpenes), n-hexadecanoic acid/palmitic acid (long chain fatty acid), phytosphingosine (amino alcohol) and rhein (anthroquinone) were identified from the rhizome extract of Amomum subulatum. The chromatograms are represented in the supplementary figures (Figs. S1 and S2).

Sirat et al. (1994) [38] isolated two labdane diterpenes, namely labda-8(17),12-diene-15,16-dial and coronarin E, from the chloroform extract of Alpinia javanica rhizomes. Furthermore, Subramanyan et al. (2020) [39] reported coronarin compound from Hedychium flavescent rhizome extract. A similar finding was reported from the ethyl acetate fraction of Amomum maximum [40]. We have identified coronarin E (0.4%) from the ethanolic extract of Amomum subulatum rhizome by GC–MS, which is in agreement with these findings. In addition to that, we found n-hexadecanoic acid (palmitic acid) in abundance (9.49%) in the extract. In literature, the presence of palmitic acid was commonly reported from various plants. Ali et al. (2019) [41] found high content of palmitic acid from Piper nigrum, Nigella sativa, Cinnamomum zeylanicum and Elettaria cardamomum extracts.

LC–MS analysis in our study revealed the presence of two major compounds rhein and phytosphingosine, in the ethanolic extract of A. subulatum. The phytochemical rhein was reported in a diverse number of plants. Cunha et al. (2017) [42] reported rhein as one of the major compounds in the ethanolic extract of Cassia bakeriana. Furthermore, Thanh-Tam Ho et al. (2019) [43] identified rhein from the plant Polygonum multiflorum. Nevertheless, to the best of our knowledge, we reported this compound (2.06%) for the first time from A. subulatum. Finally, we reported the presence of an uncommon phospholipid, phytosphingosine (3.2%) for the first time from A. subulatum rhizome. However, phytosphingosine was reported earlier from other plants. Agil et al. (2020) [44] identified 2.5% of phytosphingosine from the ethyl acetate fraction of Marsilea crenata Presl. leaves. This compound was also reported from the homoeopathic preparation of Gymnema sylvestre [45].

Pathway Analysis

Human papilloma-virus (HPV) proteins E6 and E7 play a crucial role in the molecular pathogenesis of cervical cancer. Literature indicated that these proteins interact with various cell cycle regulatory proteins to establish cervical cancer prognosis [46-48]. BCL2 is considered one of the most critical anti-apoptotic proteins related to lymphoma development [49]. This protein is also reported to relate with the manifestation of cervical cancer directly. Wang et al. (2019) [50] showed that upregulation of anti-apoptotic protein BCL2 expression through Homeobox C6 (HOXC6) transcription factor led to cervical cancer progression. Two other proteins, namely XIAP (X-linked inhibitor of apoptosis) and LIVIN, are also known to have a significant role in cervical cancer progression. XIAP is one of the members of the inhibitors of apoptosis proteins (IAP) family, and it acts by blocking the maturation and activation of caspase 3, 7 and 9 [51,52]. Overexpression of XIAP was reported with direct or indirect inhibition of caspase activity, thereby interfering with the related pathways of cell death, cell cycle, autophagy and cell migration [53]. The vascular endothelial cell angiogenesis helps to develop new blood vessels and supply oxygen and nutrients to the tumour cells. This angiogenesis was directly related to the overexpression of HPV E6 and E7 proteins [54,55]. Furthermore, it was observed that extracellular vesicles might take an essential role in cancer progression by transferring their content from a donor to cancerous cells. In cervical cancer, overexpression of anti-apoptotic proteins like XIAP and LIVIN was reported from these extracellular vesicles [56]. LIVIN, a member of IAP family, was first reported in 2001, and authors showed that silencing of this protein could trigger apoptosis specifically in the cancer cell lines [57]. We used KEGG pathways to illustrate the role of the anti-apoptotic and HPV proteins in the development of cervical cancer. Three KEGG pathways, namely cancer (map05200), apoptosis (map04210) and HPV infection (map05165), were selected (Figs. S3–S5), and finally, a combined pathway was constructed, see Fig. 1. The created diagram indicated that targeting these proteins may inhibit the overall HPV-induced cervical cancer pathway.

Fig. 1

KEGG pathways derived integrated HPV-induced canonical cervical cancer pathways and corresponding therapeutic consequences resulting from drug application leading to apoptotic cellular death or cell cycle arrest

Pharmacophore Screening for Phytochemical Analogues

Pharmacophore-based models were screened based on the descriptors of the identified phytochemicals. For phytochemicals, applications of such techniques are common in the literature [58,59]. Based on pharmacophore-based screening, Sharma et al. (2020) [60] filtered ten compounds and identified the phytochemical schinilenol as the best lead compound against the protein cyclin-dependent kinase 2 (CDK2). Friday et al. (2020) [61] designed ten analogues based on quercetin pharmacophore, and 7-(2,3-dihydroxycyclopropyl)-2-(3,4-dihydroxyphenyl)-3,5-dihydroxy-4H-1-benzopyran-4-one was selected as the best compound. Similar to these studies, based on root mean square distance (RMSD), we identified three analogues for each of the four phytochemicals by using the ZINCPharmer database. The functional descriptors selected for pharmacophore-based model search is shown in Fig. 2. Maximum number of descriptors (hydrophobic, aromatic, H bond acceptor and donor, negative ion) was attempted to cover each of the phytochemicals. The RMSD values for the small molecules against four phytochemicals, namely coronarin E, n-hexadecanoic acid, phytosphingosine and rhein, ranged within 0.235–0.683, 0.544–0.715, 0.511–0.606 and 0.044–0.435 Å, respectively (Table 2). The two-dimensional structures of the selected phytochemical analogues are shown in Fig. 3.

Fig. 2

Functional descriptors of identified phytochemicals from Amomum subulatum rhizome used as pharmacophores as generated by PubChem and ZINCPharmer

Table 2

ZINC small molecules screened based on phytochemical pharmacophores

Phytochemical	Small molecules based on phytochemical pharmacophores
	Analogues (ZINC id)	Chemical name of the analogues	RMSD (Å)	Mass
Coronarin E (CID 9,971,144)	ZINC43742965	[(2R,4aR,5R,8aS)-1,1,4a-trimethyl-6-methylidene-5-[(2-oxochromen-7-yl)oxymethyl]-3,4,5,7,8,8a-hexahydro-2H-naphthalen-2-yl] acetate (polyanthin)	0.235	425
	ZINC12231395	8-(3,5-dimethylpyrazol-1-yl)-3,7-dimethyl-1-(2-phenylethyl)purine-2,6-dione	0.630	378
	ZINC09794889	(2S)-N-benzhydryl-2-[[4-(2-fluorophenyl)-1,2,4-triazol-3-yl]sulfanyl]propanamide	0.683	433
n-Hexadecanoic acid (PCID 985)	ZINC13368240	2-(4-fluorophenoxy)-N-[4-[[2-(2-methylpropyl)-5-oxo-[1,3,4]thiadiazolo[3,2-a]pyrimidin-7-yl]methoxy]phenyl]acetamide	0.675	483
	ZINC13368607	2-(3,4-dimethylphenoxy)-N-[4-[(5-oxo-2-propyl-[1,3,4]thiadiazolo[3,2-a]pyrimidin-7-yl)methoxy]phenyl]acetamide	0.715	479
	ZINC08603481	(E)-2-cyano-N-(4-fluorophenyl)-3-[4-[(2-methyl-5-oxo-[1,3,4]thiadiazolo[3,2-a]pyrimidin-7-yl)methoxy]phenyl]prop-2-enamide)	0.544	461
Phytosphingosine (CID 122,121)	ZINC08781724	(5S)-N'-[3-chloro-5-(trifluoromethyl)pyridin-2-yl]-N'-methyl-3-[(3S)-3-methyl-5-oxo-1-phenyl-1,2,4-triazolidin-3-yl]-4,5-dihydro-1,2-oxazole-5-carbohydrazide)	0.511	498
	ZINC08964497	Ethyl (3S)-3-[(4aR,6aS,6bS,7aR)-2-acetoxy-4a,6a-dimethyl-1,2,3,4,4a,4b,5,6,6a,7a,8,8a,8b,9-tetradecahydro-6bH-naphtho[2',1':4,5]indeno[1,2-b]oxiren-6b-yl]-3-hydroxybutanoate	0.606	461
	ZINC59385771	(4S)-4-(2,4-dihydroxyphenyl)-N'-{(E)-[4-(dimethylamino)phenyl]methylene}-6-methyl-2-oxo-1,2,3,4-tetrahydro-5-pyrimidinecarbohydrazide	0.536	409
Rhein (CID 10,168)	ZINC03875972	1,4-dihydroxy-9,10-dioxoanthracene-2-sulfonic acid (rufianic acid)	0.435	319
	ZINC03977762	4,5,7-trihydroxy-9,10-dioxoanthracene-2-carboxylic acid (emodic acid)	0.044	299
	ZINC34053106	5-acetyloxy-4-hydroxy-9,10-dioxoanthracene-2-carboxylic acid	0.044	325

Fig. 3

2D structures of selected phytochemicals and associated ZINC analogues

Evaluation of Drug Likeliness

SwissADME is an open-source server and predicts absorption, distribution, metabolism and excretion properties of candidate drugs through the parameters such as physicochemical descriptors, pharmacokinetic properties, drug-like nature and medicinal chemistry friendliness [29,62]. The result was represented in a ‘Boiled Egg’ graphical interface, see Fig. 4. Boiled Egg graph plotted TPSA (total polar surface area) and WLogP (lipophilicity) in the x and y axis respectively. TPSA indicates drug properties such as gastro-intestinal absorption (GIA) and brain permeability (BP), while WLogP showed lipophilicity of those compounds. Except for two compounds, all other sixteen (16) molecules showed high gastro-intestinal absorption, and coronarin E analogue ZINC09794889 was permeable in the brain tissue. These parameters (GIA and BP) are associated with physiological drug metabolism and submitted candidates, including controls found to pass either one or both the criteria.

Fig. 4

Evaluation of drug-like properties of ligands, represented by Boiled-Egg graph; yellow yolk: blood–brain barrier permeate; white, gastro-intestinal absorption

Furthermore, other drug likeliness criteria values such as TPSA, iLogP and XlogP3, ESOL Class and ESOL Log S and p-glycoprotein substrate (PGP) and Lipinski’s violations are presented in Table 3. The TPSA values should be ranged between 20 and 130 Å [63]. iLogP and XlogP3 both are the determinants of lipophilicity. iLOGP is based on GB/SA (Generalized-Born and solvent accessible surface area) model which relies on free solvation energies in n-octanol and water. XlogP3 on the other hand is automistic method with correction factors. All the compounds were found well within the lipophilicity range (iLOGP: 0.19 to 4.37 and XlogP3: 1.48 to 6.19). All the sixteen candidate drugs and two controls were found to comply more or less within the range. Log S depicts the water solubility of the drug, and compounds having values less than 0 are highly water-soluble, as determined in our study for all submitted molecules. PGP is a membrane-bound transporter and actively removes the drug molecules from the cells through an active concentration gradient. Therefore, PGP + compounds may not be bioavailable to the cells [64]. Among eighteen compounds, only four were found to be PGP + . However, these compounds passed all the other criteria. Finally, all the molecules are found to be well within the limits of Lipinski criteria violation (violations 0–1) [65].

Table 3

Evaluation of drug likeliness by SwissADME

Molecule	TPSA	iLOGP	XLOGP3	ESOL log S	ESOL class	GI absorption	BBB permeate	Pgp substrate	Lipinski #violations
Kaempferol	111.13	1.7	1.9	− 3.31	Soluble	High	No	No	0
Ciglitazone	80.7	2.82	4.72	− 4.74	Moderately soluble	High	No	No	0
Coronarin E	13.14	3.86	6.19	− 5.55	Moderately soluble	Low	No	No	1
ZINC43742965	65.74	3.86	5.52	− 5.86	Moderately soluble	High	No	No	1
ZINC12231395	79.64	3.72	2.63	− 4.11	Moderately soluble	High	No	No	0
ZINC09794889	85.11	3.27	4.97	− 5.67	Moderately soluble	High	No	Yes	1
n-Hexadecenoic acid	37.3	3.85	7.17	− 5.02	Moderately soluble	High	Yes	No	1
ZINC13368240	123.06	3.63	3.94	− 5.11	Moderately soluble	High	No	No	0
ZINC13368607	123.06	3.94	4.13	− 5.21	Moderately soluble	High	No	No	0
ZINC08603481	137.62	3.41	3.54	− 4.94	Moderately soluble	Low	No	No	0
Phytosphingosine	86.71	4.23	4.64	− 3.68	Soluble	High	No	Yes	0
ZINC08781724	111.19	2.73	3.07	− 4.66	Moderately soluble	High	No	Yes	0
ZINC08964497	85.36	4.37	3.53	− 4.46	Moderately soluble	High	No	No	0
ZINC59385771	126.29	2.06	1.48	− 3.21	Soluble	High	No	Yes	0
Rhein	111.9	1.37	2.23	− 3.36	Soluble	High	No	No	0
ZINC03875972	140.18	0.19	2.01	− 3.42	Soluble	High	No	No	0
ZINC03977762	134.96	0.93	1.88	− 3.22	Soluble	High	No	No	0
ZINC34053106	120.8	1.79	2.34	− 3.5	Soluble	High	No	No	0

GI, gastro-intestinal; BBB, blood–brain barrier

Molecular Docking

Molecular docking is considered a promising tool to obtain a fast and preliminary understanding of the drug-like interaction properties of the compounds against target proteins. Atabaki et al. (2021) [66] showed through molecular docking experiments that phytochemical compounds from Jurinea macrocephala could inhibit HPV18 E6 oncoprotein. Furthermore, it was reported in an in silico study that the phytochemical quinizarin could prevent oncogenesis through inhibition of anti-apoptotic protein BCL2 [67]. DockThor was used here for the molecular docking experiment, and resultant binding affinities (kcal mol−1) lower than the control were considered optimal interaction values. DockThor was substantially used in the literature, Singh et al. (2021) [68] used this software to show the inhibitory potentials of phytochemicals of Schefflera vinosa extract against the target protein histone deacetylase (HDAC) of rice pathogen Magnaporthe oryzae. Recently, DockThor was used substantially to establish anti-viral properties of various medicinal plants including our previous study on SARS CoV2 [69-71]. In this study, we found that the phytochemical n-hexadecanoic acid and its analogues (ZINC13368607) had the overall high efficacy against most of the target proteins, when compared with the controls. n-Hexadecanoic acid analogue ZINC13368607 scored the highest against both the proteins XIAP and LIVIN (− 9.158 and − 8.130 kcal mol−1). Phytosphingosine analogue ZINC08964497, n-hexadecanoic acid-based compounds ZINC08603481 and ZINC13368607 showed the highest scores against BCL2, E6 and E7, respectively (Table 4). In most of the phytochemical analogue interactions, we found common amino acids compared to that of control (Table 5). Overall, it was noted that phytochemical analogues performed better than the original phytochemicals. The details of amino acid interactions of top-ranked ligands are presented in Fig. 5.

Table 4

Binding affinities (kcal mol.−1) of phytochemicals and subsequent pharmacophore modelled compounds with target proteins interactions as obtained by DockThor

Compounds		Binding affinity (kcal mol−1)
		XIAP	LIVIN	BCL2	E6	E7
Natural control: kaempferol (CID 5,280,863)		− 8.252	− 7.470	− 7.664	− 8.389	− 7.595
Synthetic control: ciglitazone (CID 2750)		− 8.055	− 7.477	− 8.493	− 7.246	− 7.854
Coronarin E (CID 9,971,144)		− 7.948	− 8.085	− 9.148	− 8.101	− 8.626
Coronarin E analogues	ZINC43742965	− 8.319	− 7.615	− 8.659	− 8.587	− 8.099
	ZINC12231395	− 8.720	− 7.952	− 7.525	− 8.013	− 8.584
	ZINC09794889	− 8.216	− 7.931	− 7.877	− 8.790	− 7.679
n − Hexadecanoic acid (CID 985)		− 7.353	− 6.411	− 6.429	− 7.855	− 6.616
n − Hexadecanoic acid analogues	ZINC13368240	− 8.496	− 7.840	− 8.752	− 9.501	− 7.906
	ZINC13368607	− 9.158	− 8.130	− 9.054	− 7.892	− 8.992
	ZINC08603481	− 9.022	− 7.749	− 9.009	− 9.537	− 8.459
Phytosphingosine (CID 122,121)		− 7.878	− 7.047	− 7.410	− 7.176	− 7.428
Phytosphingosine analogues	ZINC08781724	− 8.352	− 7.436	− 7.604	− 8.055	− 7.315
	ZINC08964497	− 8.428	− 7.318	− 9.563	− 8.452	− 8.645
	ZINC59385771	− 7.694	− 7.458	− 7.471	− 8.799	− 7.074
Rhein (CID 10,168)		− 6.434	− 7.425	− 6.682	− 8.008	− 7.679
Rhein analogues	ZINC03875972	− 6.812	− 7.320	− 6.660	− 8.129	− 7.445
	ZINC03977762	− 6.551	− 7.226	− 6.527	− 8.045	− 7.663
	ZINC34053106	− 7.183	− 7.099	− 6.925	− 7.833	− 6.761

Bold, binding affinities higher than controls; Underline, best rank ligand for each of the proteins

Table 5

Amino acid interaction between the best ranked ligands and proteins

Compounds	Proteins	Amino acid interaction
Kaempferol	XIAP	GLU75A, GLU79A, LYS83A, TRP84A
Ciglitazone		GLU75A, GLU79A, LYS83A
ZINC13368607		TRP71A, LYS72A, GLU75A, LYS83A, TRP84A, LEU109A
Kaempferol	LIVIN	ALA7A, LEU10A, HIS11A, ARG29A, ASP43A,
Ciglitazone		LEU10A, HIS11A, PRO25A, ASP43A, ALA47A,
ZINC13368607		MET1A, LYS8A, ALA7A, LEU10A, ARG46A, LYS135A, ARG136A
Kaempferol	BCL2	GLU66A, ASP71A, GLY72A, HIS115A
Ciglitazone		ASP34A, ARG38A, TYR133A, SER136A
ZINC08964497		PHE35A, ARG38A, GLY76A, TYR133A
Kaempferol	E6	LEU49C, CYS50C, VAL61C, SER73C
Ciglitazone		CYS50C, VAL61C, ARG130C
ZINC08603481		ARG7C, ARG9C, VAL30C, PHE44C, LEU49C, CYS50C, VAL52C, TYR58C, VAL61C, LEU66C, ARG130C
Kaempferol	E7	ASP13B, GLU31B, LEU34B, MET35B
Ciglitazone		ARG17A, LEU18B, CYS19B, PRO43A, SER46A
ZINC13368607		HIS24B, ILE27B, ARG28B, MET35B

Bold, common interacting amino acids of natural control kaempferol and ligands; Underline, common interacting amino acids of synthetic control ciglitazone and ligands

Fig. 5

Interacting amino acid between the best ranked ligands and proteins

Principal Component Analysis and Molecular Pharmacophore Alignment

PCA is a multivariate statistical tool used to understand the correlation among different variables and witnessed diverse applications in scientific research. Hajdari et al. (2020) [72] grouped eight plant species of Lamiaceae into four principal components based on their volatile constituents. Recently, Meftahizadeh et al. (2021) [73] showed that the Dalgan genotype of Hibiscus sabdariffa had the most suitable morphological and phytochemical traits among other genotypes. In our previous study on the dengue virus, we applied the PCA tool to substantiate the beneficial role of piperine to inhibit dengue proteins [74]. The comparative docking results (binding affinity, kcal mol−1) of 18 ligands against subsequent target proteins are presented in Fig. 5A as a heat map, and PCA was performed with these docking scores. The first (PC1) and second (PC2) components explained 65.20 and 17.0% of the variance. We observed four clusters in the PCA. While cluster 1 comprised four compounds, namely ZINC13368240, ZINC08603481, ZINC09794889 and ZINC43742965, cluster 2 had five compounds ZINC59385771, ZINC03875972, ZINC34053106, n-hexadecanoic acid and natural control kaempferol. Cluster 3 represented four compounds ZINC08964497, ZINC12231395, ZINC13368607 and coronarin E, and cluster 4 showed five ligands ZINC08781724, ZINC03977762, rhein, phytosphingosine and synthetic control ciglitazone, see Fig. 5B. PCA clusters 1 and 3 were found to have maximum top-ranking ligands (high binding affinity) (Fig. 6). With these two clusters, molecular alignment was performed by PharmaGist, and in case of both clusters 1 and 3, we have observed three common chemical features (1 aromatic and 2 H bond acceptors and one aromatic, 1 H bond acceptor and 1 hydrophobic residue, respectively), see Fig. 7. Details of the alignment are shown in Figs. S6 and S7.

Fig. 6

a Heat map binding affinity (kcal mol−1) of ligands and proteins. b Principal component analysis of docking results (kcal mol.−1)

Fig. 7

Molecular alignment based on pharmacophores, as determined by PharmaGist. a Common aligned functional groups (dotted line) of PCA cluster 1 and (A1–A4): PCA cluster 1 compounds showing aligned functional groups (dotted line). b Common aligned functional groups (dotted line) of PCA cluster 3 and (B1-B4): PCA cluster 3 compounds showing aligned functional groups (dotted line)

Protein Flexibility-Molecular Dynamic (MD) Simulation

Evaluation of protein flexibility post ligand was deemed necessary in our study to evaluate the structural stability of protein–ligand complex after docking interaction. CABS-flex 2.0 server was used to evaluate the protein–ligand complex stability in this study. Tumskiy et al. (2021) [75] applied CABS-flex 2.0 server to study the stability of SARS CoV 2 main protease in complex with a few novel inhibitors. Furthermore, Dey et al. (2021) [76] showed the stability of angiotensin-converting enzyme 2 (ACE2) receptor-severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) S protein by using this software. We selected the protein wise best ranking ligands (XIAP-ZINC13368607, LIVIN-ZINC13368607, BCL2-ZINC08964497 and E6-ZINC08603481) for the dynamic simulation study (Fig. 8). For all the complexes, we observed minimum fluctuations of residues for all of the complexes. The RMSF variation for all the complexes was found as less than 10 Å. Therefore, all the protein–ligand complexes evaluated were found to be stable in the physiological condition.

Fig. 8

The RMSF profiles of protein–ligand complexes as obtained by CABS-flex 2.0

Conclusion

We identified four compounds, namely coronarin E, n-hexadecanoic acid, phytosphingosine and rhein, from the rhizome extracts of Amomum subulatum (greater cardamom) screened against HPV-induced cervical cancer proteins by using in silico molecular docking approach. Furthermore, based on the pharmacophore modelling, three analogous molecules for each of the phytochemicals were identified. All the compounds were found to possess adequate drug-like properties. Overall, n-hexadecanoic acid and its analogues were found to efficiently inhibit target oncoproteins, namely HPVE6, HPVE7, XIAP, BCL2 and LIVIN, among sixteen molecules screened. Statistical evaluation based on principal component analysis (PCA) showed four phytochemical clusters including controls (kaempferol and ciglitazone). Pharmacophore alignment further revealed that majorly hydrophobic, aromatic and H-bond acceptor groups explained top-ranked ligands’ oncoprotein inhibitory potentials. Finally, it can be concluded that large cardamom possesses considerable therapeutic potential against HPV-induced cervical cancer; however, validation in the wet lab condition is warranted.

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 773 KB)