An Investigation of Molecular Docking and Molecular Dynamic Simulation on Imidazopyridines as B-Raf Kinase Inhibitors

In the recent cancer treatment, B-Raf kinase is one of key targets. Nowadays, a group of imidazopyridines as B-Raf kinase inhibitors have been reported. In order to investigate the interaction between this group of inhibitors and B-Raf kinase, molecular docking, molecular dynamic (MD) simulation and binding free energy (ΔGbind) calculation were performed in this work. Molecular docking was carried out to identify the key residues in the binding site, and MD simulations were performed to determine the detail binding mode. The results obtained from MD simulation reveal that the binding site is stable during the MD simulations, and some hydrogen bonds (H-bonds) in MD simulations are different from H-bonds in the docking mode. Based on the obtained MD trajectories, ΔGbind was computed by using Molecular Mechanics Generalized Born Surface Area (MM-GBSA), and the obtained energies are consistent with the activities. An energetic analysis reveals that both electrostatic and van der Waals contributions are important to ΔGbind, and the unfavorable polar solvation contribution results in the instability of the inhibitor with the lowest activity. These results are expected to understand the binding between B-Raf and imidazopyridines and provide some useful information to design potential B-Raf inhibitors.


Introduction
In spite of the progress in medicine, developing new anticancer drugs is still important because cancer still acts as a major problem of health all over the world [1].According to the recent reports [2][3][4], MAPK pathway (also called Ras-Raf-MEK-ERK pathway) is very crucial for cell survival and proliferation because this pathway can be activated easily in human cancers (up to 30%).There are three isoforms (A-Raf, B-Raf, and C-Raf) for Raf kinase in this pathway [5], and because the mutations of B-Raf kinase in human cancers is up to 7%, they has been considered as the primary activator in this pathway.There is a different mutation frequency of B-Raf kinase in various human cancers, such as colorectal cancers (10%), thyroid cancers (30%), ovarian cancers (35%), and melanoma (50%-70%) [6].Thus, B-Raf kinase has been a key target in recent cancer treatment [7][8][9].
The first B-Raf kinase inhibitor (BRI) approved by the Food and Drug Administration (FDA) is Sorafenib, which is used to treat hepatocellular carcinoma and renal cell carcinoma in clinic.Vemurafenib is the second BRI approved by FDA, which is used to treat metastatic melanoma in clinic.Furthermore, some BRIs are in various stage of clinical development, such as RAF265, GSK2118436, and SB-590885.In spite of the success of clinical efficiency of these inhibitors in cancer treatments, they still have some major side effects and can develop drug resistance.Therefore, it is still important to develop other potent and selective BRIs [10].Recently, Newhouse et al. have synthesized a series of imidazopyridines as BRIs, which show excellent potency and selectivity.These BRIs bind in a DFG-in, αC-helix out, inactive conformation of wild-type B-Raf kinase [11].In our previous work, we performed 3D QSAR, pharmacophore modeling, and virtual screening studies on this series of molecules to help design more potential BRIs [12].In order to know the interaction between this series of inhibitors and B-Raf kinase, an investigation of molecular docking, MD simulation and ∆G bind calculation on this kind of BRIs were carried out in this work, in which three inhibitors (Mol 1, Mol 2, and Mol 3) were selected (Figure 1).The reason we chose these three molecules is that the crystal structure of B-Raf kinase combined with Mol 1 is available, Mol 2 shows the highest inhibitory activity, and Mol 3 shows the lowest activity.Their activities (IC 50 ) are 61 nM, 0.76 nM, and 167 nM, respectively [11].
Int. J. Mol.Sci.2015, 16, page-page 2 treatments, they still have some major side effects and can develop drug resistance.Therefore, it is still important to develop other potent and selective BRIs [10].Recently, Newhouse et al. have synthesized a series of imidazopyridines as BRIs, which show excellent potency and selectivity.These BRIs bind in a DFG-in, αC-helix out, inactive conformation of wild-type B-Raf kinase [11].In our previous work, we performed 3D QSAR, pharmacophore modeling, and virtual screening studies on this series of molecules to help design more potential BRIs [12].In order to know the interaction between this series of inhibitors and B-Raf kinase, an investigation of molecular docking, MD simulation and ΔGbind calculation on this kind of BRIs were carried out in this work, in which three inhibitors (Mol 1, Mol 2, and Mol 3) were selected (Figure 1).The reason we chose these three molecules is that the crystal structure of B-Raf kinase combined with Mol 1 is available, Mol 2 shows the highest inhibitory activity, and Mol 3 shows the lowest activity.Their activities (IC50) are 61 nM, 0.76 nM, and 167 nM, respectively [11].In present work, molecular docking was carried out between all the imidazopyridines and B-Raf kinase.MD simulations were performed between the three inhibitors (Mol 1, Mol 2 and Mol 3) and B-Raf kinase, and the analysis of hydrogen bond (H-bond) in MD simulations was performed as well.Furthermore, the ΔGbind was computed by MM-GBSA method based on MD trajectories, and the contributions to the ΔGbind were also analyzed.

