Insight into molecular dynamics simulation of BRAF(V600E) and potent novel inhibitors for malignant melanoma.

BRAF inhibitors have changed the standard therapeutic protocol for advanced or metastatic melanoma which harbored notorious BRAF(V600E) single mutation. However, drug resistance to BRAF inhibitors happens just like other cancer treatment. In this study, we constructed the ideal BRAF(V600E)-modeled structure through homology modeling and introduced the method of structure-based docking or virtual screening from the large compound database. Through certain methods of molecular dynamics simulation, we realized that BRAF(V600E) had quite prominent difference of molecular character or structural variation from the wild-type BRAF protein. It might confer the metamorphic character of advanced melanoma for the patients who harbored BRAF(V600E) mutation. By the methods of ligand-based quantitative structure-activity relationship and molecular dynamics simulation, we further recommend that aknadicine and 16beta-hydroxy-19s-vindolinine N-oxide from the traditional Chinese medicine are potent novel inhibitors for the management of malignant melanoma in the future.

Introduction

Drug use of BRAF inhibitors has become the exciting option of treatment for malignant melanoma patients who harbor B-RAF (BRAF) Val600Glu (V600E) single mutation.1 BRAF inhibitors have been one of the first choices of treating metastatic melanoma.2 More than half of the advanced melanomas contain BRAF gene mutation.3 Dozens of BRAF mutations have been verified, but BRAF(V600E) is the most notorious one. It replaces the normal amino acid valine (V) with glutamic acid (E) at position 600.4 Normal BRAF receives upstream signal from the growth factor receptor, receptor tyrosine kinase. BRAF regulates the mitogen/extracellular-activated protein kinase (MEK) and extracellular signal-regulated kinase (ERK) pathway. BRAF-MEK-ERK pathway controls the survival of the melanocytes in nevi or melanoma lesions.5 RRAF(V600E) mutation leads to uncontrolled cell growth in the melanoma patients.6 Traditional chemotherapy is insufficient for treating advanced or metastatic melanoma.7 Developing BRAF inhibitors is the milestone for the management of melanoma harboring BRAF(V600E) mutation.8 Vemurafenib was the first BRAF inhibitor approved by the US Food and Drug Administration in 2011.9 Vemurafenib (PLX4032) is the derivate from the experimental precursor, PLX4720.10 It improves the survival rates in advanced melanoma patients.11 Many clinical trials have proved amazing efficacy with vemurafenib treatment.12 Dabrafenib is the other approved BRAF inhibitor for advanced melanoma, which harbors BRAF(V600E) mutation.13 It has a dramatic effect in treating metastatic melanoma patients.14

However, the satisfaction of treatment just remains for a short period.15 Drug resistance to single BRAF inhibitor treatment occurs vigorously.16 The patients who were responsive to previous BRAF inhibitor suffer from disease progression.17 Many mechanisms explain why the resistance to BRAF inhibitors happens.18 There are mechanisms of primary and acquired resistance to RRAF inhibitor.19 Acquired resistance or compensatory reactivation of MEK/ERK signaling pathway is one reason for resistance.20 The reactivation of MEK/ERK pathway relies on aberrantly spliced BRAF(V600E) protein. BRAF(V600E) copy number amplification is sufficient for the resistance to BRAF inhibitor.21 MEK mutation has been implicated in BRAF inhibitor resistance, too.22 Upregulation and activation of the upstream receptor tyrosine kinase and expression of mutant N-RAS block the effect of BRAF inhibitor.23 NRAS and C-RAF cooperate to drive the resistance to BRAF inhibitor.24,25 Primary resistance includes loss or inactivation of essential tumor suppressors. Both the phosphatase and tensin homolog and retinoblastoma susceptibility gene are implicated in BRAF inhibitor resistance.26 Thus, we need combined strategy and alternative treatment to rescue advanced melanoma patients.27

