In Silico Computational Studies of Bioactive Secondary Metabolites from Wedelia trilobata against Anti-Apoptotic B-Cell Lymphoma-2 (Bcl-2) Protein Associated with Cancer Cell Survival and Resistance

In the present study, the binding affinity of 52 bioactive secondary metabolites from Wedelia trilobata towards the anti-apoptotic B-cell lymphoma-2 (Bcl-2) protein (PDB: 2W3L) structure was identified by using in silico molecular docking and molecular dynamics simulation. The molecular docking results demonstrated that the binding energies of docked compounds with Bcl-2 protein ranged from −5.3 kcal/mol to −10.1 kcal/mol. However, the lowest binding energy (−10.1 kcal/mol) was offered by Friedelin against Bcl-2 protein when compared to other metabolites and the standard drug Obatoclax (−8.4 kcal/mol). The molecular dynamics simulations revealed that the Friedelin-Bcl-2 protein complex was found to be stable throughout the simulation period of 100 ns. Overall, the predicted Absorption, Distribution, Metabolism, Excretion, and Toxicity (ADMET) properties of Friedelin are relatively better than Obatoclax, with the most noticeable differences in many parameters where Friedelin has no AMES toxicity, hepatotoxicity, and skin sensitization. The ADMET profiling of selected compounds supported their in silico drug-likeness properties. Based on the computational analyses, the present study concluded that Friedelin of W. trilobata was found to be the potential inhibitor of the Bcl-2 protein, which merits attention for further in vitro and in vivo studies before clinical trials.


Introduction
The B cell lymphoma 2 (Bcl-2) family proteins that appropriately regulate apoptotic cell death are associated with the progression and development of several cancers [1,2]. The first founding member of the Bcl-2 family of apoptosis regulators was the Bcl-2 protein, which protects the cells from apoptosis by directly binding to the pro-apoptotic Bcl-2 family members, thereby playing a crucial role in cancer cell survival and resistance in response to chemotherapy or radiation therapy [3]. Most prescribed chemotherapy drugs are synthetically prepared and have demonstrated their cytotoxic effects, not only on cancer cells, but also on healthy cells. However, recent findings suggest that naturally derived secondary metabolites from plants are well known to exhibit selective cytotoxicity toward cancer cells and minimal toxicity towards normal cells [4]. Intense research into finding out the new cancer treatment targets has led to the identification of Bcl-2 overexpression as a hallmark for cancer and Bcl-2 inhibition as a possible cancer treatment approach [5].
The Bcl-2 protein inhibition is a promising therapeutic target for developing new anticancer drugs because it is a key regulator of the intrinsic apoptosis pathway involved in various cancers [3,6]. The Bcl-2 inhibitors offer a validated pre-clinical platform for testing new therapeutic drugs and have been used to treat multiple cancers. Therefore, this protein can be a potential target receptor for developing new therapeutic agents to treat various cancers. The literature reveals that numerous Bcl-2 targeting experimental drugs (such as Obatoclax, Gossypol, Oblimersen, Epigallocatechin-3-gallate, and Navitoclax) (Supplementary Figure 1) are still now undergoing clinical studies [2]. Hence, it is of further interest to identify the new Bcl-2 inhibitors from the plant-derived compounds. The molecular docking study is widely used as a powerful tool in the pharmaceutical industry, particularly in the analysis of protein structure activity relationships [7][8][9]. It plays a vital role in predicting the molecular interactions of small ligand molecules with appropriate binding sites of the targeted protein. Moreover, molecular dynamics simulations can be used to investigate the conformational changes that occur during protein-ligand interactions. In addition to the intrinsic pathway, the extrinsic pathway is activated by tumor necrosis factor (TNF), which triggers apoptosis by binding to its death receptors. The selective TNF receptors can presently be identified as one of the novel drug targets for the treatment of cancer [10].
Wedelia trilobata, also known as Sphagneticola trilobata (L.), is a traditionally used medicinal plant renowned for its therapeutic properties in the treatment of various diseases [11]. The plant is a rich source of natural medicines that can be used as hepatoprotective, anti-inflammatory, antimicrobial, antioxidant, anti-malarial, anticancer, and antidiabetic agents [12][13][14]. In the present study, the bioactive secondary metabolites from W. trilobata were explored to identify their inhibitory potential against the anti-apoptotic Bcl-2 protein associated with cancer cell survival and resistance by employing the in silico computational methods like-molecular docking, molecular dynamics simulation, and ADMET analysis). Ramachandran plot analysis of Bcl-2 protein (PDB: 2W3L) structure generated by the PROCHECK web server. The red color regions in the graph indicate protein residues in the most allowed regions; the brown color regions indicate the residues in additional allowed regions; the yellow color regions indicate the residues in generously allowed regions.

