Novel natural inhibitors targeting B-RAF(V600E) by computational study

ABSTRACT The aim of this research was to screen the ZINC15 database to select lead compounds and drug candidates which can inhibit B-RAF (V600E). In order to identify drugs potentially inhibited B-RAF (V600E), numerous modules of Discovery Studio 4.5 were employed. Structure-based screening using LibDock was carried out followed by ADME (absorption, distribution, metabolism, excretion) and toxicity prediction. CDOCKER was performed to demonstrate the binding affinity and mechanism between ligands and B-RAF(V600E). To evaluate whether ligand-receptor complexes were stable, molecular dynamics were employed. Two novel natural compounds (ZINC000100168592 and ZINC000049784088) from ZINC15 database were found binding to B-RAF(V600E) with more favorable interaction energy in comparison with the reference drug Vemurafenib. Also, they were predicted with less ames mutagenicity, rodent carcinogenicity, non-developmental toxic potential and tolerance to cytochrome P450 2D6 (CYP2D6). The molecular dynamics simulation analysis indicated that the compound-B-RAF(V600E) complexes had more favorable potential energy compared with Vemurafenib and they can exist in natural environments stably. The result of this study shows that ZINC000100168592 and ZINC000049784088 are ideal leading potential compounds to inhibit B-RAF(V600E). The findings of this study and these selected drug candidates greatly contributed to the medication design and improvement of B-RAF(V600E) and other proteins.


Introduction
Melanoma is a kind of skin cancer of great aggressiveness and it has significant morbidity and mortality [1,2]. It has been previously observed that the malignant proliferation and transformation of melanocytes appeared in many parts of the body (skin, uvea, gastrointestinal mucosa, genitourinary mucosa, and CNS) [3]. In spite of its low proportion (1%) of all skin cancers, it became the leading cause of death in this group [4,5]. In different stages, the survival rate for patients with melonoma varies, the five-year survival of stages I and II is 98%, stage III is 64%, stage IV is merely 23% [6]. Surgery is the treatment of early-stage patients, which is highly curable [7]. Before 2011, chemotherapy was the only option available for many patients with advanced and unresectable melanoma [8]. Recent years, there have been numerous adjuvant approvals such as immunotherapy and targeted therapy approaches [7,8]. However, the prognosis of the patients is not good. Recently published data from the 5-year pooled analysis of COMBI(Combination therapy)-d and COMBI(Combination therapy)-v trials of dabrafenib and trametinib report an overall survival rate of 34% [9]. One of the reasons for the poor prognosis is that certain enzymes in the body prevent cancer drugs killing tumor cells through certain mechanisms.
The RAF family comprises A-RAF, B-RAF, and C-RAF [10]. In contrast to C-RAF and A-RAF, the basal kinase activity of B-RAF is the highest, and B-RAF is the most frequently activated by somatic mutations in malignancies (approximately 7%-9%) [11][12][13]. Besides, B-RAF is the most potent activator of MEK (Mitogen-activated protein kinase) in the RAS-MAPK signaling pathway which plays an important role in cell differentiation, proliferation and survival [10,11]. Furthermore, nearly 80% of the mutations of B-RAF observed is B-RAF (V600E) mutation [12]. Thus, B-RAF(V600E) plays a vital role in cancer targeted therapy.
Through hydrophobic interactions between the activation segment and the negatively charged regulatory region of B-RAF, which block the docking of ATP and the substrate, B-RAF keeps its inactive conformation [13]. When the two key sites of the activation segment (T599 and S602 for B-RAF) phosphorylate, negative charge will be provided. Given that the negative charge disrupts these hydrophobic interactions, the protein will be enzymatically active [13,14]. Study indicates that V600E mutation can mimic the phosphorylation at T599/S602 (phosphomimetic mutation), which can provide negative charge and have the B-RAF activated. As a result, B-RAF(V600E) has potent oncogenic activity [3,13,14].
In summary, the B-RAF(V600E) mutation can activate B-RAF which is the most potent activator of MEK (Mitogen-activated protein kinase) in the RAS-RAF-MEK-ERK signaling pathway. When the pathway is continuously and abnormally activated, cell proliferation, cell differentiation and cell viability will be abnormal, which may lead to the occurrence of tumor. B-RAF(V600E) mutations are present in hairy-cell leukemia, cutaneous melanoma, thyroid carcinomas and, less commonly, in ovarian, colon, lung, and other malignancies [15]. Hence, effective complexes targeting B-RAF (V600E) needed to be developed. Vemurafenib is a kind of drug of great selectivity targeting B-RAF (V600E). Because of its improving therapeutic effect on melanoma patients with B-RAF(V600E) mutations, rapid tumor response, low toxicity profile and tolerance, it quickly received the approval of FDA in 2011 [16]. It inhibits B-RAF kinase to decrease the downstream activation of MAPK pathway, resulting in the cell cycle arrest and apoptosis of B-RAF(V600E) mutation cells increased [17]. Although the efficacy of vemurafenib has been confirmed clinically, its durability is still a serious problem since many patients develop resistance to it after 6-8 months [18]. Hence, this study aimed to screen natural compounds from natural drugs that are more effective in inhibiting B-RAF(V600E).
Natural products, as lead compounds, can be transformed into new drugs through appropriate structural modification, which is an important source of new drug research in pharmaceutical industry. Paclitaxel and fungal metabolites and/or their analogues were good examples of natural products used in anti-tumor therapy [19,20]. Paclitaxel has been widely used in the treatment of breast cancer, and a number of fungal metabolites and/or their analogues such as anguidine, aphidicolin, fumagillin, illudins, irofulven, rhizoxin, wort-mannin, plinabulin (NPI-2358, a semisynthetic analogue of Phenyl-ahistin) and sonolisib (PX-866, a synthetic analogue of wortmannin) have progressed to various stages of cancer clinical trials [21]. In recent years, there has been an increasing interest in developing drugs targeting B-RAF(V600E). This prospective study was designed to screen natural compounds to verify potential inhibitors of B-RAF(V600E).
Recently, more and more emerging technologies have been integrated into the research of cancer detection and treatment [22,23]. In order to find new potential B-RAF(V600E) inhibitor, our study carried out a virtual screening against the Natural Products database (NP) in the ZINC database by using multiple modules of Discovery Studio 4.5. We analyzed the ADME (absorption, distribution, metabolism, excretion) and toxicity properties of molecules. By docking, we also analyzed binding modes and interactions between potential compounds and B-RAF(V600E). To assess if their binding interactions are stable, we performed a molecular dynamics simulation. This study lays groundwork for future clinical trials of these identified compounds and these identified novel natural compounds with structural modifications can be potential contributors for leads in further rational drug design for targeting B-RAF(V600E). The drug screened in this study is the first to be used in the inhibitory study of B-RAF(V600E), which will generate fresh insight into finding potential compounds to inhibit B-RAF (V600E). Besides, the method we used can also be helpful for the research of other drugs.