BRAF protein has 766 amino acids and is composed of three main domains.28 The most important catalytic domain that phosphorylates consensus substrates is residues 457–717, conserved region 3. There are two lobes that are connected by a short hinge region. N-lobe (residues 457–530) is responsible for adenosine triphosphate (ATP) binding. C-lobe (residues 535–717) binds substrate proteins.29 Val471, Ala481, Leu514, Thr529, Trp531, and Cys532 form a pocket in which ATP is anchored. Asp(D)594, Phe(F)595, and Gly(G)596 compose a DFG motif, which determines the inactive or active state of BRAF protein. DFG-out or active state indicates that DFG motif moves out the pocket for ATP binding. DFG-out or active state indicates that DFG motif moves out of the pocket following ATP binding.30–32 ATP binding in N-lobe and substrate protein binding in C-lobe facilitate the phosphorylation of substrate protein. Thus, most BRAF inhibitors are designed to bind with the hinge region (residues 530–535) to prevent contact of ATP and substrate protein.33

In this study, we attempted to explore if there was any difference between wild-type BRAF and mutant BRAF(V600E) protein by computational simulation. Then we could discover the molecular character of this metamorphic protein. Investigation of the structural variation of BRAF(V600E) helped us understand the possible mechanism why the inhibitor resistance occurred in the conformational study.34 Systems biology is an indispensable science in modern biology and chemistry.35,36 Computer-aided drug design offers an alternative method in medicinal and biochemical technique of drug candidates searching. It has the advantage of screening out appropriate compounds for medicinal purpose efficiently. There are two major methods, structure-based and ligand-based, in computer-aided drug design. Structure-based method depends on three-dimensional (3D) structure and active sites of the target protein to survey the interaction, binding affinity, and steric relationship between the ligand and protein.37,38

As mentioned previously in the literature, searching BRAF inhibitors through structure-based virtual screening was a viable way to find out the potent compounds and their derivatives for further investigation from the large molecular database. In this study, we tried to identify if there was any compound better than the existing or experimental BRAF inhibitors by virtual screening from the large traditional Chinese medicine (TCM) compound database.39–42 Ligand-based methods depend on the features of the ligand when it interacts with the target protein. This technique adopts statistical algorithms to link the conformational features and biological activities when a ligand interacts with the protein. The other novelty and significance of this study is molecular dynamics (MD) simulation. It is one of the most powerful tools in the biological study. The interaction of a ligand and protein is a dynamic process. Because molecular interaction contains a large number of atoms, it is necessary to analyze the enormous data of atom deviation or fluctuation by MD simulation. It provides detailed information about dynamic changes of the ligand and conformational changes of the protein. Through MD simulation, we could explore the mystery of wild-type BRAF from BRAF(V600E) protein. In addition, we could examine if the selected potent inhibitors bound with BRAF(V600E) protein successfully and stably.

Materials and methods

Homology modeling

We obtained the sequence of human BRAF protein (residues 1–766) from the UniProt Knowledgebase (P15056, human) and 3D conformation (residues 445–723) from Protein Data Bank (ID: 2FB8, human). Homology modeling of BRAF(V600E) with crystal structure of wild-type BRAF (2FB8) was constructed by the Build Homology Models program of Accelrys Discovery Studio (DS). We further confirmed the BRAF(V600E)-modeled structure by a Ramachandran plot with Rampage program and by Verify score with Profiles-3D program in DS.43

Structure-based docking or virtual screening

We prepared structure-based docking procedure for SB-590885 and made all the small molecular compounds from TCM Database@Taiwan to dock with BRAF(V600E) protein. LigandFit module in DS was conducted to perform the docking procedure. Key position of the binding sites was set at the residues 530–535. We minimized all docking poses between the ligand and BRAF(V600E) protein by the force field of Chemistry at HARvard Molecular Mechanics. LIGPLOT module in DS was adopted to illustrate hydrogen bond (H-bond) or other binding forces between the ligand and BRAF(V600E) protein.44

Ligand-based quantitative structure–activity relationship (QSAR) models

We obtained 36 compounds and the pIC50 data of BRAF inhibitors to draw two-dimensional and 3D structures of these compounds by ChemBioDraw software.45 Calculate Molecular Property program and Genetic Function Approximation (GFA) algorithm in DS were conducted to find the appropriate molecular descriptors. These descriptors were verified by libSVM and Matlab Statistics Toolbox to construct the multiple linear regression (MLR) and support vector machine (SVM) models. We conducted Bayes Net Toolbox (BNT) for Bayesian network modeling to predict the pIC50 value of any given compound. These predicted models adopted fivefold cross-validation. We calculated the activity of square correlation coefficient (R2) by the pros and cons of each model. The highest R2 of SVM, MLR, and Bayesian network, which were chosen to be the predictive models, had sufficient capacity to predict the activity of selected candidates and the control (SB-590885).46