Results and Discussion
The Bcl-2 is an important regulatory protein of apoptosis or programmed cell death that regulates the fundamental biological process by the formation of heterodimers with selective pro-apoptotic Bcl-2 family members [1]. Preventing the targeted anti-apoptotic Bcl-2 protein in the present study is the key to stimulating the apoptosis process in cancer cells. In the present study, in silico computational methods can be exploited to build the appropriate drug targets as well as to design novel therapeutic agents to block the apoptosis regulator Bcl-2. However, the ideal criteria for each drug in the docking techniques are that it should have the highest binding affinity for inhibiting the Bcl-2 protein. The obtained PROCHECK plot showed that most of the amino acid residues of Bcl-2 protein are found within their most favored regions ( Figure 1). The Z-score generated was found within the range of its native protein, which validates that the protein model used is high quality (Figure 2). The binding sites of the structure of chimaeric Bcl2-xL (PDB: 2W3L) with phenyl tetrahydroisoquinoline amide complex is given as Supplementary Figure S2.

Molecular Docking Study
The docked results showed that among the bioactive secondary metabolites of W. trilobata, Friedelin (−10.1 kcal/mol) has the best binding affinity toward the target protein, followed by Friedelinol (−9.6 kcal/mol) in comparison with the remaining metabolites and the standard drug Obatoclax (−8.4 kcal/mol). The chemical structures of the top potential metabolites from W. trilobata and the standard drug are represented in Figure 3. The molecular docking analysis data of the bioactive secondary metabolites from W. trilobata against Bcl-2 protein are presented in Table 1. The detailed interaction (2D and 3D) between the selected docked compounds and Bcl-2 protein is shown in Figures 4 and 5.   Grandiflorenic acid −7.7 3. ent-kaura-9(11), 16-dien-19 oic acid methyl ester −7.6 4.
Wedeliatrilolactone  Generally, a good inhibitor has a sufficient number of hydrogen bonds and other hydrophobic interactions, which are very important for protein-ligand interactions [15].
During protein-ligand interaction, new hydrogen bonds are formed between them by breaking the existing hydrogen bonds with water molecules. The residues ALA59, PHE63, ALA72, VAL107, and TYR161 of the Bcl-2 protein were found to interact with the Friedelin of W. trilobata ( Figure 5). Therefore, the optimal binding features of the Friedelin of W. trilobata and Obatoclax with the Bcl-2 protein are supported for further study and validation through molecular dynamics simulation. Adewole and Ishola [16] have reported that six compounds (Campesterol, Cycloartenol, Oleanolic acid, β-Sitosterol, Stigmasterol, and Ursolic acid) identified from 47 compounds in Morinda lucida Benth were found to display higher binding affinities towards Bcl-2 (PDB: 2W3L) than the reference ligand Obatoclax. In addition to hydrogen bond interaction, the hydrophobic interaction was also predominantly involved in the strong binding of potential compounds with binding pockets of Bcl-2. In addition, Mohamad Rosdi et al. [17] have reported that the pre-clinical Bcl-2 inhibitor Obatoclax was used for comparison of the bioactive compounds from Annona muricata with the binding energy of −7.01 kcal/mol, including interacting residues of Phe-109, Phe-101, Met-112, Val-130, Leu-134, Gly-142, Arg-143, Ala-146, Phe-150, and Val-153.

Structural Activity Relationship Analysis
All the selected compounds and Obatoclax were selected against the Bcl-2 protein for further structural similarity and activity relationship analysis based on the correlation of canonical SMILES structure similarity and the binding energy of selected compounds with the help of Data Warrior software. The structurally similar compounds were found to be in a group with a similar range of binding affinity values ( Figure 6).  Table 1.