Molecular Docking
To validate the docking method and docking accuracy, Mol 1, was docked into the binding site of B-Raf kinase receptor.Both Mol 1 ligand and B-Raf kinase receptor were isolated from the complex crystal structure (PDB code: 4MBJ) [11].The root-mean-square deviation (RMSD) between the docked and the crystal structures of Mol 1 was only 1.698 Å (less than 2 Å), which is satisfactory.Figure 2 shows that the docked structure (red color) and the X-ray crystal structure (green color) are quite similar.In addition, all the 36 imidazopyridines were docked into the binding pocket of B-Raf kinase receptor successfully.The chemical structures, biological activity values and docking C_scores of the imidazopyridines are shown in Table S1 (supplementary data).Almost all inhibitors show high In present work, molecular docking was carried out between all the imidazopyridines and B-Raf kinase.MD simulations were performed between the three inhibitors (Mol 1, Mol 2 and Mol 3) and B-Raf kinase, and the analysis of hydrogen bond (H-bond) in MD simulations was performed as well.Furthermore, the ∆G bind was computed by MM-GBSA method based on MD trajectories, and the contributions to the ∆G bind were also analyzed.

Molecular Docking
To validate the docking method and docking accuracy, Mol 1, was docked into the binding site of B-Raf kinase receptor.Both Mol 1 ligand and B-Raf kinase receptor were isolated from the complex crystal structure (PDB code: 4MBJ) [11].The root-mean-square deviation (RMSD) between the docked and the crystal structures of Mol 1 was only 1.698 Å (less than 2 Å), which is satisfactory.Figure 2 shows that the docked structure (red color) and the X-ray crystal structure (green color) are quite similar.In addition, all the 36 imidazopyridines were docked into the binding pocket of B-Raf kinase receptor successfully.The chemical structures, biological activity values and docking C_scores of the imidazopyridines are shown in Table S1 (supplementary data).Almost all inhibitors show high C_score values, which are more than 5.0.The correlation between C_score values and biological activity (pIC 50 values) of 36 imidazopyridines is shown in Figure 3.The above results indicate an acceptable reliability of the docking method for the B-Raf kinase receptor and these inhibitors.In order to illustrate the interactions between B-Raf kinase and imidazopyridine, we focus on receptor-ligand interactions between B-Raf kinase and Mol 1 (the ligand of 4MBJ), Mol 2 (the most active inhibitor), and Mol 3 (the least active inhibitor).Figure 4a shows the docking mode between B-Raf kinase and Mol 1, in which four H-bonds were formed: the first one is between the >C=O of CYS 532 and the H-N1 (Figure 1) of Mol 1 (>C=O•••H-N1) with a distance of 1.89 Å and a deviated angle of 15.9°; the second one is formed by amide hydrogen of TRP 531 and N4 of Mol 1 (-N-H•••N4, 2.20 Å, 18.2°); the third one is between the carbonyl oxygen of GLY 593 and H-N3 of Mol 1 (>C=O•••H-N3, 2.13 Å, 45.0°); and the fourth one is between the amide hydrogen of PHE 595 and sulphuryl oxygen atom of Mol 1 (-N-H•••O=S, 2.04 Å, 54.8°).Due to their large angles (45.0° and 54.8°), the third and fourth H-bonds show less stability than the first and second H-bonds.As shown in Figure 4a, there are a π-π stacking contact between the imidazolepyridine ring of Mol 1 and the aromatic ring of PHE583, and a hydrophobic interaction between the benzene ring of Mol 1 and the side chain of THR529.The above observations are consistent with the previous studies [13].
The docking mode between B-Raf kinase and Mol 2 can be seen in Figure 4b, in which three H-bonds were formed: the first one is between the >C=O of CYS 532 and the H-N1 of Mol 2 (>C=O•••H-N1, 1.96 Å, 10.3°); the second one is between the amide hydrogen of TRP 531 and N4 of Mol 2 (-N-H•••N4, 2.18 Å, 9.7°); the third one is between the amide hydrogen of PHE 595 and sulphuryl oxygen atom (-N-H•••O=S, 2.00 Å, 53.7°).The angle data of H-bonds indicates that the   In order to illustrate the interactions between B-Raf kinase and imidazopyridine, we focus on receptor-ligand interactions between B-Raf kinase and Mol 1 (the ligand of 4MBJ), Mol 2 (the most active inhibitor), and Mol 3 (the least active inhibitor).Figure 4a shows the docking mode between B-Raf kinase and Mol 1, in which four H-bonds were formed: the first one is between the >C=O of CYS 532 and the H-N1 (Figure 1) of Mol 1 (>C=O•••H-N1) with a distance of 1.89 Å and a deviated angle of 15.9°; the second one is formed by amide hydrogen of TRP 531 and N4 of Mol 1 (-N-H•••N4, 2.20 Å, 18.2°); the third one is between the carbonyl oxygen of GLY 593 and H-N3 of Mol 1 (>C=O•••H-N3, 2.13 Å, 45.0°); and the fourth one is between the amide hydrogen of PHE 595 and sulphuryl oxygen atom of Mol 1 (-N-H•••O=S, 2.04 Å, 54.8°).Due to their large angles (45.0° and 54.8°), the third and fourth H-bonds show less stability than the first and second H-bonds.As shown in Figure 4a, there are a π-π stacking contact between the imidazolepyridine ring of Mol 1 and the aromatic ring of PHE583, and a hydrophobic interaction between the benzene ring of Mol 1 and the side chain of THR529.The above observations are consistent with the previous studies [13].
The docking mode between B-Raf kinase and Mol 2 can be seen in Figure 4b, in which three H-bonds were formed: the first one is between the >C=O of CYS 532 and the H-N1 of Mol 2 (>C=O•••H-N1, 1.96 Å, 10.3°); the second one is between the amide hydrogen of TRP 531 and N4 of Mol 2 (-N-H•••N4, 2.18 Å, 9.7°); the third one is between the amide hydrogen of PHE 595 and sulphuryl oxygen atom (-N-H•••O=S, 2.00 Å, 53.7°).The angle data of H-bonds indicates that the In order to illustrate the interactions between B-Raf kinase and imidazopyridine, we focus on receptor-ligand interactions between B-Raf kinase and Mol 1 (the ligand of 4MBJ), Mol 2 (the most active inhibitor), and Mol 3 (the least active inhibitor).Figure 4a shows the docking mode between B-Raf kinase and Mol 1, in which four H-bonds were formed: the first one is between the >C=O of CYS 532 and the H-N1 (Figure 1) of Mol 1 (>C=O¨¨¨H-N1) with a distance of 1.89 Å and a deviated angle of 15.9 ˝; the second one is formed by amide hydrogen of TRP 531 and N4 of Mol 1 (-N-H¨¨¨N4, 2.20 Å, 18.2 ˝); the third one is between the carbonyl oxygen of GLY 593 and H-N3 of Mol 1 (>C=O¨¨¨H-N3, 2.13 Å, 45.0 ˝); and the fourth one is between the amide hydrogen of PHE 595 and sulphuryl oxygen atom of Mol 1 (-N-H¨¨¨O=S, 2.04 Å, 54.8 ˝).Due to their large angles (45.0 ˝and 54.8 ˝), the third and fourth H-bonds show less stability than the first and second H-bonds.As shown in Figure 4a, there are a π-π stacking contact between the imidazolepyridine ring of Mol 1 and the aromatic ring of PHE583, and a hydrophobic interaction between the benzene ring of Mol 1 and the side chain of THR529.The above observations are consistent with the previous studies [13].
The docking mode between B-Raf kinase and Mol 2 can be seen in Figure 4b, in which three H-bonds were formed: the first one is between the >C=O of CYS 532 and the H-N1 of Mol 2 (>C=O¨¨¨H-N1, 1.96 Å, 10.3 ˝); the second one is between the amide hydrogen of TRP 531 and N4 of Mol 2 (-N-H¨¨¨N4, 2.18 Å, 9.7 ˝); the third one is between the amide hydrogen of PHE 595 and sulphuryl oxygen atom (-N-H¨¨¨O=S, 2.00 Å, 53.7 ˝).The angle data of H-bonds indicates that the first two H-bonds are more stable than the third one.Similar with Mol 1, Mol 2 also shows a π-π stacking contact with aromatic ring of PHE583 and a hydrophobic interaction with the methyl group of THR529.Furthermore, the bromo-phenyl group attached the imidazole ring of Mol 2 interacts with the hydrophobic pocket formed by residues GLU533, Gly534, and SER535.
The docking mode between B-Raf kinase and Mol 3 is quite similar with the docking mode of Mol 2, which is shown in Figure 4c   As shown in Table S1, the pIC 50 values of Mol 1, Mol 2, and Mol 3 are 7.215, 9.119, and 6.777, respectively, which means the inhibitory activity: Mol 2 > Mol 1 > Mol 3.However, the docking C_score values of them are 6.37, 9.74, and 7.56, respectively, which indicates the docking effect: Mol 2 > Mol 3 > Mol 1.The docking result is not in accordance with the inhibitory activity completely.Therefore, in order to further explore interactions between the B-Raf kinase and imidazolepyridines, MD simulations were carried out in the subsequent work.