MD simulation

We illustrated the trajectories of MD simulation by the program of GROningen MAchine for Chemical Simulations. Each ligand, protein as well as complex passed through four phases: minimization, heating, equilibration, and production. We drew the trajectories of root mean square deviation (RMSD), mean square deviation (MSD), gyration radius (Gyrate), and solvent-accessible surface area (SASA) for the ligand, protein, and complex. Ligands included the first two selected candidates (top1,2) and the control (SB-590885). Protein included BRAF alone (apo), BRAF(V600E) alone (mutants-apo), the corresponding BRAF(V600E) of the candidates, and the control. Complex included the candidate–BRAF(V600E) complex and the control–BRAF(V600E) complex. Total energy, root mean square fluctuation, principal component analysis, cluster analysis, and occupancy of H-bond were also conducted in this study.47

Results and discussion

We substituted the amino acid valine (V) with glutamic acid (E) at position 600 in the sequence of BRAF protein and constructed the BRAF(V600E)-modeled structure based on the wild-type BRAF template (2FB8, human). Because the most significant catalytic domains that phosphorylated consensus substrates were residues 457–717, the conformation of 445–723 amino acids was sufficient to represent BRAF and BRAF(V600E)-modeled structure. Ramachandran plot of BRAF(V600E)-modeled structure displayed that 81.8% of residues were in the favored region, 14.6% were in the allowed region, and only 3.6% were in the outlier region (Figure 1). Verify score of wild-type BRAF (2FB8) and BRAF(V600E) protein showed that almost all the residues had positive values (Figure 2).

Figure 1

Ramachandran plot of homology modeling for BRAF(V600E)-modeled structure.

Notes: Number of residues in favored region (~98.0% expected): 207 (81.8%). Number of residues in allowed region (~2.0% expected): 37 (14.6%). Number of residues in outlier region: 9 (3.6%).

Figure 2

Verify score of wild-type BRAF (2FB8) and BRAF(V600E) protein (mutant).

Note: Positive scores mean accurate structure of amino acids.

To investigate the conformational character of BRAF (V600E) protein, we needed to construct the ideal BRAF (V600E)-modeled structure and confirm it through rigorous verification. The high percentage in the favored (81.8%) and allowed (14.6%) regions calculated by Ramachandran plot and qualified values presented by Verify score indicated that the BRAF(V600E)-modeled structure was an ideal conformation model.

We assigned the experimental BRAF inhibitor, SB-590885, as the ligand-binding control compound for BRAF(V600E) protein. SB-590885 was also the compound that existed in the crystal complex of 2FB8 recruited by the Protein Data Bank. We listed the first ten compounds screened by virtual screening from the world’s largest TCM Database@Taiwan. Dock score of the ten TCM compounds was better than that of SB-590885 (Table 1). According to the result of Dock score, we selected aknadicine (top1) and 16beta-hydroxy-19s-vindolinine N-oxide (top2) as the candidates for docking analysis, QSAR, and MD validation. We illustrated the binding residues of BRAF(V600E) protein with aknadicine, 16beta-hydroxy-19s-vindolinine N-oxide, and SB-590885. Aknadicine formed H-bond with Gly534 and Ser535 of BRAF(V600E) protein. 16beta-hydroxy-19s-vindolinine N-oxide formed H-bond with Trp531 and Glu533 of BRAF(V600E) protein. The control formed H-bond with Lys473, Trp531, and Glu533 BRAF(V600E) protein (Figure 3).

Table 1

Dock score and quantitative structure–activity relationship values

Name	Dock score	Predicted activity
		Support vector machine	Multiple linear regression	Bayes Net Toolbox
Aknadicine	68.69	6.83	14.92	7.13
16beta-Hydroxy-19s-vindolinine N-oxide	59.16	6.84	14.51	17.26
beta-Peltatin	55.49	6.92	13.17	9.05
clemaphenol_A	54.34	6.85	12.10	9.10
gomisin_J	54.19	6.95	13.01	9.41
Schisanlignone B	52.21	6.89	12.38	9.02
beta-Peltatin	51.19	6.89	12.06	8.50
Ribalinium	46.63	6.74	12.96	8.05
Gomisin J	47.03	6.96	13.27	9.52
Yadanzioside N	44.97	6.84	13.32	8.78
SB-590885*	31.51	6.79	7.87	12.16