Molecular Dynamics Simulation Studies
It is well established that molecular dynamics simulations can be used as powerful tools for studying biological processes at the atomic level, such as protein conformational stability, flexibility, and behavior that occurs under different physiological conditions [18,19]. In the current study, the molecular dynamics simulations were executed to determine the structural stability of the docked complex of the Bcl-2 protein with Friedelin and Obatoclax within the time frame of 100 ns. The 2W3L complexed with Friedelin was selected based on its lower binding energy, while the 2W3L complexed with the standard drug Obatoclax was also considered in the present study to obtain a better comparison. Thus, both the selected complexes were subjected to the molecular dynamics simulations at 100 ns.
The root mean square deviation (RMSD), which is calculated for the atoms in the complexes' backbones, is an important parameter for analyzing the equilibration of molecular dynamics trajectories [20]. The measurement of the RMSD backbone for two complexes offered information regarding the conformational stability. The 2W3L-Friedelin complex maintained a stable RMSD profile throughout the rest of simulation periods, and with no over-fluctuation, demonstrating the rigid conformation. However, the 2W3L-Obatoclax complex showed a similar trend of RMSD profile until 15 ns, and thereafter it increased ( Figure 7A). The root mean square fluctuation (RMSF) analysis was carried out to estimate the average fluctuations of protein residues during molecular dynamics simulation [21]. The RMSF plots of the 2W3L-Friedelin complex were predicted with fluctuations only at the loop regions and terminal ends compared to the 2W3L-Obatoclax complex, which showed more fluctuations, thus indicating the stable interactions between 2W3L-Friedelin complex ( Figure 7B).
The radius of gyration (Rg) considers the varied masses when calculating the root mean square distances with the rotation of the central axis. As the simulation progresses, it considers the protein folding, shape, and capability of the entire trajectory at each time step. The Rg plots of the protein alone and protein-ligand complex, also equilibrated within 1.85-2.0 nm, indicating the complex stability ( Figure 7C). The region surrounding the hydrophobic core generated by the protein-ligand interaction is measured by the solvent accessible surface area (SASA). The SASA plots of protein alone and protein-ligand complex were calculated with the value of 215-235 nm 2 , indicating the consistency in the formation of the complexes ( Figure 7D). The intermolecular hydrogen bonds are ubiquitous and play a vital role in protein folding and protein-ligand interactions [15]. The stability of the hydrogen bond network formed in the protein-ligand complex was calculated throughout the simulation period of 100 ns. The total number of hydrogen bonds in complexes versus time at 300 K can be seen in Figure 7E. The 2W3L-Friedelin complex exhibited the required amount of hydrogen bonds throughout the simulation, indicating that Friedelin has stable and strong hydrogen bonds with the 2W3L protein.

Binding Free Energy Calculation
During molecular dynamics simulations, a variety of energy metrics, including electrostatic, van der Waal's, SASA, polar solvation, and binding energy, were used to estimate the strength of the interaction between the ligand and target protein. In this study, the van der Waal energy, electrostatic energy, SASA energy, and binding energy were mainly used to form the protein-ligand complex. As the numbers were positive, the polar solvation energy was predicted to have no effect on the formation of the protein-ligand complex. The 2W3L complexed with Friedelin, along with the standard drug Obatoclax, were considered for their binding energy calculation. The standard deviation of the protein-ligand complex was also computed. A higher standard deviation indicates that the data values are dispersed over a larger range, whereas a lower standard deviation indicates that the data values are closer to the mean (or expected value). However, there was no higher standard deviation in the Friedelin-receptor complex compared to Obatoclax, which has high standard deviation values, thereby indicating that Friedelin strongly binds to the 2W3L protein with higher binding affinity and stable interaction. The binding free energy calculation of the 2W3L target protein complexed with Friedelin and Obatoclax are represented in Table 2.