MD Simulations Features
Although docking analysis can provide an acceptable binding mode, the solvent effect and flexibility of protein were not fully taken into account.Therefore, MD simulations were carried out on the three docked complexes (Mol 1, Mol 2, and Mol 3 complex) to further explore the ligand-receptor interactions.
In order to evaluate the stability of the MD simulations, the properties of each complex (such as temperature, pressure, energy, and structure) were inspected during the entire MD trajectory.The fluctuations of temperature, pressure, and potential energy during the MD simulations are depicted in Figures S1-S3, respectively (supplementary data), which show that all of them are stable in the whole MD simulations process.The RMSD values of backbone atoms referring to the starting structure were used to monitor the dynamic stability of the MD trajectories.Figure 5 shows the RMSD for the Mol 1 complex, Mol 2 complex, and Mol 3 complex.For Mol 1 complex, the average RMSD fluctuations for the protein and ligand are 1.85 and 1.21 Å, respectively.The protein and Mol 1 reach to equilibrium after 4 ns.For Mol 2 complex, the protein and Mol 2 are quite stable after 6 ns, and the average RMSD fluctuations for protein and ligand are 2.12 and 0.99 Å, respectively.For Mol 3 complex, the protein, and Mol 3 reach to equilibrium after 8 ns, and the average RMSD fluctuations for the protein and ligand are up to 2.31 and 1.71 Å, respectively.The above results reveal the average RMSD fluctuations of the three ligands: Mol 2 < Mol 1 < Mol 3, which is in accordance with their inhibitory activity: Mol 2 > Mol 1 > Mol 3. 5

