Bulbine frutescens phytochemicals as novel ABC- transporter inhibitor: a molecular docking and molecular dynamics simulation study

Aim: The present in silico study aimed to evaluate the ATP-binding cassette (ABC) transporter inhibition potential of Bulbine frutescens (B. frutescens) phytochemicals. Methods: Several previous studies and databases were used to retrieve the ligands and target protein structure. The molecular docking study was performed using the Auto Dock Tools, and the GROMACS package was applied to accomplish molecular dynamics simulation. Results: Utilizing the molecular docking and simulation approach, ~25 phytochemicals were screened against the ABC transporter protein. Docking score analysis revealed that B. frutescens phytochemical 4’-Demethylknipholone 2’-β-D-glucopyranoside exhibited strong binding on the ABC transporter protein with a minimum binding score -9.8 kcal/mol in comparison to the standard ABC transporter inhibitor diltiazem (-6.86 kcal/mol). Furthermore, molecular dynamics simulation for 4’-Demethylknipholone 2’-β-D-glucopyranoside showed an acceptable root mean square deviation, radius of gyration, root mean square fluctuation, and hydrogen bond, in addition to other lead compounds. Page 2 of 13 Kushwaha et al. J Cancer Metastasis Treat 2021;7:2 I http://dx.doi.org/10.20517/2394-4722.2020.92 Conclusion: The in-silico study demonstrated that B. frutescens phytochemical 4’-Demethylknipholone 2’-β-Dglucopyranoside possesses anti-drug resistance properties and requires further testing in preclinical settings.


INTRODUCTION
Cancer is one of the major health problems worldwide. In addition to surgery, radiation, and hormonal therapy, chemotherapy is an important strategy in the management of cancer [1,2] . Cancer drug resistance is a well-known phenomenon that arises when the disease becomes refractory to chemotherapy [3] . One of the mechanisms related to drug resistance in cancer cells is the involvement of drug efflux transporters such as P-glycoprotein, also known as ABC (ATP-binding cassette) transporters. Besides their normal role, these proteins facilitate the movement of drugs and their metabolites across the cellular membrane. ABC transporters-mediated drug efflux is the cause of the relapse of disease in patients and decreases therapeutic output [4] . Therefore, ABC transporters remains an important target for therapeutic strategy against this devastating disease [4] .
Phytochemicals are known to possess anticancer activity and play an important role in the management of the disease [5] . Various properties such as reducing side-effects, efficacy in chronic conditions, costeffectiveness, and widespread usage of natural compounds give them a role as chemotherapeutic agent. Thus, currently, natural compounds are widely used in cancer chemotherapy and have a promising future [6] . Bulbine frutescens (B. frutescens) is a flowering plant that belongs to the Asphodelaceae family. It is a frequently used medicinal plant found in South Africa for the treatment of burns and skin wounds [7] . The potent anti-proliferative activity of B. frutescens was first reported in 2014 [8] . There are few published articles on B. frutescens species which encompass the isolation of its phytochemicals and their biological activities. Aqueous extracts of B. frutescens can increase glucose utilization activity in vitro [9] . B. frutescens phytochemical isofuranonapthoquinone 1 is a known inhibitor of glutathione transferase P1-1 protein involved in anticancer drug detoxification-mediated drug resistance [10] . The ethanolic and aqueous extracts of the Bulbine genus inhibit HIV-1 protease activity. Recently, we reported that B. frutescens phytochemicals (extract) have strong potential to induce apoptosis and inhibit Notch signaling pathway in breast cancer cells [11] .
The computational approach is a fast-growing stream that employs an integrative approach to characterize the selective target of a particular disease and explore the interactions involved with drugs or other molecules [12] . Computer-based drug discovery is a highly time-and cost-effective strategy. Various online and offline tools and programs are available which have been used in computer-based drug discovery and make the process easier and cost effective. In the present study, we screened the ABC transporter (PDB 3G60) inhibitory potential of B. frutescens phytochemicals and predicted their possible IC 50 using an in silico approach.

Software, servers, and tools
The present study was performed on HP Workstation with 3.40 GHz processor, 16 GB RAM, and 1 TB hard disk running the Ubuntu 16.04.5 LTS and Windows 7 operating systems. Bioinformatics tools, such as AutoDock, PyRx, and GROMACS, were used for molecular docking and simulation studies. Pymol, VMD (visual molecular dynamics), Open Babel, LigPlot + (v.1.4.5), Marvin Sketch 15.9.7.0 UCSF Chimera, and Avogadro were employed for visualization purposes. National Center for Biotechnology Information (NCBI), Protein data bank (PDB), PubChem, Swiss Model, Charmm-gui, RAMPAGE, etc. online servers/ resources were used for the retrieval, evaluation, and analysis of the data.

Identification, retrieval, and preparation of ligands
The phytochemicals present in B. frutescens were searched in the published literature [13] . The 2D/3D structures of identified phytochemicals and standard inhibitor compounds of the test protein were downloaded from the NCBI PubChem database. Open Babel molecule format converter [14] and Marvin Sketch 15.9.7.0 software packages were utilized to perform 2D to 3D conformation and .sdf to .pdb format conversions, respectively. The force field mmff94 and conjugate gradients optimization algorithm were applied to minimize the test ligands energy using PyRx-Python prescription 0.8 for 200 steps [15] .