Prediction of ADMET Properties of Ligand Molecules
The main issue associated with the further development of novel chemical compounds into therapeutic drugs is their poor pharmacokinetics, including their physicochemical properties and ADMET properties [22]. The ADMET properties are phenomena that are directly connected to how a chemical substance behaves inside the human body. Each of the properties of ADMET will depict what happens when a chemical substance interacts with different organs in the body. The predicted ADMET properties of compounds are crucial, especially for foreign chemical compounds consumed in high doses or over extended periods. Information regarding the ADMET properties of chemical compounds is primarily required to develop novel drug compounds. Computational model-based predictions serve as virtual filters for drug-like features and help medicinal chemists to design drug candidates rationally [22,23].
The molecular polar surface area (PSA) is a valuable characteristic that is strongly related to drug absorption [24]. The PSA of Friedelin and Obatoclax were >140 Å 2 , indicating that they showed a strong polarity and were not easily absorbed by the body. Lipophilicity is the most important physicochemical characteristic of substances, connected to their solubility and human intestinal membrane permeability [25]. It was noticed that Friedelin was considered as poor lipophilicity as it displayed logP ≥5 compared to Obatoclax. The result suggested that Friedelin has poor absorption and permeation abilities. Friedelin was predicted to show moderately high Caco-2 permeability and was found to absorb more easily via the human intestine than Obatoclax. Compared to Obatoclax, Friedelin was considered to show relatively low skin permeability. The results suggested that Friedelin was not a P-glycoprotein substrate but was found to be a P-glycoprotein inhibitor when compared to Obatoclax.
The predicted results also showed that the volume of distribution at steady state (VDss) and fraction unbound (Fu) in Friedelin's plasma was lower than Obatoclax. The Friedelin, showing log BB of 0.72, was considered to cross the blood-brain barrier easily compared to Obatoclax and could penetrate the central nervous system. Friedelin was the only substrate for CYP3A4 and was not an inhibitor of Cytochromes P450 (CYPs). The total clearance of Friedelin was found to be the lowest and not a renal OCT2 substrate. Friedelin may not be toxic in the AMES test, but proved to have the highest oral rat acute toxicity with LD 50 of 2.64 mol/kg. Friedelin may not be hepatotoxic compared to Obatoclax, but it did not have any skin sensitizing potency. Friedelin showed the lowest Tetrahymena pyriformis toxicity (0.3 µg/L) and minnow toxicity (−2.384 mM) compared to Obatoclax. Therefore, the predicted results indicated that Friedelin's physicochemical and ADMET properties supported the in silico computational studies compared to the drug Obatoclax (Supplementary Table S1).

Preparation of Ligands
For the present study, the ligand molecules were selected based on the three-dimensional (3D) chemical structures of 52 bioactive secondary metabolites from W. trilobata [26]. In addition, the anti-apoptotic known inhibitor Obatoclax was used as the standard drug for the present study. All of the selected phytochemicals' chemical structures were retrieved in structured data format (SDF) from the PubChem compound database [27], and those that were not available in the database were drawn using Chemaxon's chemical drawing tool (MarvinSketch version 18.30). OpenBabel version 2.3.2 was utilized to prepare the ligands, which were then exported in protein data bank (PDB) format [23]. The geometries of ligand PDB files were optimized using the PRODRG server before molecular docking.

Preparation of Bcl-2 Protein Structure
Based on the literature, in the present study, the crystal structure of chimeric Bcl2-xL and phenyl tetrahydroisoquinoline amide complex with the PDB ID: 2W3L [16] was selected as the main therapeutic target protein. The 3D X-Ray crystallographic structure of the target protein was retrieved in PDB format from Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB PDB) database [28]. For further docking studies, the protein structure was processed by removing the co-crystallized ligand and water molecules bound to it and subsequent addition of hydrogens and Kollmann charges [29]. Swiss-PDB Viewer software (version 4.10) was used to perform the energy minimization of the protein structure, with empirical force fields to generate its lower energy conformation, which suggests a more stable conformation. While modeling the protein structure, the procedure essentially optimized the conformational errors in the structure's geometry. In addition, the geometry optimization method used the steepest descent algorithm with the GROMOS96 force field [30].

Validation of Protein Structure
The model structure of the Bcl-2 protein was further validated in the PROCHECK online server using the Ramachandran plot, which inferred the reliability of the protein structure model developed in the SWISS MODEL server. Each amino acid included in the protein may be visualized in the plot, along with its highly favored, preferred, and disallowed phi (ϕ) and psi (ψ) dihedral angles in degrees [18]. The Protein Structure Analysis (ProSA) web tool was also used to check the protein structure and analyze the protein model quality. In addition, the selected protein was validated as high-quality when its Z-score fell within the range of its corresponding native protein.

Molecular Docking Analysis
The molecular docking experiment was executed using AutoDock Vina software integrated in the PyRx 0.8 tool, based on the Lamarckian Genetic Algorithm, to study the ligand interaction with the active site residues of Bcl-2 structure [31,32]. The cubical grid box dimension for the XYZ coordinates was fixed at 60 × 60 × 60 Å with 0.375 Å around the active binding sites of the 2W3L protein. The whole protein was surrounded within the dimensions of a grid box with docking exhaustiveness of 100 poses. The docking parameters were set to default, and the protein coordinates were saved in pdbqt file format. The docking results were represented in binding affinities expressed in kilocalories per mole (kcal/mol). The best docking pose was the conformation with the lowest binding energy. After docking, the best binding poses of each ligand molecule were visualized, and their interactions with the Bcl-2 protein structure were analyzed using BIOVIA Discovery Studio Visualizer [33]. The re-docking (self-docking) was performed with the best compounds and Bcl-2 protein to validate the accuracy of the docking methodology.