MD Simulations Features
Although docking analysis can provide an acceptable binding mode, the solvent effect and flexibility of protein were not fully taken into account.Therefore, MD simulations were carried out on the three docked complexes (Mol 1, Mol 2, and Mol 3 complex) to further explore the ligand-receptor interactions.
In order to evaluate the stability of the MD simulations, the properties of each complex (such as temperature, pressure, energy, and structure) were inspected during the entire MD trajectory.The fluctuations of temperature, pressure, and potential energy during the MD simulations are depicted in Figures S1-S3, respectively (supplementary data), which show that all of them are stable in the whole MD simulations process.The RMSD values of backbone atoms referring to the starting structure were used to monitor the dynamic stability of the MD trajectories.Figure 5 shows the RMSD for the Mol 1 complex, Mol 2 complex, and Mol 3 complex.For Mol 1 complex, the average RMSD fluctuations for the protein and ligand are 1.85 and 1.21 Å, respectively.The protein and Mol 1 reach to equilibrium after 4 ns.For Mol 2 complex, the protein and Mol 2 are quite stable after 6 ns, and the average RMSD fluctuations for protein and ligand are 2.12 and 0.99 Å, respectively.For Mol 3 complex, the protein, and Mol 3 reach to equilibrium after 8 ns, and the average RMSD fluctuations for the protein and ligand are up to 2.31 and 1.71 Å, respectively.The above results reveal the average RMSD fluctuations of the three ligands: Mol 2 < Mol 1 < Mol 3, which is in accordance with their inhibitory activity: Mol 2 > Mol 1 > Mol 3.