Figure 3

Scaffold (left) and binding sites (right): (A) aknadicine, (B) 16beta-hydroxy-19s-vindolinine N-oxide, (C) the control, SB-590885.

Note: Key position of the binding sites was set at the residues 530–535.

Because most BRAF inhibitors are designed to bind with the hinge region (residues 530–535) to prevent contact of ATP and substrate protein, we set the binding sites at the position 530–535 of BRAF(V600E) protein. Both the candidates and the control could bind to the region with the most important binding force, H-bond. Dock score of aknadicine, 16beta-hydroxy-19s-vindolinine N-oxide, and the control was 68.69, 59.16, and 31.51, respectively. This result demonstrated that more than ten TCM compounds could dock with BRAF(V600E) protein better than the control; thus, we selected the first two potent inhibitors, aknadicine and 16beta-hydroxy-19s-vindolinine N-oxide, as the candidates for further QSAR and MD validation.

Ligand-based QSAR models

We acquired 36 compounds and pIC50 data of BRAF inhibitors. The 36 compounds of BRAF inhibitors were randomly divided into 27 training sets and nine test sets for validation. We chose the following seven optimum descriptors for constructing GFA template model: ALogP, ES_Count_sssN, ES_Count_ssssC, LogD, Num_AromaticBonds, Minimized_Energy, and Shadow_XZfrac.

After determining GFA template model, we constructed SVM-, MLR-, and BNT-predictive models. Then we obtained the following mathematical slope equation by illustrating observed and predicted activity of these compounds (Table 2): R2 value of SVM, MLR, and BNT was 0.8576, 0.8449, and 0.8293, respectively. The high R2 value meant that the predictive activity of any chosen compound was almost the same with its observed activity (Figure 4). Predicted activity of SVM value for aknadicine, 16beta-hydroxy-19s-vindolinine N-oxide, and the control was 6.83, 6.84, and 6.79, respectively. MLR value was 14.92, 14.51, and 7.87, respectively. BNT value was 7.13, 17.26, and 12.16, respectively (Table 1). The results of these predictive values demonstrated that the two candidates, aknadicine and 16beta-hydroxy-19s- vindolinine N-oxide, might have better biological activities than the control in general.

Table 2

Observed (pIC50) and predicted activity of 36 experimental BRAF(V600E) inhibitors

Name	pIC50	Support vector machine predict	Multiple linear regression predict	BNT predict
1	7.63827	7.41024	7.6194769	7.374755
21a	9	8.49954	7.947385429	7.570344
21c	7.63827	7.94884	7.715589335	7.507701
21d	8.09691	7.63028	7.69481958	7.599865
21f	6.97469	7.47464	7.689070605	7.392667
21g	5	5.50016	5.404198765	7.451606
21j	5.18709	5.68728	5.372310431	5.18709
21k	7.3279	7.82789	7.757460256	7.443077
21m	7.4318	7.93152	7.66824277	7.444933
21n	9	8.49142	8.031087118	7.626747
21o	7.92082	7.95511	7.714453524	7.51947
21p	7.79588	8.58137	8.08604531	7.648566
21r	8.1549	8.33563	7.866688002	8.210212
21s	8.1549	8.65499	8.317593039	8.16193
21t	8.22185	8.41459	8.145699383	6.733601
21u	7.4318	7.93179	7.758683368	7.026349
21w	7.10237	7.44954	7.10237	7.760075
21x	7.85387	8.27306	8.142843629	6.19518
21y	6.19518	6.36307	5.893678679	7.639428
22b	7.63827	7.39801	8.208587719	7.61749
22c	7.76955	7.26955	7.480139135	7.23233
22e	8	7.50047	7.588321338	6.760189
22f	7.61979	7.11976	7.742020164	6.553714
23a	6.80688	6.86424	6.952963588	6.635629
23c	6.36653	6.86646	6.772286613	5.284
25a	5.284	5.36044	4.734245058	5.33724
25b	5.33724	5.83723	5.542500261	7.022547
21b	7.60206	7.55399	7.675530551	7.580189
21e	7.52288	8.02369	7.800943114	7.580609
21h	7.25181	7.91268	7.531800304	5.929145
21i	7.65758	8.00893	7.712916529	7.424633
21l	7.67778	7.96514	7.624559479	7.419267
21q	7.92082	7.98343	7.73585786	7.513427
21v	6.65758	7.15765	6.524388082	3.495288
23b	6.77728	7.16443	6.736866666	7.56727
26	7.04576	5.75116	6.973415379	10.43935