Discovery studio software and Ligand library
Discovery Studio 4.5 software (BIOVIA, San Diego, California, USA) is a suite of software for simulating small molecule and macromolecule systems; it is a new generation of molecular modeling and environmental simulation software for the life sciences field [22]. It aims to provide protein modeling, optimization, and drug design tools by applying protein structure and structural biologic computation. Numerous lead compounds and drug candidates were identified and refined through this method. LibDock module of Discovery Studio was employed for virtual screening; CDOCKER module was used for docking study; and ADME module was analyzed for pharmacologic properties. The Natural Products database in the ZINC15 database was selected to screen B-RAF(V600E) inhibitors. The ZINC15 database is a free database of commercially available compounds provided by the Irwin and Shoichet Laboratories, Department of Pharmaceutical Chemistry, University of California, San Francisco (San Francisco, California, USA) [24] Use LibDock for structure-based virtual filtering Ligand-binding pocket region of B-RAF(V600E) was selected as the binding site to screen compounds that could potentially inhibit B-RAF (V600E). Virtual screening was carried out using the LibDock module of Discovery Studio 4.5 [25]. LibDock is a rigid-based docking module. It calculates hotspots for the protein using a grid placed into the binding site and polar and a nonpolar probe. The hotspots are further used to align the ligands to form favorable interactions. The Smart Minimizer algorithm and CHARMM force field (Harvard University, Cambridge, Massachusetts, USA) were performed for ligand minimization [26]. After minimization, all ligand poses were ranked based on the ligands score. The 2.0 Å crystal structure of human B-RAF(V600E) and the inhibitor Vemurafenib were downloaded from the Protein Data Bank and imported to the working circumstance of LibDock. The molecule structure of B-RAF(V600E) was shown in Figure 1. The protein was prepared by removing crystal water and other heteroatoms around it, followed by addition of hydrogen, protonation, ionization, and energy minimization. The CHARMM force field and the Smart Minimizer algorithm were applied for energy minimization. The minimization performed 2000 steps with a root mean square gradient tolerance of 12.277, and the final root mean square gradient was 0.690. The prepared protein was employed to define the binding site. Using the ligands B-RAF(V600E) binding position, the active site for docking was generated. Virtual screening was performed by docking all the prepared ligands at the defined active site using LibDock. Based on the LibDock score, all the docked poses were ranked and grouped, and all compounds were ranked according to the LibDock score.