RMSF for Residues of the Binding Pocket
To explore the stability of the binding pocket during the MD simulations process, the root-mean-squared fluctuations (RMSF) of all the residues around the ligand at a ≤5 Å distance were calculated by the VMD software.Before the RMSF calculation, the average structures of the complexes were computed within the last 1 ns trajectory of MD simulations, and then each residue surrounding the ligand was aligned to the average structure.The residues around the ligand and their RMSF values compared with the starting structures are listed in Table 1.In all the complexes, the RMSF for each residue surrounding the ligand is lower than 1.0 Å, which means that the binding pocket is quite stable during the MD simulation.

RMSF for Residues of the Binding Pocket
To explore the stability of the binding pocket during the MD simulations process, the root-mean-squared fluctuations (RMSF) of all the residues around the ligand at a ď5 Å distance were calculated by the VMD software.Before the RMSF calculation, the average structures of the complexes were computed within the last 1 ns trajectory of MD simulations, and then each residue surrounding the ligand was aligned to the average structure.The residues around the ligand and their RMSF values compared with the starting structures are listed in Table 1.In all the complexes, the RMSF for each residue surrounding the ligand is lower than 1.0 Å, which means that the binding pocket is quite stable during the MD simulation.

H-Bonds in MD Simulations
H-bonds interaction is quite important in the binding between receptor and ligand.In our work, H-bonds were computed within the last 1 ns trajectory of MD simulations, and all the possible hydrogen acceptors were taken into consideration, such as ligands, protein, and water molecules.The distance cutoff was set to 3.00 Å (<3.00 Å), the angle cutoff (deviation from linearity) was set to 60.00 degree (<60.00 ˝), and the occupancy was set to 0.05% (>0.05%).The results of H-bonds analysis for the three systems in MD simulations are listed in Table 2.During the MD simulations, there were four hydrogen bonds formed in Mol 1 complex.The first H-bond is formed by the >C=O of CYS 532 and the H-N1 of Mol 1 with an occupation time of 54%, which is in accordance with the docking result.The second one is formed by the >C=O of THR 529 and the H-N2 of Mol 1 with an occupation time of 41%, the third one is formed by the >C=O of ASP 594 and the H-N3 of Mol 1 with an occupation time of 48%, and the fourth one is between oxygen of water (solvent) molecules and the H-N3 of Mol 1 with an occupation time of 36%.The second, third, and fourth H-bonds formed in MD are different from the docking result, which is caused by the solvent effect and movement of both receptor and ligand during the MD process.
Three hydrogen bonds were formed in the Mol 2 complex during the MD process.The first one is formed by the >C=O of CYS 532 and the H-N1 of Mol 2 (occupancy 66%), which is consistent with the docking result.The second one is formed by the >C=O of THR 529 and the H-N2 of Mol 2 (occupancy 11%), which indicates this H-bond is not important.The third one is formed by the >C=O of ASP 594 and the H-N3 of Mol 2 (occupancy 69%).The second and third H-bonds formed in MD show a difference with the docking result because of the movement of both receptor and ligand during the MD process.
Four hydrogen bonds were formed in Mol 3 complex during the MD process.The first one is formed by the >C=O of CYS 532 and the H-N1 of Mol 3 (occupancy 10%), which means this H-bond is not important.The second one is formed by the >C=O of ASP 594 and the H-N3 of Mol 2 (occupancy 97%), which means this H-bond is very important.The third and fourth H-bonds are between two amino hydrogen atoms of LYS 601 and the sulphuryl oxygen atom of Mol 3 (occupancy 42.5% and 26%).The H-bonds formed in MD are quite different from the docking result for Mol 3 complex, which indicates that there was a large movement of receptor and ligand during the MD process.It is in accordance with the result of RMSD fluctuations for Mol 3 complex.