Structural Activity Relationship Analysis
The Data Warrior software (Version 5.5.0) was used to study the structural activity relationship analysis of the compounds used in the study, with the Bcl-2 protein based on their canonical SMILES and binding energy [34][35][36].

Molecular Dynamics Simulation
The 2W3L-Friedelin complex, which showed the lowest binding energy, was further subjected to the molecular dynamics simulation compared with the 2W3L-Obatoclax complex as a reference protein-ligand complex. The molecular dynamics simulations were executed with GROningen MAchine for Chemical Simulations (GROMACS) version 2018.1 using a CHARMM36 force field to confirm the conformational stability of identified ligand molecule with the Bcl-2 protein [37,38]. The ligand topology files were generated with the help of the CGenFF server [39]. Initially, the pdb2gmx module of the GROMACS was utilized to add hydrogens to the system. The protein-ligand complex was then solvated by placing the complex in a pre-equilibrated rhombic dodecahedron (RHDO) solvent box with a minimum distance of 10 Å from the box boundaries. The dodecahedral simulation box model was filled up with the TIP3P water molecules. The solvated system was then neutralized by adding an adequate number of counter ions (Na + or Cl − ) and co-ions to maintain the complex system's electroneutrality states. The steepest descent method was employed for the energy minimization of systems to 1000 steps to remove the short contacts and atoms overlaps. The system was equilibrated for 100 ps at 300 K under constant pressure conditions in the canonical ensemble (NVT) and isothermal-isobaric (NPT) ensemble. Finally, the molecular dynamics simulations were executed at 100 ns, with an integration time of 0.002 ps under the same conditions [40].

Binding Free Energy Calculation
The outcome of the molecular dynamics simulation run for the target protein complexed with Friedelin and standard Obatoclax was subjected to the calculation of binding free energy using the Molecular Mechanics Poisson-Boltzmann (Generalized-Born) Surface Area (MM-PBSA) approach. The thermodynamic cycle is typically used to calculate the protein-ligand binding free energy based on the molecular dynamics simulations. The gmx_MMPBSA which uses the GROMACS 2018.1 trajectories was used to calculate the binding free energy for each protein-ligand complex [41]. The binding free energy was calculated by the gmx_MMPBSA program using three main components: molecular mechanical energy, polar solvation energy, and apolar solvation energy. The calculations were based on the molecular dynamics trajectories of the last 100 ns, which compute ∆G with dt 1000 frames. The free binding energy calculation using Equations (1) and (2) G Binding , binding free energy; G Complex , total free energy of protein-ligand complex; G Protein and G Ligand , total free energies of protein and ligand in solvent, respectively; ∆G, standard free energy; ∆E MM, average molecular mechanics potential energy in vacuum; ∆G Solvation, solvation energy; ∆E, total energy of bonded as well as non-bonded interactions; ∆S, change in entropy of system upon ligand binding; and T, temperature expressed in Kelvin.

Prediction of ADMET Properties of Ligands
The early prediction of ADMET profiles of chemically synthesized and eco-friendly pharmaceutical drugs has revolutionized disease management methods using the most popular open access in silico tools. The pkCSM database server was used to determine the physicochemical properties and ADMET properties of target ligand molecules, which facilitated identifying promising drug candidate molecules [23,42].

Conclusions
In summary, the molecular docking study suggested that the binding orientation of Friedelin of W. trilobata within the binding sites of the Bcl-2 protein resulted in inhibition of enzyme activity compared to other metabolites and Obatoclax (as a standard drug). The molecular dynamics simulation studies showed that the 2W3L-Friedelin complex maintained stable RMSD and Rg profiles, along with strong hydrogen bonds throughout the simulation period of 100 ns compared to the 2W3L-Obatoclax complex. The pkCSM provided the predicted ADMET properties of Friedelin, which are relatively better than Obatoclax, wherein Friedelin has no AMES toxicity, hepatotoxicity, or skin sensitization. This study concluded that Friedelin of W. trilobata is one of the good inhibitory compounds of the Bcl-2 protein. Therefore, it seems that the improved binding features of Friedelin with the Bcl-2 protein supports further consideration towards in vitro and in vivo validation before clinical trials.

Data Availability Statement:
The data presented in this study are available within the manuscript.