ADME (Absorption, Distribution, Metabolism, and Excretion) and Toxicity Prediction
The ADME module of Discovery Studio 4.5 was employed to calculate absorption, distribution, metabolism, and excretion (ADME) of selected compounds, including their aqueous solubility, blood-brain barrier penetration, cytochrome P-450 2D6 (CYP2D6) inhibition, hepatotoxicity, human intestinal absorption and plasma protein binding level [27]. TOPKAT module of Discovery Studio 4.5 was employed to calculate the toxicity and other properties of all the potential compounds, such as U.S. National Toxicology Program rodent carcinogenicity, Ames mutagenicity, developmental toxicity potential, and rat oral median lethal dose (LD50) and chronic oral lowest observed adverse effect level (LOAEL). These pharmacologic properties were fully considered when selecting proper drug candidates for B-RAF (V600E)

Analysis of ligand binding and ligand pharmacophore
On the basis of the CHARMM force field used for both receptors and ligands, CDOCKER, a molecular docking module of the software, was employed for the study of docking, providing docking results with high precision [28]. During the docking process, the receptor is held rigid, while ligands are allowed to flex, besides, the interaction energy and the CHARMM energy of the complexes can be calculated. Crystal structure of B-RAF(V600E) was obtained from the protein data bank. Our study removed crystalline water molecules in semi-flexible and rigid docking, which may result in the conformation of the receptor ligand complex to be affected by the fixed water molecules. Then, the protein was hydrogenated. For the purpose of verifying if the combination was reliable, Vemurafenib was removed from the binding site of B-RAF(V600E) and re-docked into it in order to compare the root-mean-square deviation (RMSD) of these two conformations. The binding site sphere of B-RAF (V600E) was defined as the regions that come within 5-Å radius from the geometric centroid of the ligand Vemurafenib. The structures of identified hits were prepared and docked into the binding pocket of B-RAF(V600E). Based on the interaction energy CODOCKER module provided, different positions of each ligand-B-RAF(V600E) complex were analyzed. 3D-QSAR, one of modules of Discovery Studio 4.5, was employed to generate the pharmacophores of the compounds and their pharmacophores were displayed.

Molecular Dynamics Simulation
The best binding conformations of the ligand-B-RAF(V600E) complexes among the poses predicted by the molecule docking program were submitted to MD simulation using Discovery Studio 4.5 [29]. The ligand-receptor complex was put into an orthorhombic box and solvated with an explicit periodic boundary solvation water model. In order to simulate the physiological environment, sodium chloride was added to the system with an ionic strength of 0.145. Then the system was subjected to the CHARMM force field which was used for ligand parameterization based on analogy. For the system, the following simulation protocols were applied: 1000 steps of minimization by steepest descent and conjugate gradient; 5ps-equilibration simulations in temperature of 300 K (slowly driven from initial temperature of 50 K for 2ps) and normal pressure ensemble; 25ps-MD simulation (production module) under NPT (normal pressure and temperature). The Particle Mesh Ewald (PME) algorithm was used to calculate long-range electrostatics, and the linear constraint solver (LINCS) algorithm was adapted to fix all bonds involving hydrogen. With initial complex setting as a reference, a trajectory was determined for root mean-square deviation(RMSD), potential energy, and structural characteristics through the Discovery Studio 4.5 analysis trajectory protocol in Discovery Studio 4.5.

Virtual Screening of Natural Products Database Against B-RAF(V600E)
Ligand-binding pocket was a critical regulatory site of B-RAF(V600E). When the ligand pocket was connected with some inhibitor, V600E could not give B-RAF negative charges to break the hydrogen interactions, thus preventing V600E from activating B-RAF, therefore, this pocket region was selected as a reference site. A total of 17,799 purchasable natural named product molecules were taken from the ZINC15 database. Molecular structure of B-RAF(V600E) was selected as the receptor protein. Vemurafenib, one of B-RAF(V600E) inhibitors, was chosen as a reference compound to evaluate the binding ability of other compounds. 7763 compounds were identified to bind with STING stably by the LibDock algorithm. The top 20 compounds with higher LibDock scores than were listed in Table 1.