Binding Free Energies
The calculated ∆G bind of the three complexes are shown in Table 3.It can be seen that the calculated ∆G bind values of the three complexes are consistent with their activities; their ∆G bind values being: ∆G bind (Mol 2) < ∆G bind (Mol 1) < ∆G bind (Mol 3), and their inhibitory activity pIC 50 values being: pIC 50 (Mol 2) > pIC 50 (Mol 1) > pIC 50 (Mol 3).The detailed contributions to the ∆G bind also can be obtained from Table 3.Firstly, the van der Waals interaction contribution (∆E vdw ) is the most important to the ∆G bind for each complex.For the Mol 2 complex, two benzene rings, and one imidazolepyridine ring show hydrophobic interaction with residues surrounding Mol 2, which is in accordance with the highest ∆E vdw (´59.02kcal¨mol ´1).For the Mol 3 complex, one benzene ring, one imidazole ring, and one piperidine ring show hydrophobic interaction with residues surrounding Mol 3, which is consistent with its middle ∆E vdw (´56.25 kcal¨mol ´1).For the Mol 1 complex, one benzene ring and one imidazole ring show hydrophobic interaction with residues surrounding Mol 1, which accords with the lowest ∆E vdw (´52.94kcal¨mol ´1).The result also indicates that the bromo-phenyl group in Mol 2 has stronger hydrophobic interaction with the residues GLU533, Gly534, and SER535 than piperidine group in Mol 3 does.Secondly, the electrostatic contribution (∆E ele ) is important to the ∆G bind for each complex as well, which is in accordance with the results of MD simulations because all three ligands formed no less than three H-bonds with the surrounding residues.Thirdly, the unfavorable polar solvation contribution (∆G GB ) affects the ∆G bind greatly, which shows the main differences for the three complexes.The reason why Mol 3 shows the lowest activity and lowest stability is that it has the highest unfavorable polar solvation contribution (∆G GB = 61.94),which is more than Mol 1 (∆G GB = 53.66)and Mol 2 (∆G GB = 56.45).Finally, the favorable non-polar solvation contributions (∆G SA ) to the ∆G bind for each complex are similar.

Preparation of Protein and Ligands
The SYBYL 7.3 software (Tripos Inc., St. Louis, MO, USA) package installed on Linux workstations was used to prepare the protein and ligands.Among the imidazopyridines, the crystal structure of B-Raf kinase combined with Mol 1 can be obtained from the protein data bank (PDB code: 4MBJ) [11], so B-Raf kinase receptor and Mol 1 were isolated from the complex.The protein extracted from the complex was treated by removing all of the substructures, removing all of the water molecules and adding hydrogen atoms.Without any conformation change, Mol 1 isolated from the complex and hydrogen atoms were added and geometrically optimized with three steps: (i) optimization using Steepest Descent with Gasteiger-Marsili charges and Tripos force field; (ii) optimization using conjugate gradient; and (iii) optimization using BFGS [14].The structures of all other imidazopyridines were built by modifying Mol 1, and geometrical optimizations were carried out by using the above procedure.

Molecular Docking
Molecular docking process between ligands and the receptor was carried out by using the Surflex-Dock module of SYBYL [15].In this program, a computational representation of the intended binding site (ProtoMol) was used to dock ligands into the binding site of a receptor automatically [16].In the present work, protomol_bloat was set to 0, protomol_threshold was set to 0.50 Å, and other parameters were set to default values.After docking, 10 conformations were present for each ligand, and the obtained final conformation was chose according to the following conditions: (i) the orientation of the docked conformation is in accordance with that of the ligand in crystal complex; and (ii) the conformation owns the highest C_score value.In the Surflex-Dock, the structures of ligands are flexible and the structure of the receptor is rigid.

