Potent Bioactive Compounds From Seaweed Waste to Combat Cancer Through Bioinformatics Investigation

Introduction

After terrestrial plants and microorganisms, marine organisms are considered to be a prominent source of drug discovery motivation. The massive diversification of organisms in the marine ecosystem serves as a reservoir of a wide range of natural compounds (1, 2). New marine molecules and their biological potential are discovered every year, allowing these compounds to be used in functional foods and the development of new nutraceuticals and therapeutic drugs (3, 4). In Asian countries, particularly China, Korea, and Japan, marine algae have long been consumed as food, whereas in Western countries, algae are mostly used as a source of food additives (5), natural pigment (6), and minerals (7). The consumption of edible algae, as well as the establishment of algae-based companies, has recently dramatically developed (8). Macroalgae or seaweeds and microalgae are the two types of dominant photosynthetic aquatic organisms. Seaweeds (Marine algae) are a type of photosynthetic non-flowering plant-like multicellular algae found in the ocean that are rich in proteins, vitamins, minerals, and a vast number of bioactive compounds possessing significant antibacterial, anticancer, antifungal, and antiviral properties (9, 10). Around 6,000 species of seaweeds are found in deep oceans, up to 180 m deep. They are sometimes seen in shallow coastal waters. Seaweeds are classified into three major categories based on their pigmentation: red (Rhodophyta), green (Chlorophyta), and brown (Phaeophyceae; phylum: Ochrophyta) algae (11). Sulfated polysaccharides and bioactive compounds are abundant in algae that grow in saline seas or near river mouths, transforming seaweeds into potential pharmaceutical agents (12).

The increasing demand and growth in seaweed farming for seaweed consumables encourages producers to enhance their degree of seaweed processing and cultivation (13). This situation leads to an uncontrolled and uneven quality of the seaweeds' nutritive composition and pharmacological characteristics. In turn, it results in large amounts of seaweed waste usually being dumped in landfills as garbage. Industrial rejection, uneven harvest time, drying processes, and packaging waste are scattered at the processing site (14, 15). Thus, the seaweed waste can be collected and can be potentially applied in both the food and pharmaceutical industries. It has already been established and widely known that diverse groups of seaweeds contain different bioactive compounds possessing a wide range of pharmacological properties. Thus, the seaweed waste can be collected and used in pharmaceutical industries. There needs to be a way to use this potential seaweed debris for sustainable development.

Cancer is a severe issue that threatens the health of all human communities (16), and natural products have played an important role in the history of anticancer medication development (17–20). A large number of experiments in recent years have shown potent anti-cancer properties of several seaweed-derived compounds by inhibitions of tumor development, adhesion, invasion, and metastasis (10, 21, 22). Also, in recent years, much emphasis has been placed on exploring the anti-proliferative activity of bioactive compounds from edible seaweed. With this in mind, we investigate the anticancer properties of compounds from some species of edible seaweed. The phloroglucinol compound isolated from Ecklonia cava has good anti-proliferative activity in human breast cancer cells (MCF 7). Similarly, extracts from the brown algae Laminaria japonica inhibit human hepatocellular carcinoma (BEL7402) and murine leukemic cells from proliferating (P388). Fucoidan derived from the brown algae has shown anti-cancer activity in human colon cancer cells (HCT116) by activating the apoptosis signals of cancer cells and thereby inducing apoptosis. Furthermore, several seaweed compounds have been tested against different types of glioblastoma cells with promising results in suppressing cell survival while causing no adverse effects on the cells (23, 24).

Histone modifications, specifically histone acetylation/deacetylation, are an important epigenetic regulatory system that is engaged in a variety of human cancers. Histone deacetylation is a key regulator of gene transcription that involves the removal of acetyl groups. Histone deacetylases are a family of enzymes that regulate this mechanism (HDACs). HDACs are categorized into four categories: class I (HDAC 1, 2, 3, and 8), class II (IIa: HDAC 4, 5, 7, and 9; IIb: HDAC 6 and 10), class III (SIRT1, 2, 3, 4, 5, 6, and 7), and class IV (SIRT1, 2, 3, 4, 5, 6, and 7; HDAC11). HDAC dysregulation is correlated to a poor prognosis in a variety of human cancers, making these enzymes a promising therapeutic target (25, 26). HDAC inhibition has been shown to induce cell cycle arrest, tumor angiogenesis inhibition, differentiation of some transformed cell lines, and/or apoptosis in tumor cells, indicating its potential as a therapeutic target for cancer treatment (27). HDAC 2 overexpression plays a significant role in cancer progression. For example, deregulation of p53, an important tumor suppressor protein, is frequently found in tumors, and HDAC 2 deacetylation is one of the reasons, indicating its significance in specificity. Thus, developing a potent HDAC 2 selective inhibitor was useful for a potent therapeutic strategy for cancer progression and development (28).