ADME and Toxicity Prediction
Pharmacological properties of all selected ligands and Vemurafenib were first predicted by ADME module of Discovery Studio 4.5, including brain/ blood barrier (BBB), human intestinal absorption, aqueous solubility, cytochrome P450 2D6 (CYP2D6) binding, hepatotoxicity and plasma protein binding properties (PPB) ( Table 2). The aqueous solubility prediction (defined in water at 25°C) indicated that eight compounds were soluble in water. For human intestinal absorption, five compounds had a good absorption level and two compounds and Vemurafenib had a moderate absorption level. Eight compounds and Vemurafenib were found to be highly bound with plasma protein and the rest were just opposite. Seventeen compounds and Vemurafenib were predicted to be non-inhibitors of cytochrome P450 2D6 (CYP2D6), which was one of the important enzymes involved in drug metabolism. For hepatotoxicity, 11 compounds were predicted as nontoxic, which was similar to Vemurafenib Safety was also fully investigated in this study. To examine safety of the selected compounds, different toxicity indicators of the compounds and Vemurafenib, including Ames mutagenicity (AMES), Rodent carcinogenicity (based on the U.S. National Toxicology Program (NTP) dataset) and developmental toxicity potential (DTP) properties, were predicted using TOPKAT module of Discovery Studio 4.5 (Table 3). Results showed that 12 compounds were non-mutagenic. Considering all the results above, ZINC000100168592 and ZINC000049784088 were identified as ideal lead compounds, which were not CYP2D6 inhibitors thereby without hepatotoxicity. Moreover, they were predicted with less ames mutagenicity, rodent carcinogenicity and developmental toxicity potential compared with other compounds, which also strongly a Aqueous-solubility level: 0 (extremely low); 1 (very low, but possible); 2 (low); 3 (good) b Blood Brain Barrier level: 0 (Very high penetrant); 1 (High); 2 (Medium); 3 (Low); 4 (Undefined) c Cytochrome P450 2D6 level: 0 (Non-inhibitor); 1 (Inhibitor) d Hepatotoxicity: 0 (Nontoxic); 1 (Toxic) e Human-intestinal absorption level: 0 (good); 1 (moderate); 2 (poor); 3 (very poor) f Plasma Protein Binding: 0 (Absorbent weak); 1 (Absorbent strong) suggested their prospective application in drug development. In summary, ZINC000100168592 and ZINC000049784088 were identified as safe drug candidates and selected for following research (Figure 2).

Analysis of Ligand Binding
To study the two selected compounds and Vemurafenib's ligand-binding mechanisms with B-RAF(V600E), the CDOCKER module was employed. In comparison with ZINC000100168592 and ZINC000049784088, the CDOCKER potential energy of the reference ligand Vemurafenib was significantly higher. These indicated that the binding affinity of B-RAF(V600E) with ZINC000100168592 and ZINC000049784088 was higher than with Vemurafenib (Table 4).  (Tables 5 and 6).

Molecular Dynamics Simulation
A molecular dynamics simulation module was built to assess the stability of the ligand-B-RAF (V600E) complexes under natural environmental circumstances. The original conformations were obtained from the CDOCKER module through the molecular docking experiment. Figure 5 shows potential energy and RMSD curves chart of each complex. The time when the trajectories of each complex reached equilibrium was 22.5ps. Potential energy and RMSD of these complexes reached stable state as time goes by. In summary, ZINC000100168592 and ZINC000049784088 could interact with B-RAF(V600E), and the complexes existed stably in the natural environment.