MD Simulations
The AMBER 12 software package was used to carry out all the MD simulations [17].The initial structures of Mol 1, Mol 2, and Mol 3 complexes for the MD simulations were obtained from the docked results.The FF12SB AMBER force field was taken in the protein, and charges were added to the protein by using the software database.The general AMBER force field (GAFF) was taken for ligands [18], and AM1-BCC method was applied to assign their partial charges because of the lack of partial charge parameters for ligands in GAFF force field [19].The Antechamber suite in the AMBER 12 package was used to generate the topology files and atomic charges of ligands [20].The Tleap module of the AMBER 12 was used to produce the topology and coordinate files of the whole system.The whole system was dipped into a water box of TIP3P with a margin distance of 10 Å [21].In order to neutralize the charge of the system, a proper number of chloride ions were added.To deal with the long-range electrostatic interactions, the particle mesh Ewald (PME) was adopted during the MD simulations [22], and the cut-off distance of non-bonded interactions was set to 10 Å.The bonds involving hydrogen were constrained by the SHAKE algorithm [23].
Firstly, two stage energy minimizations were performed on each system: the algorithms (10,000 steps of the steepest descent and 10,000 steps of the conjugate gradient) with restrain were performed in the first stage; the same algorithms without restrain were further used in the second stage.Secondly, each system was heated from 0 to 300 K within 50 picoseconds (ps), gradually.Next, the system was equilibrated up to 500 ps at 300 K and constant pressure.Finally, a production process of 10 ns was performed in the constant temperature and pressure (NTP) with a step of 2 fs.The trajectories were recorded each 10 ps and the stability of the system was checked by the RMSD of the backbone.Trajectory analysis was carried out by using the CPPTRAJ [24].

Calculation of Binding Free Energy
The MM-GBSA method in AMBER 12 was used to compute the binding free energies (∆G bind ) of the receptor-ligand complexes [25].All the 100 snapshots of the simulated structures within the last 1 ns trajectory of MD simulations were extracted to perform the ∆G bind calculations.In MM-GBSA, ∆G bind is calculated as follows: ∆G bind " ∆G complex ´p∆G receptor `∆G ligand q (1) where ∆G complex , ∆G receptor and ∆G ligand are the free energy of the complex, receptor, and ligand, respectively.They can be obtained by the following equations: ∆G " ∆E gas `∆G sol ´T∆S gas (2) ∆E gas " ∆E ele `∆E vdw (3)
Int. J. Mol.Sci.2015, 16, page-page 3 C_score values, which are more than 5.0.The correlation between C_score values and biological activity (pIC50 values) of 36 imidazopyridines is shown in Figure 3.The above results indicate an acceptable reliability of the docking method for the B-Raf kinase receptor and these inhibitors.

Figure 2 .
Figure 2. Comparison between the docked and X-ray crystal structures of Mol 1 (red: docked structure; green: crystal structure).

Figure 2 .
Figure 2. Comparison between the docked and X-ray crystal structures of Mol 1 (red: docked structure; green: crystal structure).

Figure 2 .
Figure 2. Comparison between the docked and X-ray crystal structures of Mol 1 (red: docked structure; green: crystal structure).

Figure 3 .
Figure 3.The correlation between C_score values and pIC 50 values of 36 imidazopyridines.

Figure 4 .
Figure 4. Docking modes between B-Raf kinase and Mol 1 (a), Mol 2 (b) and Mol 3 (c).As shown in TableS1, the pIC50 values of Mol 1, Mol 2, and Mol 3 are 7.215, 9.119, and 6.777, respectively, which means the inhibitory activity: Mol 2 > Mol 1 > Mol 3.However, the docking C_score values of them are 6.37, 9.74, and 7.56, respectively, which indicates the docking effect: Mol 2 > Mol 3 > Mol 1.The docking result is not in accordance with the inhibitory activity completely.Therefore, in order to further explore interactions between the B-Raf kinase and imidazolepyridines, MD simulations were carried out in the subsequent work.

Table 1 .
Residues of the binding pocket and their RMSF values (Å).

Table 1 .
Residues of the binding pocket and their RMSF values (Å).

Table 2 .
H-bonds analysis for Mol 1, Mol 2, and Mol 3 in MD simulations.

Table 3 .
Binding free energy (kcal¨mol ´1) for the three complexes.

Activity Mol 1 Complex Mol 2 Complex Mol 3 Complex
∆E gas : molecular mechanics energy in gas phase; ∆E ele : electrostatic energy; ∆E vdw : van der Waals potential energy; ∆G sol : solvation free energy; ∆G GB : polar salvation free energy; ∆G SA : non-polar solvation free energy; ∆G bind : binding free energy; IC 50 : half maximal inhibitory concentration; pIC 50 : ´logIC 50 .