Figure 4

Ligand-based quantitative structure-activity relationship models. (A) SVM: support vector machine; (B) MLR: multiple linear regression; (C) BNT: Bayes Net Toolbox.

Note: Green lines mean 95% confidence interval.

Molecular character

First, we compared two differences of molecular character in the process of MD. One was wild-type BRAF and BRAF(V600E) protein, the other was BRAF(V600E) protein to which both the candidates (aknadicine and 16beta-hydroxy-19s-vindolinine N-oxide) and the control bound.

We drew the trajectory of RMSD to investigate the deviation degree of each ligand (ligand RMSD) and the ligand–BRAF(V600E) protein complex (complex RMSD). The average ligand or complex RMSD value of the two candidates was smaller than that of the control (Figure 5). These results demonstrated that either aknadicine or 16beta-hydroxy-19s-vindolinine N-oxide had smaller ligand and ligand–protein complex deviation than the control. It was evident that both the candidates could bind to BRAF(V600E) protein stably.

Figure 5

Ligand (A) and complex (B) root mean square deviation (RMSD).

Note: Aknadicine and 16beta-hydroxy-19s-vindolinine N-oxide have smaller ligand deviation than the control.

The trajectory of MSD was drawn to observe the deviation degree of wild-type BRAF protein alone (apo), BRAF(V600E) protein alone (mutants-apo), and the corresponding BRAF(V600E) protein of each ligand. Mutants-apo had larger average protein MSD value than apo. Aknadicine corresponding protein had smaller average protein MSD than the control, but the 16beta-hydroxy-19s-vindolinine N-oxide corresponding protein had larger average value than the control (Figure 6A). This result demonstrated that BRAF and BRAF(V600E) protein alone had different molecular characters. It was evident that although when 16beta-hydroxy-1dolinine N-oxide bound to BRAF(V600E) protein, the larger protein deviation did not affect the stability of the ligand and its complex, which had been described in ligand and complex RMSD.

Figure 6

(A) Protein mean square deviation (MSD) and (B) protein gyration radius (Gyrate).

Note: BRAF alone (apo) and BRAF(V600E) protein alone (mutants-apo) have different molecular character.

We drew the trajectory of Gyrate to evaluate the compact degree of apo, mutants-apo, and the corresponding protein of each ligand. Apo had larger average Gyrate value than mutants-apo. The average protein Gyrate of the two candidates was smaller than that of the control (Figure 6B). This result demonstrated that BRAF and BRAF(V600E) protein alone had different character again. It was evident that even single mutation of BRAF protein led to quite different molecular character.

The trajectory of SASA was drawn to show water contact surface of the ligand and its corresponding protein. Apo had larger average protein SASA than mutants-apo. The average ligand SASA value of the two candidates was smaller than that of the control (Figure 7). These results demonstrated that there was another different molecular character between BRAF and BRAF(V600E) protein alone. It was evident that both the candidates had smaller volume or less hydrophobic side chains.

Figure 7

Protein (A) and ligand (B) solvent-accessible surface area (SASA).

Note: Aknadicine and 16beta-hydroxy-19s-vindolinine N-oxide have smaller volume or less hydrophobic side chains than the control.

We illustrated the graph of total energy to assess the most stable energy needed for apo, mutants-apo, either the candidates or the control when they bound to BRAF(V600E) protein. The most stable total energy for apo or mutants-apo was almost the same. Total energy needed for the control bound with BRAF(V600E) protein was prominently higher than both the candidates (Figure 8). These results demonstrated that BRAF(V600E) protein alone did not change its stability compared to wild-type BRAF protein alone. It was evident that both the candidates could bind to BRAF(V600E) protein more easily than the control.

Figure 8

