There are several reports on the anti-cancer activity of zerumbone, however, the reports lack the actual protein target. Moreover, the study on its analogues is the least concern and there are a very few report on the anti-cancer activity of its analogs. This investigation is an attempt to determine the inhibitory action of zerumbone against TNF-α release in suppression of cervical cancer cell proliferation (HeLa). Further, the investigation also focused on the computational studies on its analogs as the inhibitor of TNF-alpha using molecular docking, molecular dynamics (MD) simulation and density functional theory (DFT) approaches. The outcome of the present study reveals that the zerumbone analogues may be a potent anti-cancer compound as compared to zerumbone. From the molecular docking analysis, it was observed and substantiated that the alpha-beta unsaturated carbonyl scaffold is the main driving force for its anti-cancer activity.
Keywords: TNF-α, docking, MD simulation, DFT, ADME-Toxicity
Introduction
Cervical cancer accounts for the fourth-most common cause of cancer and the fourth most common cause of death from cancer in women.1 90 % of cervical cancer is caused by Human papillomavirus (HPV) infection and around 70 % of cervical cancers have occurred in developing countries.2,3 Cervical cancer is also the most common cause of cancer death in developing countries with the low-income group.4 The tumor necrosis factor alpha (TNF-α) is a vital cytokine involved in inflammation, immunity, and cellular organization.5 It is associated with the promotion and proliferation of cervical cancer cells.6 It is also associated with an increased risk for the development of invasive cervical cancer (ICC).7 Anti-TNF therapy or inhibition of TNF-α has shown some meaningful results in cancer therapy. There are reports of certain biopharmaceuticals acting on TNF-α for the treatment of cervical cancer, renal cell carcinoma and breast cancer.8 On the other hand, zerumbone, a class of sesquiterpenoid which is found in the subtropical ginger Zingiber zerumbet Smith is known for its distinct anti-inflammatory and anti-carcinogenesis activities.9-11. There are also reports of Z. zerumbet traditionally utilized in Southeast Asia for various purposes such as anti-inflammatory medicine, chemoprevention and chemotherapy strategies.12
Recent in vitro and in vivo studies have revealed that zerumbone exhibit anti-proliferative effect on human cervical cancer (HeLa) cell line.13 There are also reports which suggest that zerumbone used in combination with cisplatin were able to suppress cervical cancer in BALB/c mice.14 There are also histological reports which revealed that zerumbone inhibited the cervical dysplasia from developing into more severe dysplasia.15 In fact, there are several reports on various important proteins that can be inhibited by zerumbone for controlling cancer cell growth.16-18
There are several reports on zerumbone however, the protein targets, binding pocket and the molecular level events which might have occurred between zerumbone and the target protein has not been fully studied or understood properly. Nevertheless, the molecular interaction between an enzyme’s active site and a ligand can be investigated using molecular docking tools and techniques. In fact, docking tools are helpful in understanding the insights of the molecular events that occur in the binding site of a target protein. The importance of molecular docking helps in supplementing the experimentally obtained data. In this investigation, using Molegro Virtual Docker (MVD 6.0), molecular dockings studies were set up to evaluate the binding mechanism of zerumbone and its analogs with the target enzyme TNF-α. Further, the protein-ligand complex of the top dockings poses was carried forward for molecular dynamics simulation for 20 nanoseconds using Gromacs 5.02. Lastly, ADME-Toxicity (Absorption, Distribution, Metabolism, Excretion and Toxicity) parameters were calculated for the top docking hits followed by the DFT study of the docking hits.
Results and discussion
The compound zerumbone (2,6,9,9-tetramethylcycloundeca-2,6,10-trien-1-one) contains one active moiety, alpha-beta unsaturated carbonyl groups. Due to the presence of this moiety zerumbone is found to be active against various cancer cells leading to apoptosis19. Besides zerumbone, the rhizomes of Z. zerumbet also contain humulene and camprene which lacks this particular scaffold. It should be noted that the alpha-beta unsaturated carbonyl group is a notable moiety of zerumbone which play a vital role in interactions with unidentified drug targets19,20. Furthermore, Murakami et al. reported that alpha humulene, which lacks the alpha-beta unsaturated carbonyl group has not been able to disrupt 12-O-Tetradecanoylphorbol-13-acetate (TPA) which is a tumour promoter.19 Hence, these two compounds were consistently inactive compared with zerumbone20. There are several reports of compounds that have the alpha-beta unsaturated carbonyl group such as zerumbone which exhibit versatile biological activities including inhibition of tumors cell growth, induction of differentiation, apoptosis, heat shock protein etc19-21. Thus, these reports indicated that the alpha-beta unsaturated carbonyl group is the driving force to act against various cancer cell lines and protein targets.19
The cell line assay of zerumbone revealed the cell death of HeLa cells after 24 h of zerumbone treatment (Fig. 1). In fact, the tumor necrosis factor (TNF) α protein plays a dual role in promoting and inhibiting cancer depending largely on the pathway initiated by the binding of the protein to its receptor22. Excessive release of the TNF-α protein play a role in a number of pathological conditions, inflammations and cancer. There is release of TNF-α and other proteins of TNF pathway in cancerous cells, therefore, blocking the release of this protein can lead to inhibition of cell proliferation. Our study showed that 5 µM of zerumbone was enough to inhibit the release of TNF-α in cervical cancer (HeLa) cells (Fig. S1 of supplementary material). Takada et al.17 in their experiments shows that 25 µM of zerumbone were enough to suppress TNF-α induced activation of NF-κB in human lung carcinoma (H1299) cell lines. Takada et al.17 had pointed out that the constitutive expression of the protein can be blocked by zerumbone. They found that zerumbone had no direct effect on the NF-κB inhibition but it can inhibit TNF induced activation in human lung carcinoma cell lines (H1299) pre-treated with zerumbone. Rivas et al.23 further proposed that TNF-α and NF-κB inhibition was essential for combating breast cancer. Their experiments revealed that the blocking of tumour necrosis factor receptor 1 (TNFR1) or TNFR2 with specific antibodies at concentrations of 2–10 µM were able to impair the TNF-α signaling as a breast cancer tumor promoter. Zerumbone was active due to the presence of alpha-beta unsaturated carbonyl groups in its structure and zerumbone analogs which possessed this functional group were carried forward for computational studies which include molecular docking and MD simulation studies.
The docking scores of the analogs which were sorted based on MolDock score, Rerank score, and interaction energy are presented in Table 1. The 2D structure of the docked compounds at the active site of TNF-α is shown in Table ST1 of supplementary material. From, Table 1 it is observed that the analogs bind at the active site of TNF-α with favourable docking scores which represent favourable binding affinity. Five analogs were found to be more potent than zerumbone which bind at the active site of TNF-α based on the docking score and binding affinity (Table 1 and Fig 2). CID101417968, CID16086592, and CID102074066 were among the top docking hits which have effective docking scores and interaction energy than zerumbone. The binding affinity of the studied compounds also suggest that CID101417968 possessed a stronger binding affinity of -45.88 kJ/mol compared to -22.88 kJ/mol of CID16086592 and -19.54 kJ/mol of zerumbone (Fig. 2). Moreover, the target enzyme studied in this investigation is a signaling pathway that interferes with the TNF-α promotion and proliferation of cervical cancer cells. Among these potent compounds, CID101417968 and CID16086592 formed hydrogen bonding interaction at the active site of the enzyme (Table 2 and Fig S3 of supplementary material) suggesting a stronger binding affinity than zerumbone. While zerumbone showed hydrophobic and steric interaction with Tyr59, Ser60, Tyr119, Leu120, Gly121 and Gly122 (Fig 3C). However, no hydrogen bond interaction was observed between TNF-α and zerumbone. The overall steric and hydrophobic interaction of the top docking hits is shown in Fig. 3.
Earlier, Fatima et al. reported a similar molecular docking study which suggests that zerumbone binds at the same binding pocket.22 They also reported that zerumbone exhibited high binding affinity from molecular docking studies with TNF-α which propose its favorable binding affinity.22 Furthermore, the X-Ray diffraction experimental study conducted by He et al24 reported that an experimental ligand Compound-307 (CID5327044) also forms molecular interaction with Leu57, Ile58, Tyr59, Ser60, Gln61, Tyr119, Leu120, Gly121, Gly122 and Tyr151 at the binding pocket of the enzyme.25 Interestingly, the potent docking hits also bind at the same active site of binding cavity. Surprisingly, the top docking hit CID101417968, showed a similar molecular interaction with Compound-307 forming similar steric/hydrophobic interactions with Leu57, Ile58, Tyr59, Ser60, Gln61, Tyr119, Leu120, Gly121 and Tyr151. Whereas CID16086592 also formed a similar interaction with CID101417968 except for Gln61 (Fig 3A and B). Therefore, the analogs of zerumbone are likely to possess strong anti-cancer property since they have 90% structural similarity with zerumbone. Moreover, these compounds also possessed the alpha-beta unsaturated carbonyl groups and they bound at the same binding site as that of zerumbone. Additionally, Fig. S4 of supplementary material depicts the electrostatic interaction map of the docking hits at the binding site of the TNF-α. It is evident from the electrostatic map that these compounds possessed strong electrostatic interaction. Based on the available reports on TNF-α, there is no experimental data available for zerumbone binding to TNF-α.25 Hence we could assume that the there is a fair affinity between the analogs and TNF-α. Our docking results highlight that the analogs of zerumbone can possibly inhibit the enzyme better than zerumbone. Nevertheless, the binding of the zerumbone at the active site of TNF-α is also consistent with the data reported by of Fatima et al.22
The molecular dynamics simulation analysis depicts the convergence of the root mean squared deviation (RMSD) backbone of the protein-ligand complex within 20 ns MD simulation production (Fig 4A). The MD simulation observed a small deviation depicting a stable conformation which indicates the stable system during the 20 ns MD simulation run. The average drift of the 2AZ5-CID101417968 complex off to 0.18 nm, 2AZ5-CID16086593 at 0.23 nm, 2AZ5-zerumbone at 0.24 nm and the TNF-α enzyme at 0.22 nm. In the radius of gyration (Rg) analysis some minor drifts were observed in the 2AZ5-zerumbone complex in the 4.8 ns, 10.5, 14.7 ns and 14.9 ns respectively. But in the top docking hits, the complex looks very stable without any major drifts and deviations (Fig 4 B). The curves of solvent accessible surface area (SASA) indicate that all the systems both the protein and the protein-ligand complex are stable during the 20 ns MD simulation run. The SASA calculates the hydrophobic properties of molecules and their match or mismatch in receptor–ligand complexes. In our case, all the systems look stable as evidence from Fig 5C. Overall, the RMSD, SASA and Rg analysis, correspond to the relative global account of the general tertiary structure of the TNF-α enzyme and the protein-ligand complex which represent the system are stable during the studied simulations (Fig. 4A, B, and C). The root mean squared fluctuation (RMSF) which highlights the flexible regions of the system indicates that the values higher than 0.25 nm are characteristics of amino acids residues belonging to flexible regions (Fig 5A). Figure 5B shows that the patterns of the backbone hydrogen bonds formed are closely related except for CID101417968 docked complex which showed additional contacts.
The ADME-Toxicity analysis of the compounds indicate that the analogs of zerumbone fall within the permissible range of the drugs available in the market as evident from predicted ADME-Toxicity servers (Table 3). The investigated analogs of zerumbone have been classified as a mid-structure drug without any violation of Lipinski rule of five and qualified for CMC like rule. The ADME-Toxicity analysis also revealed that the top docking hits have more or less similar pharmacokinetics behaviour with zerumbone. Thus, the analogs of zerumbone can target the TNF pathway and the TNF enzyme for inhibition with similar binding affinity.
Finally, in the DFT study, the frontier orbital energies, viz., HOMO and LUMO were used to calculate the band gap energy which was subsequently employed to ascertain the strength and stability of the molecular interactions (Fig 6). The band gap energy of the top dockings hits was calculated which indicates a low band energy gap. In fact, a low band energy gap is an indication of a reactive compound while a wide energy gap implies that the activity is not sufficient at the active site of a protein receptor.26 The band gap energy of CID101417968 and CID16086592 is lower than zerumbone which imply a reactive nature at the active site of the TNF-α. Thus indicating that these docked compounds will definitely have a strong binding affinity than zerumbone.
Conclusions
As a concluding remark, the authors would like to point that analogs of zerumbone might be potent anti-cancer agents targeting the TNF-α enzyme. In addition, the analogs also have the active site scaffold, i.e. alpha-beta unsaturated carbonyl groups which are responsible in zerumbone for apoptosis of cancer cells. Additionally, the molecular interaction analysis reveals that these docking hits have similar binding mode with the experimental compound-307 and their conformational stability are confirmed by MD simulations. Nevertheless, these results require further investigation and hence conclude that these zerumbone analogs may be promising anti-cancer agents.
Experimental
Compound
The plant materials were collected from Moreh, Manipur, North-East India, 286 m from sea level, longitude 94.34217° E and latitude 24.35172° N in March 2012. The plant was identified by the taxonomist of the institute and had given the accession number as IBSD/Z-42-23. Fresh rhizomes were collected and washed thoroughly with tap water. These were cut into 5–6 mm slices and put into the Clevenger type oil extractor. Zerumbone was isolated from this oil and identified by using FT-IR, Proton and 13 C NMR and mass spectrophotometer. This isolated zerumbone was compared with the standard found in the market and was used for further studies.
Cell Cultures and TNF-α assay
All the cells were purchased from ATCC (Manassas, VA). HeLa cells were cultured in DMEM medium supplemented with 10 % FBS having 100 units/ml penicillin and 100 μg/ml streptomycin solution in a humidified atmosphere of 5 % CO2 at 37°C. All the culture reagents were purchased from Invitrogen (Carlsbad, CA). Cells were seeded in 24-well tissue culture plates, at a density of 5 × 104 cells/well and incubated for 16 h to adhere the cells. Then the cells were treated with DMSO, cisplatin (Sigma Aldrich, USA) or increasing concentrations of Zerumbone (5, 10 and 20 μM) for 24 h. After 24 h, the media was collected and stored at -20oC until use. Zerumbone was dissolved in DMSO at a stock concentration of 10 µM. Cisplatin was dissolved in water at a stock concentration of 25ug/mL and cells were treated with a concentration of 50 µg/mL. Collected media along with cells were taken and centrifuged at 5000 rpm for 15 minutes. The pellet was discarded; the supernatant was used for TNF-α assay. The TNF-α level was determined using a commercially available ELISA kit (Millipore, USA) according to the instruction of the manufacturer. A range of zerumbone concentrations (5, 10 and 20 μM) was tested for release TNF-α level in the cervical cancer cell (HeLa). Cisplatin, a known anti-cancer drug was used as positive control.
Chemical similarity search
The analogs of zerumbone were obtained by performing a chemical similarity search. The similarity search was set at 95% Tanimoto similarity parameters.27 The geometries of zerumbone along with the retrieved compounds from the similarity search were optimized and converted to three-dimensional SYBYL mol2 file format using ChemOffice 2010.
TNF-α docking set up
The crystal structure of TNF-α (PDB Code: 2AZ5) was retrieved from the RCBS Protein Data Bank.24 The crystal structure has 148 amino acids in length and has X-Ray resolution of 2.1 Å. The binding cavity was predicted and the coordinates were set at (X: -21.69, Y: 69.17, Z: 37.41) having a restriction sphere of 14 Å radius (Fig. S1 of supplementary material). Additionally, the residues Leu57, Ile58, Tyr59, Ser60, Gln61, Tyr119, Leu120, Gly121, Gly122, Val123, Tyr151, Phe152 surrounding the vicinity of the active site pocket were set flexible.
Flexible docking of zerumbone analogs
Rigid docking is no more acceptable and hence flexible docking approached was employed for the docking simulations. In fact, the rigid molecular docking considers the molecules as objects that are rigid which do not orientate their spatial shape during the process. Nevertheless, setting all the possible conformational changes flexible and its scoring requires heavy computational cost and not recommendable. However, flexible docking which selects a small region for its conformational changes is considerable. Hence, the amino acid residues surrounding the vicinity of the binding pocket were set flexible for both the setups. The flexible side chains of the selected amino acids (Leu57, Ile58, Tyr59, Ser60, Gln61, Tyr119, Leu120, Gly121, Gly122, Val123, Tyr151, Phe152) was set with a tolerance of 1.10 having a strength of 0.80. The tolerance of a potential refers to the size of the region between a ligand atom and a receptor atom where the interaction energy is optimal. While the strength is a factor which is multiplied onto all interaction energies for the side chain such as the atomic pairwise steric interactions, hydrogen bondings, and electrostatic interactions. The 36 analogues of zerumbone retrieved from the NCBI PubChem database were imported in the Molegro Virtual Docker (MVD) and the bond flexibility of all the imported compounds was set for molecular docking simulations. The scoring function was set for MolDock Grid Score having a grid resolution of 0.30 Å with further evaluation for Internal ES (ElectroStatic) and Internal Hbond (Hydrogen bond).
Docking computation
The docking algorithm was set at a maximum iteration of 1,500 having a maximum evolution size of 50 and a minimum of 30 runs was performed for each compound. The pose generation energy threshold was set to a maximum of 30 and a minimum of 10. The simplex evolution step was set with a maximum of 300 having a neighbour distance factor of 1. The poses were ranked based on the Tabu clustering and the RMSD threshold for the multiple cluster poses was set at 2.00 Å with an energy penalty of 100. The best pose of each compound was selected for the subsequent ligand–protein interaction analysis.
Molecular docking simulation was carried out using Molegro Virtual Docker 6.0.1 which is based on a differential evolution algorithm and the MolDock scoring function. The MolDock Score is some sort of improvement which was proposed by Gehlhaar et al28 based on Piece Wise Least Potential (PLP) scoring functions. In MolDock scoring function, the Escore, is defined by the following energy terms:
Escore = Einter + E intra
Where Einter= ∑_iligand ∑_iprotein[EPLP (rij) + 332.0 (q_i q_j)/〖4rij〗^2 ]
While the PLP describes the electrostatic interactions between charged atoms. Itis a Coulomb potential with a distance-dependent dielectric constant given by D(r) = 4r.
Binding affinity calculation
The binding affinity of the docked poses was calculated using Molegro Data Modeller. It was calculated from the coefficients for the binding affinity terms which were derived using multiple linear regressions from the dataset of 200 structurally diverse complexes (taken from the PDB data bank) with known binding affinities expressed in kJ/mol. The MLR equation for calculating the binding affinity is given below.
Affinity = -19.0155 * C0 + 3.3813 * "Cofactor (hbond)" – 0.594128 * Csp2 + 0.464056 * HBond – 0.061912 * "E-Intra (vdw)" + 0.953672 * "E-Solvation" + 0.483368 * HeavyAtoms – 1.00763 * N + 3.0229 * Nplus + 1.61426 * OH – 3.10696 * OS – 3.9493 * halogen
Molecular dynamics simulation studies
All MD simulations were carried out on a Fujitsu CELCIUS R670 Workstation equipped with Intel Xeon processor having 28 GB RAM operated with Ubuntu 16.04.03 (Xenial Xerus). MD simulation for the protein and the protein-ligand complex was carried out using GROMACS version 5.1.2. For the individual enzyme, (PDB ID: 2AZ5), it was processed with OPLS-AA/L all-atom force field (2001 amino acid dihedrals) and for the protein-ligand complex, it was processed with GROMOS96 43a1 force field.29,30.
For the protein-ligand complex, initially, the coordinates of the ligands were separated using the “grep” command and the topology of the ligand were generated using the PRODRG 2.5 server31. All the MD simulations were carried out using periodic boundary conditions. Solvation was accomplished by configuring using the spc216 solvent model and the solvated structures were energy minimized using the steepest descent minimization, terminating when the maximum force is greater than 1000.0 kJmol-1nm-1. Next, the system was equilibrated with NPT ensemble phase at constant temperature (300 K) and pressure (1 bar) with a time step of 2 fs. Further, the systems were progressed with NVT simulation for 1 ns (nanoseconds), and the minimized structure was equilibrated with a timescale of 20 ns (nanoseconds). The molecular dynamics simulation was performed on both the protein and the protein–ligand binding complex for 20 ns to understand the dynamic behavior of the protein and their stability.
Analysis of the trajectories was performed using the GROMACS commands gmx –rms (root mean square deviation), -hbond, -rmsf (root mean square fluctuations), -gyrate (radius of gyration) and –sasa (solvent accessible surface area).
Assessment of ADME-Toxicity and risk parameters
Molecular docking studies usually result with the prediction of the most active bioactive compounds with favourable docking scores and binding affinity. But it does not tell the property of the compound to be effective as drug molecule. Hence, ADME-Toxicity parameters were assessed for the assessment of absorption, distribution, metabolism, and excretion (ADME). These parameters are a crucial step in the discovery process. In fact, ADME-Toxicity prediction is mostly considered in which compounds are numerous but have little or limited access to the physical samples. It is also evident that poor human pharmacokinetics (ADME–toxicity) is the main reason for the majority of the drug failure.29 BMDRC-PreADMET [https://preadmet.bmdrc.kr/] was used to compute the drug-likeness, ADME and Toxicity parameters. While SwissADME [http://www.swissadme.ch/] was used to compute the pharmacokinetics parameters.
DFT Calculations
Finally, the top docking hits were carried forward for the DFT calculations. The conformation of the top docking hits and zerumbone were exported and the molecular orbital study was analyzed using DFT calculations. Frontier molecular orbitals were calculated using program US-GAMESS.32 The calculations were set-up using “DFT/ Becke3-Lee-Yang-Parr correlation function (B3LYP) protocols with 6-31G basis set, water medium, and initial MO guess set as Huckel.33 The resulting molecular orbitals and the calculated energies have been visualized using ChemOffice. The computation for calculating the frontier orbital energies of the Highest Occupied Molecular Orbital (HOMO) and the Lowest Unoccupied Molecular Orbital Energy (LUMO) was also carried out. The band energy gap ΔELUMO-HOMO was also computed for the top docking hits. These orbital energy values play an important role in terms of electron donor and acceptor properties of a molecule and can thus be utilized to understand the reactivity of a molecule in the active site of the TNF-α enzyme.