Retrieval and preparation of the receptor
The 3D structure of ABC-transporter protein or P-glycoprotein (PDB ID 3G60) was retrieved from RCSB-PDB [16] with a resolution of 4.40 Å. The protein structure was loaded into UCSF Chimera for molecular docking preparation [17] . The protein model was cleaned and optimized by removing the ligands and heteroatoms such as water. The protein structure was minimized by using UCSF Chimera software with a 100-step steepest descent and a 10-step conjugate gradient method. The step size of both methods was 0.02 Å.

Molecular docking
Auto Dock Tools 1.5.6 (ADT) was used for the docking of the protein and ligands [18] . Gestgeiger partial charges were allocated and torsions were introduced to the test ligands. Using ADT, Kollman charges, polar hydrogen atoms, and solvation parameters were applied. The Lamarckian Genetic Algorithm was applied to discover the place for active binding. A maximum of 27,000 GA operations were made for each of the thirty independent runs on a single population of 150 individuals. Operator loads for cross overrate, elitism, and gene mutation rate were set as default parameters.

Molecular dynamics simulation
The GROMACS (version 2018) package was used to perform molecular dynamics simulation analysis of the dynamic behavior of membrane-embedded P-glycoprotein complexes with AMBER Lipid14 force field. The receptor topology was constructed by pdb2gmx. Online server CHARMM-gui was used to insert the dipalmitoylphosphatidylcholine (DPPC) membrane. The ligand potential energy function was derived using AMBER GAFF parameters. SPCe water model was applied to represent the water molecules explicitly. To neutralize the positive charge on the receptor, the respective number of negative ions was added. Bad contact inside the system was discarded by energy minimization through the steepest descent method, and the long-range interactions were computed by particle mesh Ewald at every step. The energy minimized system underwent 100 ps (NVT) and 1ns (NPT) equilibration run at constant T (323 K) and pressure (1 bar). We employed modified Berendsen thermostat and Parrinello-Rahman barostat for NVT and NPT equilibrations. In addition, to avoid the hot solvent-cold solute problem, temperature coupling was applied by the indexing system into water and non-water components. On the completion of two equilibration phases, the production run was conducted for 10 ns after taking away the position restraints. Throughout the data collection phase, the pressure and temperature were maintained with the same methods as in the NPT and NVT phases, respectively. Moreover, the P-glycoprotein ligand complex coordinates were saved every 1 ps. Structural analysis results were processed by using modules of GROMACS package. The appropriate GROMACS command was used for the root mean square fluctuation (RMSF), Root Mean Square Deviation (RMSD), radius of gyration, peptide-peptide hydrogen bond, peptide-water hydrogen bond, and SASA calculations. xmgrace software was used to visualize these GROMACS data calculations.

Molecular Mechanics/Generalized Born Model and Solvent Accessibility analysis
Prime Molecular Mechanics/Generalized Born Model and Solvent Accessibility (MM-GBSA) was performed to calculate the protein and protein-ligand binding energies. The analysis utilized the VSGB solvent model, OPLS_2005 force field, and rotamer search algorithms. The glide pose viewer file was utilized to compute the total binding free energy using the Prime MM-GBSA simulation aid. The MM/ GBSA analysis provided relative binding affinity of the target protein and target protein-ligand complex (kcal/mol). High negative MM/GBSA binding energies indicate stronger binding [19] .

Binding mode analysis of ligands
The phytochemicals present in B. frutescens were searched from the published literature and their twodimensional structures are depicted in Figure 1.

Receptor-ligand energy calculation through MM-GBSA assay
The energy calculation for the unbound receptor and receptor-ligand complex was estimated by using MM-GBSA analysis, and the results are depicted in Figure 3A. The results show that binding free energy of the protein-ligand complex was highly negative (-41,338.2 kcal/mol) in comparison to the unbound protein (-41,241.2 kcal/mol).

Prediction of pIC50
The IC 50 of lead B. frutescens phytochemicals were predicted on the basis of reported experimental IC 50 of P-glycoprotein standard inhibitors. The predicted IC 50 (pIC 50 ) values for the lead phytochemicals were between 15.08 and 54.67 µM. The results are comparable to the range of IC 50 of standard compounds (0.015-127.30 µM) [ Figure 3B]. A weak positive correlation (R 2 = 0.2709) was observed among the pIC 50 and binding energy of the test compounds.