Total energy. SB-590885 had higher peak total energy than aknadicine and 16beta-hydroxy-19s-vindolinine N-oxide.

Structural variation

Second, we compared the detailed structural or conformational variation for wild-type BRAF, BRAF(V600E) protein, both the candidates, and the control when they bound to BRAF(V600E) protein.

The number matrix of root mean square fluctuation was illustrated to identify the structural similarity from the view of individual residue. The number “1” indicated the highest similarity when we compared two proteins. The relative value between BRAF(V600E) alone (M-apo) and wild-type BRAF (apo) was 0.1225 (Table 3). This result demonstrated that the structural variation of BRAF(V600E) protein was quite prominent.

Table 3

Number matrix of root mean square fluctuation (RMSF) for wild-type BRAF alone (apo), BRAF(V600E) alone (M-apo), and the corresponding BRAF(V600E) protein for aknadicine (top1), 16beta-hydroxy-19s-vindolinine N-oxide (top2), and the control (SB-590885)

	Apo	M-apo	Top1	Top2	Control
Apo	1	0.1225	0.1158	0.13	0.1113
M-apo	0.1225	1	0.8192	0.7798	0.7117
Top1	0.1158	0.8192	1	0.8021	0.7049
Top2	0.13	0.7798	0.8021	1	0.7442
Control	0.1113	0.7117	0.7049	0.7442	1

Note: The number “1” means the highest similarity when we compare two proteins.

Abbreviation: M-apo, mutated apo.

We performed principal component analysis to discuss the similarity of two eigenvectors (PC1 and PC2) based on the backbone of BRAF, BRAF(V600E) protein alone, and the corresponding BRAF(V600E) protein of each ligand. Apo and mutants-apo had certain degree of difference in the eigenvectors. Comparison of the candidates or the control with mutants-apo, the eigenvectors only had little difference. There still was prominent difference in the eigenvectors of apo and either the candidates or the control (Figure 9). These results demonstrated that the structural variation of BRAF(V600E) was very meaningful and may be the reason why the patient had poor prognosis of advanced melanoma harboring just single mutation.

Figure 9

Principal component analysis.

Notes: Apo (blue color) and (A) mutants-apo (red color), (B) top1 (red color), (C) top2 (red color), (D) the control (red color). Mutants-apo (blue color) and (E) top1 (red color), (F) top2 (red color), (G) the control (red color). The control (blue color) and (H) top1 (red color), (I) top2 (red color).

Cluster analysis was displayed to show the representative structure of BRAF (apo), BRAF(V600E) (mutants-apo), and the corresponding protein for aknadicine, 16beta-hydroxy-19s-vindolinine N-oxide or the control. We chose 12.14 ns of apo and 17.04 ns of mutants-apo for calculating the diameter of the pocket in which ATP was anchored. The representative structure after 15 ns was chosen for analyzing the binding position between the ligand and its corresponding protein. We chose 18.9, 15.92, and 16.52 ns of aknadicine, 16beta-hydroxy-19s-vindolinine N-oxide, and the control to illustrate the schematic diagram, respectively (Figure 10).

Figure 10

Cluster analysis for choosing the representative structure.

As mentioned in the introduction, Val471, Ala481, Leu514, Thr529, Trp531, and Cys532 formed a pocket in which ATP was anchored. Asp(D)594, Phe(F)595, and Gly(G)596 composed a DFG motif, which determined the inactive or active state of BRAF protein. We illustrated the schematic diagram to explore the diameter change of the pocket in which ATP was anchored in the DFG-out or active state. The diameter of apo at 0 ns and 12.14 ns was 12.818 Å and 15.404 Å, respectively. The diameter of mutants-apo at 0 ns and 17.04 ns was 12.536 Å and 10.518 Å, respectively (Figure 11). This result demonstrated that the ATP pocket of apo became larger during MD, but the pocket of mutants-apo became smaller in contrast. We speculated that the diameter change of ATP pocket amplified the phosphorylation and the metamorphic character of BRAF(V600E) protein.

Figure 11

Schematic diagram of the diameter change of ATP pocket for (A) wild-type BRAF at 0 ns and 12.14 ns or (B) BRAF(V600E) protein at 0 ns and 17.04 ns.

Note: The main chain of Asp(D)594, Phe(F)595, and Gly(G)596 (DFG-motif) was painted as green color.

