Molecular modeling, dynamics simulations, and binding efficiency of berberine derivatives: A new group of RAF inhibitors for cancer treatment

RAF kinases are a family of enzymes in the MAP kinase pathway that contribute to the development of different types of cancer. BRAF is the most important member of RAF kinases. BRAF mutations have been detected in 7% of all cancers and 66% of melanomas; as such, the FDA has approved a few BRAF inhibitor drugs to date. However, BRAF can activate CRAF leading to resistance to BRAF inhibitors. Berberine (BBR) is an alkaloid that is widely distributed in different plant species. Several studies have been carried out on the anti-cancer effects of BBR but direct targets of BBR are unknown. In this study, interactions of BBR derivatives against BRAF and CRAF kinases were modeled and predicted using an in silico-based approach. To analyze and identify the residues important in BRAF docking, we modeled interactions of ATP, the universal substrate of BRAF, and found that Lys483 and Asp594 are the most important residues involved in both ATP and BBR binding [(The average score = -11.5 kcal/mol (ATP); Range of scores = -7.78 to -9.55 kcal/mol (BBR)]. In addition to these polar residues, Trp530 and Phe583 are also applicable to the molecular docking of BRAF. We also observed that Asp593 was excluded from the enzyme cavity, while Phe594 was included inside the cavity, making the enzyme inactive. Finally, three alternatives for BBR were identified with dual RAF inhibition effects [The best scores against BRAF = -11.62 kcal/mol (BBR-7), -10.64 kcal/mol (BBR-9), and -11.01 kcal/mol (BBR-10); the best scores against CRAF = -9.68 kcal/mol (BBR-7), -9.60 kcal/mol (BBR-9), and -9.20 kcal/mol (BBR-10)]. Direct effects of BBR derivatives against BRAF and CRAF kinases had not yet been reported previously, and, thus, for the first time, we report three cycloprotoberberines as lead compounds against RAF kinases.


Ligand selection and molecular docking
Using PubChem database, we virtually screened all compounds similar to berberine (Pub-Chem CID: 2353). We found 1,544 compounds, which were sorted based on chemical properties, such as molecular weight, H-bond donor, H-bond acceptor, formal charge, and total formal charge ( Table 1).
Lipinski's rule of five was verified for 1,544 derivatives, among which 485 ligands passed the test [24]. These 485 compounds were chosen and their SDF files were downloaded and prepared for docking. In addition to BBR derivatives, vemurafenib and ATP-Mg were used as controls in the present study. ACD/ChemSketch 12.01 was used for in silico ligand generation, and pdb-formatted files were prepared for all compounds. To prepare and optimize the ligands for docking, polar hydrogen atoms were inserted, torsional degrees of freedom (nTDOF) were determined, and Gasteiger charges were calculated for all generated ligands. Final suitable berberine derivatives and controls used are depicted in Fig 1. The top 10 ligands ranked based on thermodynamics parameters were shown here.
Five crystal structures were downloaded for the BRAF kinase domain (PDB IDs: 1uwj, 1uwh, 3c4c, 3og7, and 3psd) from the RCSB protein data bank, and all generated ligands were docked to all BRAF crystals; however, a single CRAF kinase domain was found on the RCSB protein data bank (PDB ID: 3omv) [25]. Water molecules and original ligands were removed from the downloaded pdb files ( Table 2), and the crystals were optimized and energetically minimized by Swiss-PdbViewer 4.1.0. Optimized crystals were used only for docking the generated ligands. AutoDock 4.2 was used to dock all ligands [26] (the grid and docking configuration set ups are available in the online edition of this article). To set up the grid box, points were adjusted to 60 in each dimension and optimized coordinates of the CA atoms of Lys482 (PDB IDs: 1uwh and 1uwj), Lys483 (PDB IDs: 3c4c, 3og7, and 3psd), and Lys375 (PDB ID: 3omv) were determined as the center of the grid box in the corresponding docking.

Ligand topology and molecular dynamics simulations
Hydrogen atoms were added to ligands evaluated for molecular dynamics simulation. We carried out semi-empirical quantum mechanical calculations (QM) for ligands with atom numbers of more than 40 (BBR, BBR-7, BBR-9, BBR-10, sorafenib, and vemurafenib) to create ligand topology files as input files required for molecular dynamics simulations. Bonded parameters, such as bond length, bond angles, and dihedral angles, as well as the initial partial atomic charges, were calculated for the united-atom version of the 53a6 GROMOS force field. Accordingly, the parameters and topology files were generated and optimized using the Automated Force Field Topology Builder (ATB) version 2.2 [30,31]. All molecular dynamics (MD) simulations were carried out using GROMACS v.5.0.4 under the Ubuntu Linux platform version 15.0.4 with a GROMOS 53a6 force field [32]. Water systems were generated such that each docked protein was placed in the center of a triclinic box with a minimum distance between the solutes and the edge of the box of 1.0nm. Accordingly, all solutes were solvated using the SPC water model [32]. Adding sodium and chloride ions neutralized the systems. All configurations used were based on previously published protocols [2,33,34].

Molecular visualization
To demonstrate inter-molecular interactions (e.g., hydrophobic, h-bonds, halogen bonds, and π/aromatic interactions), the Accelrys Discovery Studio Visualizer version 4.1 (DS) was applied. In addition to DS, intermolecular hydrogen-bonds were also checked using the Lig-Plot + v. 1 [26]. Using the UCSF Chimera and DS, we added all hydrogen bonds and any editing required for ligand topology generations. Molecular dynamics trajectories were visualized and plotted using Grace 5.0.0.

Molecular docking of BBR derivatives on BRAF
To analyze the effects of BBR derivatives on BRAF, we selected five pdb-formatted crystal structures. Each PDB structure was built under slightly different X-ray crystallography conditions. Therefore, we repeated, compared, and calculated the means of binding energies and inhibition constants in five independent in silico experiments ( Table 3).
Three active BRAFs (PDB IDs: 3og7, 3psd, and 3c4c) and two inactive BRAFs (PDB IDs: 1uwh and 1uwj) were used to compare the effects of each ligand on different BRAF conformations. We also focused on active BRAFs and, specifically BRAF V600E , the well-known BRAF mutation (PDB ID: 3og7).
In Table 3, docking information is described. The number of conformational clusters was dependent upon the size of the ligands, the number of torsional degrees of freedom, n TDOF , and the conformation of the active site of the enzyme (Table 1). Each ligand requires free positive torsional energy to fit suitably with the BRAF active site. For molecules with higher n TDOF , torsional free energy is greater than the compounds with the lower n TDOF , which leads to diversity in ligand conformations entering the active site in separate iterations. BBR, BBR-4, and BBR-5 have the lowest torsional free energy (torsional ΔG = +0.60 kcal/mol), while BBR-2 has the highest torsional free energy (torsional ΔG = +1.79 kcal/mol).
Estimated free energy of binding is the sum of the following energies:

Torsional Free Energy
For the most potent BBR derivatives, the estimated free energy should be less than -9.00 kcal/mol. As torsional free energy is a positive parameter, it affects the total binding energy that should normally be at least possible or negative. To identify the optimal BRAF inhibitors, inhibition constants (Ki) were served in the form of pKi. pKi equals-LogKi, in which the unit of Ki (nM) is converted to molar (M). After calculating the mean Ki and pKi, the optimal BRAF inhibitors were chosen based on the highest pKi. BBR and most of its derivatives are effective in targeting BRAF, with pKis greater than 6; however, we generated three novel BBR derivatives, which were not reported previously in the PubChem database.
After comparing the results of docking and molecular properties, we added a trifluoromethyl benzyl group to BBR at three positions (9, 10, and 13 for BBR-9, BBR-7, and BBR-10, respectively) [35]. This group has two features: 1. Hydrophobicity is increased.
2. Halogen bonds were generated; the chance of interacting with lysine is increased.
After performing molecular docking algorithms for these new compounds, we found that pKi of BBR-9, BBR-7, and BBR-10 were increased to 6.81, 6.60, and 6.57, respectively, whereas the pKi for BBR was 6.08. Therefore, these three novel compounds are novel BRAF inhibitors and, as such, their results were compared to the molecular docking of vemurafenib and sorafenib as positive controls.

Molecular interactions
To identify the residues important to BRAF docking, we modeled the interactions of ATP, the universal substrate of BRAF (Fig 2). ATP binds to Lys483 and Asp594 [9]. In addition to these polar residues, Trp530 and Phe583 are other residues applicable to the molecular docking of BRAF. In our comparison between active and inactive forms, Asp593 was found to be located outside of the enzyme cavity and Phe594 was found to be inside the cavity, making the enzyme inactive (PDB ID: 1uwj).
The active site of the enzyme was chosen (selective approach) using the preliminary docking of ATP-Mg and molecular structure of the kinase domain. BBR-9 was predicted as a BRAF inhibitor, which is bound strongly to the kinase pocket. Residues Lys483, Phe583, Trp531, and Asp594 were recognized as the crucial residues located at the kinase domain of BRAF.
To validate whether BBR-9 is a selective BRAF-kinase inhibitor or non-selective, blind docking (randomized approach) was also performed. The minimum inhibition constants belonged to the conformations randomly docked to the kinase cavity to the residues involved in the selective approach. In Fig 2, BBR binding to both forms of BRAF were compared and show the similarity of BBR interactions in both BRAF conformations. BBR was bound to both active and inactive forms with similar residue involvement. In almost all iterations, Lys483 (or Lys482 of an 3c4c) and inactive (1uwj) BRAFs. LYS plays a key role in both forms of BRAF. Lys482 (1uwj) and Lys483 (3c4c) indicate H-bond and cation-π interactions; however, in 3c4c (BRAF active form), Phe594 is out of the active site. TRP531 and PHE583 are the major amino acids involved in the ππ hydrophobic interactions (π-π and cation-π interactions). ASP594 and PHE595 are located in and out of the active site, respectively, to show the active form. The results predicts that BBR binds to both active and inactive forms of BRAF.
https://doi.org/10.1371/journal.pone.0193941.g002 inactive form) was the most important residue involved in ATP-BRAF and BBR-BRAF interactions, and lysine, with its cation group ( + NH3), is involved in H-bond and cation-π interactions. Aromatic groups of BBR play a major role in hydrophobic interactions. Two aromatic residues, Trp531 (or Trp530 of inactive BRAF) and Phe583 (or Phe582 of inactive BRAF), increased binding stability by providing π-π interactions with BBR aromatic rings.
The three novel compounds generated had the lowest binding energies and inhibition constants and the highest pKi ( Table 3). Binding of BBR-7, BBR-9, and vemurafenib (Fig 3) to the continuously active form of BRAF V600E (PDB ID: 3og7) was compared. We concluded that Lys483 is the most important residue, while Trp531 and Phe583 form π-π interactions.
We added additional aromatic group to different positions of BBR to generate BBR-7, BBR-9, and BBR10; therefore, the binding energies were considerably reduced, as Phe583 and Trp531 formed additional π-interactions. Thus, the additional intermolecular π-π interaction makes these novel compounds better inhibitors than BBR. In addition to the intermolecular π-interaction, BBR-9 formed intramolecular π-π interactions (Fig 3E), which renders it the best-fitted compound with a more rigid docking structure compared with the other BBR derivatives.

Molecular docking of BBR-9 on CRAF
To find dual RAF inhibitors, the optimal BRAF inhibitors were selected based on pKi values and molecular docking analysis was performed for CRAF crystal (PDB ID: 3omv) with 50 iterations. Vemurafenib and ATP-Mg were also used as in silico positive controls for both BRAF and CRAF kinases. As pKi is a target-specific parameter, we observed that BBR-9 hds the highest pKi of 6.58; pKi was calculated based on the mean Ki of 50 iterations ( Table 4). Table 4 shows the comparative docking results of BRAF and CRAF. The structure of CRAF binding cavity is nearly the same as BRAF with Lys, Asp, and Phe. In Fig 3F and 3G, BBR-9 shows interactions with Lys375 and Asp486 of CRAF, and shows that Phe475 forms a face-to-face ππ interaction with the bended aromatic ring of BBR-9.
Based on Table 4, we predicted that BBR-9 and BBR-10 may be candidates as dual RAF inhibitors stronger than BBR. We then validated the stability of docking using MD analysis.

Computational validation by MD analysis
By comparing the results of docking between BBR derivatives and several protein structures, we predicted that BBR-9 and BBR-7 would be BRAF inhibitors and BBR-9 would be a CRAF inhibitor. To show the stability and system energies of complexes, 10 ns MD simulations were conducted using Gromacs. RMSDs of around 0.3nm, with RMS fluctuations of less than 0.2nm, showed the stability of the interactions (Figs 4 and 5). During MD simulations, vemurafenib was used as a positive control. As the experimental property of vemurafenib was already reported, this comparison was used to predict our novel candidates, BBR-7 and BBR-9, as BRAF inhibitors. Modifications to protein motifs were evaluated using the radius of gyration (Rg). Rg provides information on total protein volume distribution in spherical states and indicates molecular shape over time. For BRAF/BBR-9 interaction (Fig 4G-4L), Rg was higher than CRAF/BBR-9 (Fig 5).
To validate the interactions between RAF kinases and BBR-9, radial distribution functions (RDF or g(r)) were computed. RDF demonstrates the distance between each molecule's center of mass [36]. The g(r)s of BBR-9 and BBR-7 for BRAF were higher than that of BRAF bound to vemurafenib (Fig 4); however, the average distances (horizontal axis) between inhibitors and targets for all ligands were approximately 2Å, which shows the length of effective interactions. To confirm molecular docking, we computed the average number of H-bond and thermodynamic parameters of RAFs/ligand complexes ( Table 5).

BRAF binding efficiency indices of BBR derivatives
Based on docking results, Binding Efficiency Indices (BEIs) were calculated to optimize the thermodynamics of docking based on the physicochemical properties of each ligand. BEIs are enzyme-specific indices for each molecule. To calculate efficiency indices, we combined molecular properties and docking information to provide quantitative parameters of the specificity of each compound for the BRAF kinase. Then, chemical and docking interactions of different ligands were compared, and bivariate plots were created to determine those that most potently antagonized BRAF. To calculate the difference of BEIs presented in Table 6, we used formulas and definitions described by Abed-Zapatero [37][38][39][40].
Valuable chemical properties of drugs used in efficiency index calculations include molecular size, such as the number of atoms [non-hydrogen (nha) and electronegative] and   lipophilicity or polarity (e.g., TPSA, LogP etc.) of ligands affecting drug absorbance and distribution. For better absorbance, the TPSA (topological polar surface area) should be lower than 100Å 2 . The MW-based logarithmic BEI (mBEI) and atom-based logarithmic binding efficiency index (nSEI) for BBR-9 were 6.48 and 8.35, respectively. Using different molecular properties (MW, TPSA, and the number of atoms), Abad-Zapatero reported different BEIs based on different logarithmic and nonlogarithmic numbers of atoms; however, nBEI, mBEI, and NSEI are three parameters supported by Abad-Zapatero, by which the best compounds against a certain target can be predicted. This novel prediction is a combination of chemical structure and thermodynamics (docking) results [41]. While we generated bivariate plots based on nBEI versus NSEI, we found that BBR-9, BBR-7, and BBR-10 were located in the top right corner of each plot, indicating they were the most potentially predicted BRAF antagonists (Fig 6). In contrast to vemurafenib and sorafenib, of which the LLE index is higher, BBR-9, BBR-7, and BBR-10 have much higher NSEI than LLE, meaning that the highest polarity of these ligands, especially BBR-9, may have a higher impact on the interactions between ligand and polar residues such as Lysine.

Discussion
BRAF mutations have been detected in 7% of all cancers and in 66% of melanomas and 12% of colorectal cancers. As such, the FDA has approved BRAF inhibitors for cancer treatment. However, BRAF can alternatively activate CRAF leading to resistance to BRAF inhibitors. Here, for the first time, we analyzed anti-BRAF activity of BBR and demonstrated that it targets BRAF as a part of its anti-cancer activities.
RAF kinases are located upstream of the MAP kinase pathway and downstream of receptor tyrosine kinases (RTKs) (e.g., EGFR and Her2) and RAS. RTKs and RAS are mutated in 90% of pancreatic cancers, 30% of breast cancers, more than 50% of carcinomas, 30% of acute myelogenous leukemia, 30% of liver cancers, 55% of follicular thyroid cancers, 60% of undifferentiated papillary cancers, 35% of lung adenocarcinomas, and 10% of kidney cancers [13]. Potent RAF inhibitors can significantly block signals transduced by over-activated RTKs and RAS [42].
BBR attenuates the MAPK/ERK/RAF pathway, affecting cell proliferation in many types of cancer [43]. In addition, in vitro and in silico studies have shown that BBR reduces levels of ERK and p38 MAPK [44]; BBR also silences BRAF/ERK signaling in melanoma cells [8]. In the present study, we predicted that BRAF and CRAF were additional targets of BBR that affect ERK/p38 as downstream phosphate acceptors of BRAF, as BBR is bound to the position that ATP is normally bound to. Since BBR has cytotoxic effects, identifying new BRAF inhibitors with less cytotoxicity is of great importance to cancer treatment [1]. Position 13 of BBR plays a major role in reducing MAPK pathway activity. HWY 289 and HWY 336, two BBR derivatives with added aromatic branches on position 13, suppress cell proliferation in fission yeast (Schizosaccharomyces pombe) by inhibiting the MAPK cascade. The minimal inhibitory concentration (MIC) of HWY 289 and HWY 336 are 29.52μM and 11.83μM, respectively [35]. In contrast to BBR which does not completely inhibit proliferation of wild type S. pombe, HWY 289 and HWY 336 completely block the proliferation of this yeast [35]. After docking hundreds of BBR derivatives obtained from the PubChem database, we generated BBR-10 with the addition of [4-(trifluoromethyl) methyl benzene] to position 13. After docking BBR-10 against BRAF crystals, it was bound to the BRAF active site with an average Ki of 0.269μM, while the Ki of BBR was 0.832μM, three times higher than that of BBR-10. These results identified BBR-10 as a BBR-derived lead candidate for BRAF inhibition.
The structure-activity relationship (SAR) showed that replacing the methoxyl group at position 9 with an ester moiety could significantly increase the anti-tumor activity of BBR [45]. In addition to BBR-10, we added a trifluorobenzel group to methoxy group located at the positions 9 (BBR-9) and 10 (BBR-7). The trifluorobenzel structure is similar to the active part of sorafenib. Docking results of BBR-7 and BBR-9 against BRAF and CRAF kinases showed that they are both more effective than BBR-10, as BBR-9 was bound to BRAF and CRAF with Ki values of 0.155μM and 0.265μM, respectively. Our results also showed that BBR-9 might be effective as a dual RAF inhibitor. The direct effects of BBR derivatives against BRAF and CRAF kinases had not yet been reported, but, here, for the first time, we report three cycloprotoberberines as lead compounds against RAF kinases.

Conclusion
In the present study, we sought to identify drug candidates that simultaneously inhibit BRAF and CRAF kinases. We predicted that BBR-9 would be a strong dual RAF inhibitor. In addition to BBR-9, BBR-7 and BBR-10 were also identified as attractive candidates against RAF kinases. In vitro and in vivo evaluation of these inhibitors may establish them as lead compounds to treat cancer. In addition to direct RAF inhibitory of BBR derivatives, these compounds may act as indirect inhibitors of over-activated RTKs and RAS for different types of cancer.