Discussion
Melanoma is an aggressive skin cancer with significant morbidity and mortality [1]. Oncogenic mutations in B-RAF are found in approximately 40% of patients with cutaneous melanoma and activate the MAP kinase pathway [10]. Nearly 80% of the mutations of B-RAF observed is B-RAF(V600E) mutation [15]. Once the mutation of V600E occurs, it will lead to the continuous abnormal activation of B-RAF, which will lead to the disorder of the regulation of RAS-MAPK signaling pathway and increase the possibility of cancerization. Thus, B-RAF(V600E) plays a vital role in cancer targeted therapy. Despite the great progress in the drug design and development of B-RAF(V600E), only one drug, Vemurafenib, which is the reference drug of this study, received the approval of FDA in 2011 [30]. However, a high incidence of eventual resistance limited the effect of Vemurafenib [31]. Hence, screening more compounds targeting B-RAF(V600E) is urgent.
In order to screen out more potential drug candidates that can inhibit B-RAF(V600E), four modules of the software (Discovery Studio 4.5), including LibDock, ADME/TOPKAT, CDOCKER and Molecular Dynamics Simulation, were employed to screen and analyze the structural biological properties of novel potential compounds, respectively. Molecular conformation, pharmacological properties, binding affinity and stability were also fully analyzed to determine superiority of the selected compounds. In total,  In order to analyze the pharmacologic properties of these obtained compounds to reduce the waste of resources and improve the success rate of drug research in later drug development, toxicity predictions and ADME were processed. Results showed that ZINC000100168592 and ZINC000049784088 were identified as ideal lead compounds. They were all soluble in water and   also had a good absorption level. Meanwhile, they were non-inhibitors of cytochrome P450 2D6 (CYP2D6), which indicated they did not have hepatotoxicity. Additionally, these three compounds were also predicted with less ames mutagenicity, rodent carcinogenicity and developmental toxicity potential compared with other compounds, which also strongly suggested their prospective application in drug development. On the other hand, the rest drug in the list also had potential application in drug development even though they possessed toxicity, since specific groups and atoms could be added to reduce its toxicity. Considering all the results above, ZINC000100168592 and ZINC000049784088 were selected as ideal lead compounds and further analysis would be carried out.
To conduct more accurate docking analysis of candidate compounds and investigate chemical bonds and ligand-binding mechanisms of B-RAF (V600E) with candidate compounds, the CDOCKER module was employed. Computation results illustrated that CDOCKER the interaction energy of the reference ligand Vemurafenib is higher than ZINC000100168592 and ZINC000049784088, which explains why the binding affinity of ZINC000100168592 and ZINC000049784088 with B-RAF(V600E) is higher than with Vemurafenib. Then, the molecular structure inspection was used to analyze the chemical structure of ZINC000100168592 and ZINC000049784088. Because there were more chemical bonds between these 2 compounds and B-RAF(V600E) in comparison with Vemurafenib, the compounds could bind more stably with B-RAF(V600E), which may strengthen their inhibition of B-RAF(V600E) and thus play a role in improved killing of tumors.
Finally, molecular dynamics simulation was employed to evaluate whether they can maintain stability in the natural environment. The calculation results of potential energy and RMSD of these ligand-B-RAF(V600E) complexes illustrated that the time at which the trajectories of complexes reached equilibrium is 22.5 ps. With time passing, potential energy and RMSD of these complexes reached a stable state which manifested that these two complexes could keep stable in the natural environment. Prospectively, the refinement and modification of the drugs could be proceeded in order to make ligands and receptors bind in a more stable state on the basis of the results above. Because of the provision of high stabilization and binding affinity of the selected drugs in this research, development of inhibitive durg of B-RAF(V600E) was strengthened. Table 6. π-π interaction, π-Alkyl interaction, Alkyl interaction, π-cation interaction and Unfavorable Doner-Doner interaction, parameters for each compound and mTORC1 residues. In spite of the fact that overall research on oncology drugs has not progressed well, development and design of them are still a research hotspot. This study indicated that screening more ideal lead compounds is the most important in current development and design of oncology drugs. In order to screen the structure of novel natural compounds and select potential drug candidates, numerous modules of the software (Discovery Studio 4.5) were employed in this study. Results demonstrated that the binding affinity, molecule conformation and stability of each selected compound were better than the reference compound Vemurafenib, besides, these two compounds may have huge therapeutic potential in B-RAF(V600E)-related cancer.
In summary, this study aimed to find more potential drug candidates which can inhibit B-RAF(V600E) form a natural compounds database. Even though this study was conducted by elaborate design and precise measurements was performed, we still admitted that there were still some limitations in this study. Without being improved and refined, no drug can be marketed. To make these two compounds more perfect as drug candidates, some groups and atoms which can influence the pharmacological properties of the drugs are needed to modify. More experiments need to be performed to validate our results and more indicators regarding to drugs safety, such as MTD (Maximum Tolerated Dosage) and AB(Aerobic Biodegradability) should also be assessed in our future study. These limitations are directions and focus of our further study.

Conclusion
This study performed a series of computeraided structural and chemical analysis technology (including virtual screening, molecule docking, ADME, toxicity prediction) to screen and identify the ideal-leading compound with functions to inhibit B-RAF(V600E). Two compounds, ZINC000100168592 and ZINC000049784088, were predicted as potential inhibitors targeting B-RAF(V600E). The experiments carried out by Discovery Studio4.5 confirmed that these two compounds could bind tightly with B-RAF(V600E) in the natural environments simulated by the MD simulation module. In spite of the limitations, this study lays the groundwork for future clinical trials of these two compounds. Moreover, these novel natural compounds with structural modifications can be potential contributors for leads in further rational drug design for targeting B-RAF(V600E).