Abbreviation: ATP, adenosine triphosphate.

Occupancy of H-bond between the candidates and the control with BRAF(V600E) protein was shown to explain the dynamic change of H-bond during the period of MD. The control formed more H-bonds with BRAF(V600E) protein than the candidates (Table 4). However, the schematic diagram of binding position between the ligand and its corresponding protein showed that the control had larger deviation than the candidates (Figure 12). This result was consistent with ligand RMSD and demonstrated that H-bonds of the candidates were less than the control, but they were relatively stable when they bound to BRAF(V600E) protein. We speculated that the larger deviation conferred drug resistance to the control or other existed BRAF inhibitors.

Table 4

Occupancy of H-bond between aknadicine, 16beta-hydroxy-19s-vindolinine N-oxide, and SB-590885 with BRAF(V600E) protein

Aknadicine		16beta-hydroxy-19s-vindolinine N-oxide		SB-590885
H-bond interaction	Occupancy	H-bond interaction	Occupancy	H-bond interaction	Occupancy
HD22:ASN580/O15	1.50%	HE2:HIS539/N13	32.90%	HN:GLU586/N21	45.10%
HD22:ASN580/O19	4.80%	HE2:HIS539/O29	1.60%	HE2:HIS585/N21	4.20%
HE2:HIS539/O15	0.50%	HG1:SER535/N13	3.70%	HE2:HIS539/N21	2.00%
HG1:SER536/O15	32.20%	HG1:SER535/O29	3.60%	HE2:HIS539/O8	27.90%
HG1:SER536/O19	22.40%	HN:GLY534/O30	0.80%	HE2:HIS539/N9	29.90%
HN:SER536/O15	6.30%	HN:SER465/O25	0.60%	HN:SER536/O30	1.70%
HG1:SER535/O15	43.90%	HH22:ARG462/O25	0.10%	HE1:TRP531/O8	1.50%
HG1:SER535/O16	9.20%	HE22:GLN461/O25	2.40%	HE1:TRP531/N9	0.90%
HG1:SER535/O17	4.50%	HE22:GLN461/O30	0.50%	HE1:TRP531/O30	0.80%
HG1:SER535/O19	39.00%	OE1:GLN461/H52	0.50%	HZ3:LYS473/O8	5.30%
HE1:TRP531/O15	0.10%	OE1:GLU533/H52	5.30%	HZ3:LYS473/N9	2.00%
HE22:GLN461/O15	1.20%	OE2:GLU533/H52	7.40%	HZ3:LYS473/N32	1.80%
HE22:GLN461/O19	1.10%	O:ILE463/H51	3.80%	HN:SER465/O30	1.40%
O:ILE463/H49	22.10%	ND1:HIS539/H51	18.80%	HN:SER465/N32	1.20%
O:SER465/H49	2.70%	OE1:GLN461/H22	0.50%	HE22:GLN461/N9	6.70%
OG:SER535/H49	5.70%	O:ARG462/H22	2.40%	HE22:GLN461/N32	0.40%
O:CYS532/H34	2.50%	O:ILE463/H22	14.80%	O:GLU533/H62	53.40%
				O:GLY534/H62	6.70%
				OG:SER535/H62	25.40%
				O:GLY534/H61	8.80%
				OE1:GLN461/H60	6.50%
				O:ARG462/H60	0.40%
				OE2:GLU533/H60	0.60%

Figure 12

Schematic diagram of binding position between the ligand and its corresponding BRAF(V600E) protein. (A) Aknadicine, (B) 16beta-hydroxy-19s-vindolinine N-oxide, and (C) the control.

Note: The two candidates are relatively stable than the control.

Conclusion

Even BRAF(V600E) protein harbored just a single mutation, it had quiet different molecular character from the wild-type BRAF protein as determined by MD simulation, such as MSD, Gyrate, and SASA analysis. As regards the individual residue, BRAF(V600E) had a quite prominent structural or conformational variation than the wild-type BRAF protein by the matrix of root mean square fluctuation and principal component analysis. On the basis of structure-based virtual screening, ligand-based QSAR models, and MD simulation, we recommend aknadicine and 16beta-hydroxy-19s-vindolinine N-oxide to be potent compounds for developing novel inhibitors in the future.