Molecular dynamics simulation
P-glycoprotein is a membrane-bound protein, thus the protein was incorporated in a DPPC membrane during the process of MD simulation run [ Figure 4]. The Protein-ligand complex was subjected to MD simulation to analyze the dynamics of the binding mechanism of the ligand at the P-glycoprotein active site under clearly expressed water condition. The conformational changes analysis of the protein-ligand complex was analyzed during the last 5 ns of the 10-ns MD production run. RMSD is a fundamental method to analyze structural stability. A decreased RMSD of the complex with respect to the initial structure was observed [ Figure 5A]. The RMSD plot revealed that binding of ligand stabilized the P-glycoprotein structure such that the complex showed less deviation and conformational changes    Figure 5A]. The radius of gyration (R g ) is used to analyze the compactness of the protein-ligand complex and unbound protein structure in a biological system. The results show that the Rg for P-glycoprotein-ligand complex decreased along the simulation time, showing an increase in the compactness of the structure [ Figure 5B]. The RMSF of the P-glycoprotein and ligand complex was plotted to visualize the average fluctuation of all the amino acid residues during the simulation process [ Figure 5C]. The RMSF plot indicates that amino acid residue fluctuations are present in the protein during ligandbound state at several times. The results depict that interaction of ligand and protein brings both subunits (two arms of the protein) closer to each other and leads to the reduction in the distance between them [ Figure 6A-D].
Peptide-water hydrogen bond examination was conducted to determine the atomistic interaction between the protein and water in the presence and absence of the ligand (4'-Demethylknipholone 2'-β-Dglucopyranoside) during the simulation period. In the protein-ligand complex, there is a decrease in hydrogen bond number between protein and water compared to protein alone [ Figure 7A]. Peptidepeptide hydrogen bond analysis was accomplished to find the interaction between the protein residues. During the analysis, overall, it was found that ligand-protein interaction significantly increased the number of hydrogen bonds [ Figure 7B]. The SASA of the protein was calculated during the MD simulation in  Figure 7C]. The analysis indicates the folding states of protein and its stability upon ligand binding.

DISCUSSION
The failure of cancer therapy largely arises due to the multidrug resistance potential of the cancer cells which occurs through the overexpression and/or hyperactivation of the ABC transporters. The alteration in molecular mechanisms and gene regulation often promotes overexpression of P-glycoprotein (ABC transporter) in tumors. Overexpression of ABC transporters limit the anticancer drug response in cancer patients by reducing the accumulation of the chemotherapeutic drug in cancer cell/tumor through an efflux mechanism. Plasma membrane glycoprotein (Pgp) has been known to employ resistance to a variety of cytotoxic drugs including anticancer drugs such as vinblastine, doxorubicin, indinavir, ritonavir, and paclitaxel [4] . Chemotherapy failure in different types of cancer occurs due to MDR, which enables failure of cancer treatment in ~90% of cancer patients [22] . Development of ABC transporter inhibitors might be one tactic to overcome multidrug resistance through sensitization of cancer cells for chemotherapeutic drugs. P-glycoprotein inhibition may also increase the drug load in the cells and increase excessive toxicity.
In the present in silico study, we employed B. frutescens phytochemicals against ABC transporter protein.
Following the in-silico analysis, compound 4'-Demethylknipholone 2'-β-D-glucopyranoside was identified as a potent natural inhibitor of P-glycoprotein transporter. This is the first report on the Knipholone group of compounds acting as potent inhibitors for ABC transporter. Aller et al. [23] revealed that hydrophobic, aromatic, polar, and charged residues of the drug binding pocket provide poly-specificity to the P-glycoprotein toward the drug. They also showed that amino acid residues Phe332, Ile336, Phe339, Tyr949, Phe974, Val978, Ala981, Leu64, Leu335, Leu971, and Phe728 play a crucial role in the substrate/ inhibitor-transporter interaction. They showed that V978 plays an important role in binding with the standard inhibitors (verapamil and QZ59) of the protein. Interestingly, in the present study, we found that the lead molecule (4'-Demethylknipholone 2'-β-D-glucopyranoside) showed interaction with key amino acids including V978. The IC 50 of the identified B. frutescens phytochemical inhibitors were also predicted.
MM-GBSA analysis also revealed that lead B. frutescens phytoconstituent 4'-Demethylknipholone 2'-β-D-glucopyranoside forms a tight complex with ABC transporter and was stable in comparison to unbound ABC transporter protein. Previous studies have also utilized MMGBSA analysis to evaluate the ligand binding potential with the ABC transporter [24] . The stabilization of ligand-protein complex can be analyzed thorough MD simulations In addition, this approach evaluates conformational changes and stability of the three-dimensional structure of the protein. The present MD simulation study showed tight binding potential of the ligand, which provided stability to the ABC transporter-ligand complex. Our study corroborates the previous study performed by Shahraki et al. [25] , which showed that 1,4-Dihydropyridines inhibits ABC transporter by providing stability to the ABC transporterligand complex. The major limitations of the study include the isolation and characterization of the lead molecule 4'-Demethylknipholone 2'-β-D-glucopyranoside from the B. frutescens plant. In conclusion, the present study highlighted that the compound 4'-Demethylknipholone 2'-β-D-glucopyranoside possesses remarkable inhibitory potential against ABC transporter. Further, preclinical in vitro and in vivo validation is required to establish the B. frutescens phytoconstituents as natural ABC transporter inhibitors. Altogether, the identified lead inhibitor in the present study might serve as a starting molecule in the chemotherapy resistance drug discovery in the future.

Authors' contributions
Conceptualized and designed the study and prepared the draft manuscript: Kumar S Performed the experiments and analyzed data: Kushwaha