This study implemented top–down and structure-based high-throughput screening methodologies to find novel HDAC 2 inhibitors from marine edible seaweed compounds by using the comprehensive marine natural products database (CMNPD). CMNPD is a free online database containing marine natural products (MNPs) from the ocean to facilitate drug discovery. Thousands of compounds with various physicochemical and pharmacokinetic properties, biological activity data, taxonomy, and geographic distribution of source organisms are contained in the CMNPD database. We focused on marine edible seaweed compounds to investigate novel HDAC 2 inhibitors as potential anticancer compounds. We used various in silico methods such as drug-likeness prediction, molecular docking, molecular dynamics simulation, MM-GBSA, free energy landscape, PCA, DCCM, ADME calculation, and DFT study.

Materials and Methods

Preparation of Ligands and Receptors

A total of 17 edible seaweed species were selected for our investigation (Table 1). Out of this, 193 edible seaweed compounds from the Comprehensive Marine Natural Products Database (CMNPD) (https://www.cmnpd.org) were used in this study. The Canonical SMILE format of the compounds from the database (Supplementary File S1) was retrieved, the structure was built using the structure builder option of Chimera, and geometry optimization was performed by using freely available UCSF Chimera version 1.14 software (29) before conducting experiments. The obtained structures in 3-D format were saved as.pdb files.

The crystal structure of HDAC 2 (PDB ID: 4LY1) protein was obtained from the protein data bank (PDB) repository (www.rcsb.org) with a resolution of 1.57 Å. This structure was selected from the protein data bank due to a higher resolution than the other structures in the database. The PDB file of HDAC 2 structure was processed and refined with the Protein preparation wizard [Version 2021-2, Schrodinger (30, 31)]. Ions, cofactors, and water molecules were deleted and hydrogens were added to the heavy atoms. Selenomethionines were converted to methionines. Het states were generated using Epik at pH 7.0 ± 2.0. OPLS_2005 force field (32) was used for protein minimization.

Drug-Likeliness and PAINS Filtering

SwissADME (http://www.swissadme.ch/) was used to screen potential bioactive compounds from edible seaweed by calculating the drug-likeness properties that constitute Lipinski (33, 34), Ghose (35), Veber (36), and PAINS (37) as per our previous study (38).

Molecular Docking Profiling of Seaweed Compounds

A grid was generated by taking into account the co-crystal ligand that was found in the active site of the selected protein target HDAC 2 (PDB ID: 4LY1). A grid box was generated (x = 22.19; y = −18.4; z = 0.63) at the centroid of the active site. The compounds were sequentially docked into the catalytic pocket of HDAC 2 (PDB ID: 4LY1). The glide program and its virtual screening workflow (VSW) process, including three docking protocols, namely high throughput virtual screening (HTVS), standard precision (SP) module, and extra precision (XP) module, were applied. Each ligand was docked to the crystal structure of HDAC 2 using the HTVS, resulting in one pose. Though the SP docking protocol provides a good scoring function retaining the good scoring states, about 10% of the total edible seaweed compounds were shifted from HTVS to SP, which aids in the detection of false-positive results. Moreover, about 10% of SPs total edible seaweed compounds were processed to XP, where XP provides the best scoring states (39).

Quantum Polarized Ligand Docking (QPLD)

In protein–ligand molecular docking, QPLD was used to determine the right electrical charges. QPLD is a docking method that analyzes ligand/protein interactions by combining the QM and the MM approaches. The interaction energy of the protein/newly produced ligand pose was estimated using the QSite programe in Schrödinger, with the protein being treated with the molecular mechanical (MM) approach and the ligand pose being treated with the quantum mechanical (QM) method (40). For the ligand pose within the protein environment, the Qsite algorithm created a new set of atomic partial charges. Using the Glide software with the XP-scoring function, the ligand configuration with QM-generated partial charges was re-docked to the protein (41). The same grid file that was generated by the Glide grid generating tool was utilized for QPLD, and the same ligand was selected in the ligand to be the docked option.

MD Simulations and Free Energy Landscape Analysis

The molecular dynamics simulations of the docked complex of HDAC 2 protein with the selected edible seaweed compound, cinnamyl dihydrocinnamate, were carried out using the Desmond 2020.1 from Schrödinger, LLC (42) to envisage the dynamic behavior, intermolecular interaction, and stability of the docked complex (43–45). The complex was prepared using the system builder platform by solvation with the simple point-charge (SPC) explicit water model in the orthorhombic simulation box. The solvated complex system was neutralized with a suitable number of Na⁺/Cl⁻ counter ions and a salt concentration of 0.15 M to mimic the physiological conditions (46). The receptor–ligand complex system was designated with the OPLS-2005 force field, and an explicit solvent model with the SPC water molecules was used in this system (47–51). The HDAC 2 (protein)-cinnamyl dihydrocinnamate (ligand) complex system was initially equilibrated using an NVT ensemble for 100 ns. Following the previous phase, a short run of equilibration and minimization was carried out using an NPT ensemble for 12 ns. The NPT ensemble was set up using the Nose-Hoover chain coupling scheme (52, 53) with a temperature of 37°C, a relaxation duration of 1.0 ps, and a pressure of 1 bar. A time step of 2 fs was chosen. With a relaxation time of 2 ps, the Martyna–Tuckerman–Klein chain coupling system (54) barostat method was used for pressure control. The particle mesh Ewald technique (55) was used to calculate long-range electrostatic interactions, with the radius for Coulomb interactions set at 9Å (55). The bonded forces were calculated utilizing the RESPA integrator with a time step of 2 fs for each trajectory. Finally, to monitor the binding stability of the protein–ligand docked complex, the root mean square deviation (RMSD), the radius of gyration (Rg), root mean square fluctuation (RMSF), the number of hydrogen (H-bonds), plot of ligand interaction in the binding cavity, stacked bar chart plot of protein–ligand contact analysis, total contacts timeline analysis of MD trajectory, and analysis of the torsional degree of freedom for the rotatable bonds of the selected ligand were calculated by the results obtained after MD simulations of protein–ligand docked complex of 100 ns of simulation time.

The free energy landscape of protein folding on HDAC 2 protein–cinnamyl dihydrocinnamate complex was measured using geo_measures v 0.8 (56). The MD trajectory vs. RMSD and radius of gyration (Rg) energy profile of protein folding was recorded in a 3-dimensional plot using the matplotlib python package utilizing Geo_measures, which includes a sophisticated library of g_sham.

Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) Calculations

The binding free energies of the HDAC 2 protein (receptor) and ligands viz., cinnamyl dihydrocinnamate complex were investigated using the molecular mechanics-generalized born surface area (MM-GBSA) methodology using the prime module's Python script thermal mmgbsa.py (57, 58). OPLS 2005 force field, VSGB solvent model, and rotamer search methods were used to calculate the binding free energy (59). The binding free energy calculation upon binding of ligands with the receptor HDAC 2 was calculated using the following equations:

where ΔG_bind = binding free energy, G_complex = free energy of the complex, G_protein = free energy of the target protein, and G_ligand = free energy of the ligand.

Principal Component Analysis (PCA) and Dynamic Cross-Correlation Matrix (DCCM) of the MD Simulation Trajectories

Principal component analysis analysis was performed to recover the global movements of the trajectories during a 100-ns simulation of cinnamyl dihydrocinnamate complex with HDAC 2. A covariance matrix was generated to calculate the PCA. The movements of trajectories of 10 alternative conformational modes of the main component were calculated for conformational analysis of the cinnamyl dihydrocinnamate in a complex bound with HDAC 2 protein, and a comparison of the first highest mode (PC 2) was investigated. To investigate domain correlations, the dynamic cross-correlation matrix (DCCM) was built across all Cα-atoms during a 100-ns MD simulation for HDAC 2-compound bound complex. The DCCM was carried out using Schrodinger's trj essential dynamics script (30).

Density Functional Theory (DFT) Calculations

The Jaguar module in Schrödinger was used to perform DFT calculations. The basis set 6-31G^* was used to do a complete geometry optimization using the Hybrid DFT model with Becke's three-parameter exchange potential and Lee-Yang-Parr's correlation functional (B3LYP) technique. DFT energetics depicts the 3D electronic states of molecules to assess the transfer of electrons and reactivity in a specific environment. The quantum chemical properties of the promising compound, like molecular orbital (HOMO), the lowest unoccupied molecular orbital (LUMO), and the energy gap were calculated along with the molecular electrostatic potential surface [MESP; (60, 61)].

In silico ADMET Property Prediction of the Compound

Using the Schrodinger suite's QikProp module, the absorption, distribution, metabolism, and excretion (ADME) properties of the chosen compound were determined (62, 63). SASA, FOSA, FISA, PISA, WPSA, QPlogPo/w, QPlogS, QPlogkhsa, QplogBB, QPPCaco percent human oral absorption, rule of five, and rule of three were compared to a manual-recommended range of values for ADME attributes (62). QikProp's descriptors were used to make ADMET predictions and calculate drug-likeness parameter calculations.

The toxicity profile of the selected bioactive compound, cinnamyl dihydrocinnamate, from edible seaweed waste was predicted using pkCSM online webserver (http://biosig.unimelb.edu.au/pkcsm/) (64). The selected compound structure was drawn using the Chemdraw drawing tool and was converted into SMILES format for the prediction studies. The possible predicted toxicity risks evaluated were AMES toxicity, hepatotoxicity, skin irritation, max. tolerated dose (human), oral rat acute toxicity (LD50), oral rat chronic toxicity (LOAEL), and hERG I, II inhibitors (Figure 1).

Result and Discussion

Drug Likeness Virtual Screening Analysis

The edible seaweed compounds (n = 193) extracted from the CMPND database were predicted to evaluate the pharmacokinetics and drug-likeness using the Swiss ADME online webserver (65). The Swiss ADME results give a detailed and extensive physicochemical profile of the edible seaweed compounds. We have filtered all the edible seaweed compounds (n = 194) by applying Lipinski's rule of five, Ghosh filter, Veber filter, and PAINS filter (Supplementary File S1). From the total edible seaweed compounds (n = 61), molecules are passed (Supplementary File S2) through the filter without any violations. These compounds showed drug-likeness properties. Christopher A. Lipinski created the Lipinski's rule of five (Ro5), a rule of thumb for determining drug-likeness or to determine if a compound with specific pharmacological activity has chemical and physical properties that would make it a likely orally active drug in humans (34, 66). Pan assay interference compounds (PAINS) are physicochemical filters that can predict the pharmacokinetic features of drugs. The PAINS filter assists in determining if a compound is a biological assay response or not (37). Afterward, to screen the HDAC 2 protein inhibitor ligands, we have performed HTVS molecular docking by using the Schrödinger Glide program.

Glide SP, XP, and QPLD Study

The result of virtual screening revealed the identification of one compound from the edible seaweed as a potential inhibitor as it exhibited non-covalent interactions with the amino acid residues from the catalytic residue of the HDAC 2 protein. The edible seaweed compound dataset containing 61 compounds was initially screened using SwissADME. The high throughput virtual screening module then filtered the 61 edible seaweed compounds. The lower energy score was utilized to rank the top compounds from the dataset. The lower the binding energy, the higher the binding efficiency and, as a result, the stronger the inhibition. The glide scores of seven compounds showing high binding affinity are given in Supplementary Table S3. The XP Gscore was found to be higher in cinnamyl dihydrocinnamate (CMNPD5348) with an energy of −11.829 kcal/mol. Finally, the top hit compound, cinnamyl dihydrocinnamate, is again docked into the binding cavity of HDAC 2 receptor using QPLD. The top hit compound was then subjected to re-dock with HDAC 2 using QPLD for evaluating relative binding interactions and the strength of binding by accurate charge calculation through hybrid quantum mechanics and molecular mechanics method (QM/MM), respectively. The selected compound showed the highest binding affinity toward HDAC 2 with a QPLD DG score of −11.190 kcal/mol (Figures 2A,B). The compound also interacts with the cofactor Zn²⁺ of HDAC 2 protein. Overall, the selected hit compound, cinnamyl dihydrocinnamate's, position midway between the active site channel and the cavity, as well as its direct binding with the catalytic Zn²⁺ ion is optimal for inhibiting HDAC 2 protein, which will prevent deacetylation by interacting with the acetylated lysine residues as a substrate (67). The amino acid residues of the binding cavity within 4 Å that interacts with the selected hit compound were Pro 34, Met 35, Pro 37, Arg 39, Ile 40, Asp 104, Ala 141, Gly 142, Gly 143, Leu 144, His 145, His 146, Tyr 209, Phe 210, Leu 276, Gly 305, Gly 306, Gly 307, Tyr 308, and Asp 401. The hydrophobic zones in the neighborhood of the active site of HDAC 2 is composed of amino acid residues, namely Pro 34, Met 35, Pro 37, Ile 40, Asp 104, Ala 141, Leu 144, Tyr 209, Phe 210, Leu 276, and Tyr 308. Additionally, there were positively charged amino acid residues, namely Arg 39, His 145, and His 146. Binding mode analysis revealed that the compound forms one hydrogen bond interaction with the amino acid residue Tyr 308.

Molecular Dynamics (MD) Simulation Study

Molecular dynamics (MD) simulation studies were carried out to determine the stability of cinnamyl dihydrocinnamate bound HDAC 2 protein complex. A simulation of 100 ns displayed stable conformation while studying the root mean square deviation (RMSD) values. The Cα atoms of HDAC 2 protein were used to analyze the RMSD, which was then plotted against the simulation time, as presented in Figure 3A. The cinnamyl dihydrocinnamate when bound to HDAC 2 showed RMSD (Cα atom of HDAC 2) values between 0.7 and 1.7 Å with an average value of 1.39 Å. In the case of ligand, following the initial equilibration fluctuation, the RMSD for cinnamyl dihydrocinnamate stayed in the range of 1–2.6 Å till the end of the simulation. The overall RMSD of the lead compound cinnamyl dihydrocinnamate is within the acceptable range during the 100-ns simulation, confirming the protein–ligand complex stability. Each amino acid residue plays a critical role in the stability of the protein–ligand complex in the MD simulation study. The RMSF is used to evaluate the flexibility of each amino acid residue and how much it moves or changes throughout a simulation period. The RMSF was employed to investigate the fluctuation of the protein–ligand docking complex in the function of time. The RMSF value was assessed from the MD simulation trajectory and is shown in Figure 3B. If the atoms in the active site and the main chain fluctuated minimally, it indicated that the conformational change was minimal, implying that the reported lead compound was firmly bound within the cavity of the target protein binding pocket (61). The RMSF results revealed that the Cα atoms of HDAC 2 protein bound to cinnamyl dihydrocinnamate have a mean RMSF value of 0.077 Å, indicating fewer fluctuations in the complex structure. The RMSF plot of the cinnamyl dihydrocinnamate-HDAC 2 complex revealed that the lead compound cinnamyl dihydrocinnamate made contact with 13 amino acids of the HDAC 2 protein, namely, Tyr 29 (1.1 Å), Met 35 (0.721 Å), Leu 144 (0.355 Å), His 145 (0.34 Å), His 146 (0.353 Å), Gly 154 (0.567 Å), Phe 155 (0.555 Å), Asp 181 (0.361 Å), His 183 (0.341 Å), Phe 210 (0.545 Å), Asp 269 (0.427 Å), Leu 276 (1.487 Å), and Tyr 308 (0.674 Å). The RMSF values of most residues are <2 Å with the exception of the loop regions, C and N terminal, showing that the residue conformation is relatively stable during the simulation. The RMSF plot above undoubtedly shows that the HDAC 2 protein residues bound to cinnamyl dihydrocinnamate remained stable throughout the simulation. From the results, it can be observed that cinnamyl dihydrocinnamate shows the least fluctuation throughout the 100-ns simulation time. From the 2d ligand-interaction diagram, it was observed that the hydroxyl group of the lead compound cinnamyl dihydrocinnamate interacted with three residues of the binding cavity of HDAC 2 protein via Zn²⁺ interaction with the amino acid residues: Asp 181 (100%), His 183 (100%), and Asp 269 (100%; Figure 3C). The amino acid residue His 145 (98%) also interacted with the hydroxyl group of the lead compound (Figure 3C). The binding characteristics of HDAC 2 with the ligand cinnamyl dihydrocinnamate were analyzed by studying the hydrogen bond interactions. The number of intermolecular hydrogen bonds in MD simulation trajectories of cinnamyl dihydrocinnamate bound HDAC 2 complex is shown in Figure 3D. The HDAC 2 complex with cinnamyl dihydrocinnamate forms an average of one hydrogen bond during most of the simulation time. This result confirmed the strong inhibition of HDAC 2 by the compound cinnamyl dihydrocinnamate during MD simulation. The study showed similar interpretations with the molecular docking results of one hydrogen bond within the active site residues in HDAC 2 protein. Hydrogen bonding is significant in drug design research since it influences drug metabolism, adsorption, and specificity (68). In addition, hydrophobic interactions, ionic or polar interactions, and water bridge-hydrogen bonded protein–ligand interactions mediated by a water molecule contribute to the docked complex's stability during simulation. As a result, the intermolecular contact generated between the HDAC 2 protein and the docked compound was retrieved from the respective MD trajectories using the Desmond module's default parameters. Figure 4A demonstrates the protein–ligand contacts. The dominant contribution to the interaction is given by the Asp 181, His 183, and Asp 269 amino acid residue, where ionic interactions form a 2.0 and 1.0 interaction fraction. The other interactions are given by other residues that help the protein–ligand interaction namely, Met 35, Leu 144, His 145, His 146, Gly 154, Phe 155, and Phe 210, and Tyr 308, where the combination of H-bonds, hydrophobic, and water bridges with interaction fractions are ~0.1–1.0. The protein–ligand's total contacts timeline diagram is prepared to study the intermolecular interactions between HDAC 2 residues and the cinnamyl dihydrocinnamate compound (Figure 4B) after performing 100 ns of MD simulation. The top panel displayed a maximum of 8 specific contacts formed between the protein and the ligand throughout simulations. The bottom panel shows that the residues Asp 269, Phe 210, His 183, Asp 181, Phe 155, His 145, Leu 144, and Met 35 showed crucial interactions with the ligand cinnamyl dihydrocinnamate, which is represented by a darker shade. The amino acid residue Met 35, Leu 145, and Phe 210 also showed interaction during our previous molecular docking findings, within the active site residues of the HDAC 2 pocket. To better comprehend the conformation strains, the torsional dynamics for the rotatable bonds present in the lead compound, cinnamyl dihydrocinnamate were studied and depicted in Figure 5. The plot depicted the lead compound simulation trajectory from 0 to 100 ns, as well as the conformational evolution of rotatable bonds (RB). The bonds were color-coded and rotatable bond torsions were illustrated in a two-dimensional (2D) plot. Throughout the simulation, dial charts indicated the torsion's conformation. The time evolution of the simulation is plotted radially outward, with the initiation of the simulation at the center of the radial map. The torsions' probability density is depicted in bar plots, which depict the information received from the dial plots. The kcal/mol numbers on the left Y-axis of the graphs are in kcal/mol. Cinnamyl dihydrocinnamate has nine rotatable bonds in total, and the dial plots demonstrated that the rotatable bonds are rotated between 180° and 180° (Figure 5). This demonstrated that the ligand cinnamyl dihydrocinnamate had greater flexibility in binding and maintaining a stable conformation in the HDAC 2 protein active site residue.

Free Energy Landscape Analysis

The free energy landscape of (FEL) of achieving global minima of Cα backbone atoms of protein concerning RMSD and radius of gyration (Rg) is displayed in Figure 6. Cinnamyl dihydrocinnamate bound to receptor HDAC 2 achieved the global minima (lowest free energy state) at 1.5 Å RMSD and Rg 29 Å (Figure 6). Due to its great stability and optimal conformation in the bound state of cinnamyl dihydrocinnamate, the FEL was designed for understanding the deterministic behavior of the receptor to the lowest energy state possible. Therefore, FEL is the indicator of protein folding to attain minimum energy state, and that is aptly achieved due to the bound state of ligand cinnamyl dihydrocinnamate. Further free energy landscape (FEL) of cinnamyl dihydrocinnamate-bound receptor complexes exhibited a deep basin over areas of increased free energy with the deep blue color locations (Figure 6) representing the local energy minima and actively promoting the stable conformations, similarly suggested by Singh et al. (69).

PCA and DCCM of the MD Simulation Trajectories

Principal component analysis of the MD simulation trajectories for cinnamyl dihydrocinnamate bound to HDAC 2 was analyzed to interpret the randomized global motion of the atoms of amino acid residues. This analysis interprets the more flexible scattered trajectories owing to the distortion of the protein structure. The internal coordinates' mobility into three-dimensional space in the spatial time of 100 ns was recorded in a covariance matrix, and rational motion of each trajectory is interpreted in the form of orthogonal sets or eigenvectors. The Cα atoms of HDAC 2 from molecular dynamics simulation bound to cinnamyl dihydrocinnamate displayed unordered clustering in PC1 and PC2 modes (Figure 7A). Following this PC2 and PC3 modes displayed better order eigenvalues for the trajectories (Figure 7B). The eigenvectors displayed a less scattered but more positive correlation cluster of the last 50 trajectories (Figure 7C). Most of the trajectories finally settled in PC9 and PC10 (Figure 7D), where the global motion was centered toward the origin of the plot. The global movement of the MD trajectories toward positive correlation eigenvector indicated the compact and converged structure of HDAC 2 complex with cinnamyl dihydrocinnamate. The high periodic global motion was observed along the PC9 and PC10 planes (Figure 7D) due to the grouping of trajectories in a single cluster at the center of the PCA plot. Centering of the trajectories in a single cluster indicates the periodic motion of MD trajectories due to stable conformational global motion. Ordered orientation signifies that the ordered global motion of the trajectories arose due to a stable converged structure. DCCM was generated to analyze the correlative motion of structural domains to attain a stable conformation of the compound after HDAC 2 binding from MD trajectories (Figures 8A,B). The matrix plot's different color codings represent different levels of correlation. The blue color represents a high level of positive correlation, while the red color suggests less or negative correlation among structural domains. The correlation scores on the central mean line (blue) displayed two distinct blocks in cinnamyl dihydrocinnamate (Figure 8A) bound HDAC 2, having a high correlation of movement of the amino acids in HDAC 2 (Figure 8B). In HDAC 2, bounded with the ligand cinnamyl dihydrocinnamate, α helix, β pleated sheet, and loop amino acid residues (100–240) displayed high positive correlation (blue) in three α-helix (α1, α2, α3), one β1 sheet (β1), and five loops (L1, L2, L3, L4, and L5) (Figure 8B). Thus, the DCCM analysis confirmed the previously obtained RMSF result. By studying the dynamics cross-correlation matrix, a similar type of work was earlier reported (59).

MM-GBSA Calculations

To assess the binding free energy of ligands with HDAC 2 protein molecules, the MM-GBSA technique is commonly used (30, 58). The binding free energy of HDAC 2–cinnamyl dihydrocinnamate complex, as well as the impact of other non-bonded interactions energies were calculated (Table 2). With HDAC 2, the ligand cinnamyl dihydrocinnamate has a binding energy of −11.46 ± 4.91 kcal/mol. The non-bonded interactions like ΔG_bindLipo, ΔG_bindCovalent, ΔG_bindHbond, ΔG_bindCoulomb, ΔG_bindSolvGB, and ΔG_bindvdW govern the binding free energy of ΔG_bind. Across all types of interactions, ΔG_bindvdW, ΔG_bindLipo, ΔG_bindCoulomb, and ΔG_bindHbond contributed the most to the average binding energy (Table 2) in the protein–ligand interaction. But ΔG_bindSolvGB and ΔG_bindCovalent have exhibited unfavorable energy contributions and so opposed binding. Thus, the MM-GBSA calculations well-validated the binding energy estimates obtained from molecular docking.

ADME Property Analysis

The QikProp tool was used to evaluate the ADME properties of the compound cinnamyl dihydrocinnamate in this study. In this work, the physicochemical properties and biological roles are discussed. The ADME analysis is a crucial method for determining the efficacy of a potential pharmaceutical compound. The physicochemical and biological properties investigated included molecular weight, # Stars, SASA, FOSA, FISA, PISA, WPSA, Dipole, Donor H-bond, Acceptor H-bond, QPlogPo/w, QPlogS, QPlogkhsa, QplogBB, QPPCaco, percent human oral absorption, rule of five, and rule of three effects. The values for these descriptors are shown in Table 3. The expected ADME values for the compound cinnamyl dihydrocinnamate are all within the approved range, indicating better drug likeliness and pharmacokinetic characteristics.

Toxicity Prediction of Cinnamyl Dihydrocinnamate

The predicted values related to the AMES toxicity and hepatotoxicity revealed that the cinnamyl dihydrocinnamate was nontoxic and non-mutagenic (Table 3). The liver is associated with energy exchanges as well as xenobiotic and drug biotransformation. Damage to the liver leads to disruption in normal metabolism and can possibly lead to liver failure (70). The hepatotoxicity description predicted that the compounds would not be toxic to the liver. Inhibition of potassium channels encoded by hERG-I and II resulted in fatal ventricular arrhythmia. Inhibition of hERG-I and II resulted in the withdrawal of many compounds from the pharmaceutical industry (71). The selected compound has shown no inhibition of hERG-I and II. However, the selected compound predicted skin sensitization. The compound may induce an allergic reaction when applied dermally. The toxicity analysis server predicted LD50, the lowest-observed-adverse-effect level (LOAEL) and the maximum tolerated dosage range of the selected compound, and the predicted scores are presented in Table 4. It is worth noting that the oral toxicity of the compound evaluated by this tool corresponds to the toxicity risk observed in animal models. The maximum recommended human tolerated dose (MRTD) (log mg/kg/day) bestowed the toxic dose threshold of chemicals in humans in phase I clinical trials. It is considered low when ≤ 0.477 log mg/kg/day and high >0.477 log mg/kg/day. Our findings showed that the compound has high MRTD.

Density functional theory. The electronic structures of compounds are linked to their pharmacological activities. The energy and distribution patterns of their orbitals, specifically the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO), which are commonly referred to as frontier molecular orbitals (FMOs), are also crucial to explore (72). Frontier molecular orbitals were used to analyze the compounds' electrical characteristics and their electron transport potential. The biological activity and molecular characteristics of the substances are determined by these orbitals. To reveal insight on the electronic structural properties of our selected compound, cinnamyl dihydrocinnamate, we have optimized the compound structure utilizing the density functional theory (DFT) calculation using B3LYP/6-31G^** basic set (single point energy calculation). The HOMO is the highest-energy (electron-rich) orbital that has the potential to provide electrons, whereas the LUMO is the lowest-lying empty orbital (lack of electron) that has the capacity to accept electrons (73). Figures 9A,B illustrate the energy of HOMO and LUMO, their energy gap (E), and the molecular electrostatic potential surface (MESP) of the studied compound. The LUMO and HOMO energies influence how a compound interacts with other species, as well as its chemical reactivity and kinetic stability. Low-band energy gaps (ΔE), which signify stronger compound reactivity with the receptor, were indicated by the smallest difference between LUMO and HOMO. Since the LUMO and HOMO are responsible for the charges exchanged in a chemical reaction, a lower band energy gap (ΔE) indicates more reactivity of the compound. The decrease in energy gap values corresponds to a rise in the order of reactivity (60, 74–76). A Lower energy gap was observed with the selected compound, cinnamyl dihydrocinnamate (CMNPD5348) (ΔE = – 0.320118 eV) (Figure 9A). In cinnamyl dihydrocinnamate, HOMO localizes on electronegative atoms, and highly delocalized HOMO suggests that electrons can move around the molecule more freely, resulting in better intramolecular charge transfer. The presence of negative HOMO–LUMO values in our chosen compound indicates high stability, which is required to establish a stable protein–ligand interaction. More energy gap led to a weak affinity of the protein–ligand interaction. Furthermore, the energy gap is a good method for identifying the most active molecules. The intermolecular charge transfer and bioactivity of a compound are significantly influenced by the energy gap between the HOMO and LUMO energies (61, 74). The molecular electrostatic potential (MEP) plot can be used to gain insight into a compound's three-dimensional structural and topological properties. It identifies the reactive sites for electrophilic and nucleophilic attack in a compound for binding to protein in protein–ligand interaction, as well as the charge distribution (positive and negative) in the selected compound (72). The various charges are represented by different colors. The red color denotes negative potential (negative regions), which has a propensity for attracting protons; the blue color denotes positive potential (positive regions), which has an attraction for repelling protons; and the green color denotes zero potential. Potential increases were denoted as red<orange<yellow<green<blue. To make a hydrogen bond with the protein, the negative region is crucial (60, 74, 77). In our selected compound, the blue color (positive region) is situated around all surfaces of the compound, which indicates their potentiality to act as a H-bond donor in protein–ligand interactions, whereas highly negative regions were situated at the oxygen atom. This region may be responsible for the protein's hydrogen-bonding interactions. These results corroborate our molecular docking finding. In our selected compound, the blue color (positive region) indicates their potentiality to act as the H-bond donor in protein–ligand interactions. These results corroborate our molecular docking finding.

Conclusions

Overexpression of HDAC 2, a class I HDAC isoform, has been discovered in a range of cancers. In a growing number of research articles, data, and proof, the role of HDAC 2 in cancer development and progression has been established. The primary goal of designing HDAC specific inhibitors is to avoid non-selective interactions with other HDAC classes, which can result in substantial side effects. In this study, we have employed a computational-based drug design workflow to recognize potent HDAC 2 inhibitors. To investigate potential HDAC 2 inhibitors as lead compounds from edible seaweed compounds, we used a combination of in silico methodologies such as drug-likeness, molecular docking, molecular dynamics simulation, MM-GBSA, PCA, DCCM, free energy landscape, and density functional theory based virtual screening. The lead compound cinnamyl dihydrocinnamate (CMNPD5348) interacted with HDAC 2 protein efficiently according to our findings. Molecular dynamics simulations and binding free energy calculations revealed that the docking complex has a high level of binding stability and a low binding free energy. At 100 ns of molecular dynamics simulation time, the metal ion Zn and ligand interaction was also observed, which is crucial for inhibiting the HDAC 2 protein's enzymatic activity. The hit compound also showed high pharmacokinetic characteristics in the drug-likeness and ADME analysis, indicating that it could be a promising lead candidate. Cinnamyl dihydrocinnamate has been identified in edible seaweed Caulerpa racemosa, which is mainly found in shallow water locations around the world. In this study, the compounds were discovered to be effective HDAC 2 inhibitors. This discovery is expected to lead to the development of efficient HDAC 2 inhibitors for cancer therapy. The fishing industry generates a lot of seaweed debris, as well as other substances. This study will motivate researchers to look into and extract beneficial compounds from edible seaweed waste for biomedical research and therapeutic applications. This study will also provide an important step in the exploration of bioactive compounds from seaweed waste. In addition, the outcomes of in-vitro and in-vivo analysis will be necessary to validate the molecular docking, molecular dynamics simulation, and DFT investigations.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding authors.

Author Contributions

KB, IA, SP, BR, TS, and DB: conceived and designed the experiments. KB, IA, SP, BR, TS, HP, AG, DB, and MA: writing—original draft preparation. KB, IA, SP, HE, WW, and MA: formatting and editing according journal guidelines. KB, IA, BR, TS, HP, AG, DB, ZA, HE, and WW: writing—review and editing. All authors have read and agreed to the published version of the manuscript.

Conflict of Interest

SP is employed by NatNov Bioscience Private Limited.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

KB is thankful to MHRD, Government of India and Department of Bioengineering and Technology, Gauhati University for supporting the TEQIP-III fellowship. The authors would like to acknowledge their respective institutes/universities for research and financial support.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnut.2022.889276/full#supplementary-material

References