Introduction

The mitogen activated protein kinase (MAPK) signaling cascades regulate cell growth, proliferation, and survival in mammalian cells (Samatar and Poulikakos, 2014; Lavoie et al., 2020). In the well-studied Ras-Raf-MEK-ERK pathway, the GTP-loaded RAS contacts RAF and induces its dimerization; the newly formed RAF dimer phosphorylates MEK which in turn phosphorylates ERK, which then phosphorylates a number of downstream proteins and regulate their functions (Lavoie et al., 2020). Mutations of BRAF, a kinase within the RAF family, are present in about 8% of human tumors, most commonly melanoma and colorectal cancers, with the mutation V600E accounting for about 90% of them. It is believed that the wild type BRAF signals as a dimer, the BRAF^V600E is able to signal as a monomer (Poulikakos et al., 2010, 2011; Karoulia et al., 2017). The first generation BRAF^V600E inhibitors, including the current FDA-approved inhibitors, Vermurafenib, Dabrafenib, and Enco-rafenib, inhibit the monomeric BRAF^V600E; however, drug resistance led to only short-term cancer remission in patients (Poulikakos et al., 2011; Peng et al., 2015; Monaco et al., 2021; Yen et al., 2021; Adamopoulos et al., 2021). In the adaptive drug resistance mechanism, RAF dimerization renders the monomer-selective inhibitors ineffective (Poulikakos et al., 2011; Peng et al., 2015; Monaco et al., 2021; Yen et al., 2021; Adamopoulos et al., 2021). To overcome the resistance, inhibitors that are either dimer selective or equipotent to both monomers and dimers have been developed and entered clinical development (Adamopoulos et al., 2021; Cook and Cook, 2021). Understanding the molecular mechanism of dimer selectivity would be valuable for the rational design of RAF inhibitors. Although several MD studies have examined the conformational dynamics of BRAF^V600E monomer (Maloney et al., 2021), BRAF^V600E in complex with the monomer-selective inhibitors (Tse and Verkhivker, 2016), and wild type RAF dimerization (Zhang et al., 2021), the topic of RAF dimer selectivity has not been explored.

The kinase domain of the BRAF monomer has a typical kinase structure: a primarily β-sheet N-terminal domain connected to a helical C-terminal domain by a flexible hinge (Figure 1). Like other kinases, the catalytic activity of BRAF depends on the conformation of two motifs: the αC-helix, which contains the conserved residue Glu501, and the DFG motif on the activation loop (a-loop), which contains the conserved ATP-binding (via magnesium) residue Asp594. In the active state, both the αC helix and DFG adopt the IN conformation, dubbed CIDI. In this state, the αC helix is positioned inward such that αC-Glu501 and the catalytic Lys483 form a salt bridge; meanwhile the DFG motif is also IN, meaning DFG-Asp594 is near Lys483 often in a salt-bridge distance. An inactive conformation can be achieved if either or both the αC helix and DFG motif adopt an OUT state. Specifically, αC-out involves an outward movement of the αChelix, while DFG-out involves the sidechains of the DFG Asp594 and Phe595 exchanging regions, i.e., Phe595 facing the ATP binding site and Asp594 facing the αC-helix. In the BRAF dimer, the two protomers are arranged side by side and the dimer interface involves the C-terminal end of the αC helix (Figure 1). Current monomer-selective BRAF^V600E inhibitors bind in the αC-out conformation, whereas the dimer-selective or equipotent inhibitors bind in the αC-in conformation (Supplemental Table 1). Thus, the αC conformation has been the center of attention in numerous structural and biochemical studies to understand RAF signaling and inhibitor activities (Rajakulendran et al., 2009; Thevakumaran et al., 2015; Karoulia et al., 2016).

Table 1.

Figure 1.

In a recent study, Gavathiotis and coworkers discovered a modification to the dimer-compatible inhibitor Ponatinib which can increase the dimer selectivity by more than three fold (Cotto-Rios et al., 2020). The novel inhibitor, named Ponatinib hybrid inhibitor 1 (PHI1), extends the head-group of Ponatinib by replacing the methylpiperazine with the 4-(2-aminoethyl) morpholino group. Remarkably, PHI1 showed more potent inhibition of the second protomer in the BRAF^V600E dimer; in contrast, Ponatinib and equipotent inhibitors, e.g., LY3009120 or LY, AZ-628, and TAK-632, are non-cooperative (Cotto-Rios et al., 2020). The co-crystal structure of BRAF^V600E in complex with PHI1 (PDB: 6P7G) (Cotto-Rios et al., 2020) revealed that the morpholine group extends the ligand-kinase interaction from the type-II pocket (occupied by all DFG-out inhibitors) to the center of αC helix, allowing a hydrophobic interaction with Asn500 next to the αC-Glu501 (Figure 1). The co-crystal structures show that this interaction is not available with the shorter Ponatinib (PDB ID: 6P3D) (Cotto-Rios et al., 2020) or equipotent inhibitors, e.g., LY3009120 (LY, PDB ID: 5C9C) (Peng et al., 2015). Gavathiothis and coworkers noticed that PHI1 stabilizes the αC helix in a slightly different IN conformation as compared to Ponatinib and hypothesized that the additional interaction with Asn500 is a key to the dimer selectivity of PHI1, as it may be unfavorable in monomer binding (Cotto-Rios et al., 2020). Shortly after, a biochemical study supported by the molecular dynamics (MD) simulations suggested that restriction of the αC helix movement is the basis for the difference between dimer-selective and equipotent inhibitors (Adamopoulos et al., 2021); however, the detailed mechanism remains elusive.

Prompted by the open questions regarding dimer selectivity and binding cooperativity of BRAF^V600E inhibitors, we carried out a series of all-atom molecular dynamics (MD) simulations to investigate the conformational dynamics of the monomeric and dimeric BRAF^V600E in the presence and absence of one or two dimer-selective (PHI1) or equipotent (LY) inhibitor(s). Analysis of the simulation data which ammounts to 135 μs aggregate time uncovered the atomic details of the remarkable conformational allostery in BRAF^V600E dimerization and inhibitor binding. Supported by the co-crystal structure analysis of the published monomer-selective, dimer-selective, and equipotent inhibitors, an atomically-detailed mechanism emerged that explains the monomer or dimer selectivity and binding cooperativity of BRAF^V600E inhibitors. The mechanism also led us to propose an empirical method based on the co-crystal structure for assessing the dimer selectivity of BRAF^V600E inhibitors.

Results and Discussion

Analysis of the co-crystal structures suggests the h-bond formation with αC-Glu501 as a key requirement for dimer binding

To understand the preference of BRAF^V600E inhibitors for the monomer vs. dimer form, we first examined all published co-crystal structures in complex with the monomer-selective and dimer-compatible (i.e., dimer-selective and equipotent) inhibitors (see Supplementary table 1 for a complete list). We first noticed that the monomer-selective inhibitors, e.g., Vemurafenib (VEM, PDB ID: 5JRQ) (Grasso et al., 2016), do not occupy BP-III, whereas most dimer-compatible inhibitors do. This can be explained by the observation that the monomer-selective inhibitors bind in the DFG-in, whereas most dimer-compatible inhibitors bind in the DFG-out conformation–BP-III is occupied by Phe595 in the DFG-in conformation, so the pocket is only available in the DFG-out conformation (Figure 2a). Note, the equipotent inhibitor SB590885 (PDB ID: 2FB8) (King et al., 2006) does not occupy BP-II or BP-III, as it binds in the DFG-in conformation (Supplementary table 1).

Figure 2.

The co-crystal structure analysis revealed an important distinction between the monomer-selective and dimer-compatible inhibitors, namely, the former binds in the αC-out whereas the latter binds in the αC-in conformation. The interaction fingerprints showed that while most monomer-selective inhibitors make a hydrophobic contact with Leu505 next to the conserved RKTR motif at the end of the αC helix, only the dimer-compatible inhibitors interact with αC-Glu501 by donating a h-bond (e.g., from an amide group in PHI1 and LY) to the carboxylate sidechain of Glu501 (Figure 2b). Glu501 rests above BP-II in the DFG-out conformation (called BP-II-out) and may interact with the catalytic Lys483 (see later discussion), which makes up a part of BP-I. Interestingly, even though SB590885 binds in the DFG-in conformation, it can also donate a h-bond to Glu501 through an oxime hydroxyl group (PDB ID: 2FB8) (King et al., 2006). This h-bond stabilizes the salt-bridge between the catalytic Lys483 and Glu501 such that the αC helix position is further inward (according to the KLIFs definition, see later discussion) as compared to the co-crystal structures in complex with other dimer-compatible inhibitors (Supplementary Table 1).

All monomer-selective and dimer-compatible inhibitors interact with the DFG-Asp594 although with subtle differences. In PHI1 (PDB ID: 6P7G) and LY (PDB ID: 5C9C), the amide carbonyl occupying the BP-II accepts a h-bond from the backbone amide of Asp594, while in VEM (PDB ID: 4RZV) (Grasso et al., 2016) the sulfonamide occupying the BP-II donates a h-bond to the backbone amide of Asp594 (orange in Figure 2a). One unique property of PHI1 is the ability to donate a h-bond to the backbone carbonyl of HRD-His574 through an amino nitrogen next to the morpholine head-group. This region is classified as BP-IV by KLIFS, although the sidechain of His574 is a part of BP-II (in PDB 6P7G) and makes a h-bond with the backbone of the DFG-1 Gly593. Among the other dimer-compatible inhibitors, only Ponatinib (PDB ID: 6P3D) (Cotto-Rios et al., 2020; Adamopoulos et al., 2021) makes a similar h-bond with the backbone of His574 through the methyl pyrazine headgroup.

In addition to analyzing the co-crystal structures, we also tested the inhibition of ERK1/2 phos-phorylation in two melanoma cell lines by PHI1, LY, or VEM (Figure 2c,d). SKEML239 expresses monomeric BRAF^V600E, while SKMEL239-C4 expresses dimeric BRAF^V600E(Cotto-Rios et al., 2020). Among the three inhibitors tested, PHI1 is the only compound to be more potent against SKMEL239-C4 versus SKMEL239 (IC50 of 256 nM vs. 1.5 μM). By contrast, LY has similar potency (27 nM vs. 15 nM) while VEM is more potent against SKMEL239 (3 μM vs 35 nM). This data confirms that PHI1 is dimer-selective, LY equipotent, and VEM monomer-selective.

In light of the above finding and given the central location of Glu501 on the αC helix, we hypothesized that the ability to form a h-bond with Glu501 is required by dimer-compatible inhibitors, as the h-bonding would restrict the αC helix to the αC-in conformation as observed in the co-crystal structures of all dimer-compatible inhibitors. This restriction was also suggested as a key for dimer selectivity in the recent study by Poulikakos and coworkers (Adamopoulos et al., 2021). However, the crystal structures do not provide an explanation for why the ability to induce the αC-in conformation enables the inhibitor to favor dimeric BRAF^V600E over monomeric BRAF^V600E. Thus, to test the hypothesis regarding the role of h-bond with the Glu501 and to dissect the mechanism of dimer selectivity, we conducted a series of MD simulations of the monomeric and dimeric BRAF^V600E in the absence and presence of two dimer-compatible inhibitors (see below).

Overview of the MD simulations of the monomeric and dimeric BRAF^V600E

The dimer interface of BRAF^V600E contains two histidines, His477 and His510. His510 forms a h-bond with His477 of the opposite protomer, while His477 is also in a potential salt bridge distance from Asp595 of the opposite protomer (Supplemental Figure 2). In a preliminary simulation, where all histidines were set to be neutral and in a tautomer state determined by inspection of the X-ray structure (His477 was set to HID; all others set to the AMBER (Case et al., 2020) default HIE), we found the BRAF^V600E dimer dissociated within a few hundred nanoseconds. To rigorously determine protonation states, we applied the all-atom continuous constant pH molecular dynamics (CpHMD) titration (Harris et al., 2022), which revealed that His477 is most likely in the charged HIP state while His510 is most likely in the neutral HIE state at neutral pH (Supplemental Figure 3,4).

Figure 3.

Figure 4.

Based on the CpHMD determined protonation states, we carried out a series of fixed-charge MD simulations of the monomeric and dimeric BRAF^V600E in the ligand-free state (apo) or in complex with the PHI1 or LY inhibitor in each protomer (holo). To investigate the cooperativity of inhibitor binding, MD simulations were also conducted where only one protomer is complexed with the PHI1 or LY inhibitor (mixed). Each simulation lasted 5 μs and was repeated three times for statistical significance; in total, 135 μs trajectory data was collected (Table 1) and the last 3 μs of each repeat was used for analysis.

Dimerization restrains and shifts αC inward while increasing the flexibility of DFG

In order to understand why an inhibitor prefers binding with a dimer or monomer BRAF^V600E, it is important to understand the difference in the conformation and dynamics between the apo monomeric and dimeric BRAF^V600E. We focus on the αC helix and DFG motif due to their flexibility and importantly specific interactions with the inhibitors (Figure 2). Following KLIFS (Kanev et al., 2020), the αC position is characterized by the distance between Ile582 on β7 (representing a stable reference point) and the center of mass of the Cα atoms of Asn500, Glu501, and Val502 (representing the center of the αC helix); a distance below 19.6 Å defines the αC-in while a distance above defines the αC-out states. We also examined the salt-bridge formation between the αC-Glu501 and catalytic Lys483; a minimum sidechain distance below 4.5 Å is an alternative way to define the αC-in states (Tsai et al., 2019; Sultan et al., 2018). These two definitions are consistent and offers complementary information (see later discussion). The holo PHI1-bound structure (PDB: 6P7G) has both protomers resolved with the αC positions of 19.1 and 19.0 Å, suggesting that the αC helix is in but close to the boundary (19.6 Å) with αC-out according to the KLIFS definition (Kanev et al., 2020). Unlike in the co-crystal structures of dimer-compatible inhibitors, the simulations of the apo monomer and dimer revealed that the αC helix mostly samples the αC-out state. Compared to the apo monomer, the αC position is not only more restrained but also shifted inward by about 1 Å in the apo dimer, as seen from the increase of the peak height and the left-shift of the peak position in the probability distributions, from 23.2 to 22.0 Å (Figure 3a). The flexibility of the αC position in the apo BRAF^V600E is consistent with a previous MD study (Maloney et al., 2021) Enabled by the αC inward movement, the probability of salt-bridge formation between Glu501 and Lys483 is increased by two-fold in the apo dimer (∼25%) relative to the apo monomer (∼12%, Figure 3b). The enhanced but nonetheless unstable Glu501–Lys483 salt bridge indicates that dimerization primes the αC-helix for adopting the αC-in state, e.g., upon interacting with a dimer-compatible inhibitor. In contrast to the αC helix, the motion of the DFG motif is significantly enhanced, as evident from the significant widening of the probability distribution of the DFG pseudo dihedral (Figure 3c), defined by the Cα atoms of Ile592 (DFG-2), Gly593 (DFG-2), Asp594 (DFG-Asp), and Phe595 (DFG-Phe) (Möbitz, 2015). Based on a cutoff of 140°, the DFG pseudo dihedral has been found to discriminate between the DFG-in and DFG-out states of kinases (Möbitz, 2015; Tsai et al., 2019). Accordingly, the distributions indicate that the DFG motif samples the DFG-out state in both apo monomer and dimer, with the DFG pseudo dihedral of ∼210°; however, the dimeric DFG is capable of occasionally sampling the DFG-in state due to the increased flexibility (Figure 3c). While this does suggest dimerization loosens the DFG motif, our simulations do not appropriately model the DFG-out/-in transition as the DFG-in state is only occasionally sampled.

PHI1 and LY binding induces the αC-in state to varying degrees and shifts DFG out

Having understood how dimerization modulates the conformational dynamics of the αC helix and DFG motif, we proceeded to explore conformational changes induced by the dimer-compatible inhibitors PHI1 and LY. Interestingly and as expected, both inhibitors further restrain the motion of the αC helix, with its position sampling a narrower range of 4 Å, as compared to 7 Å in the apo dimer (Figure 3d). Importantly, the αC position is shifted inward by at least 2.7 Å in the holo relative to the apo dimer, and PHI1 induces a larger shift, to 18.3 Å as compared to 19.3 Å in the presence of LY (Figure 3d). The inward shift of the αC helix by the two inhibitors is also reflected in the stabilization of the Glu501–Lys483 salt bridge, which is promoted in the presence of LY (60% vs. 25% in the apo dimer) and is completely locked in the presence of PHI1 (Figure 3e). Although the DFG motif is also significantly restrained through inhibitor binding, the DFG pseudo dihedral in the holo dimer is shifted outward by 80° in complex with either PHI1 or LY (210° in the apo dimer vs. 290° in the holo dimer, Figure 3f).

H-bond formation with Glu501 is critical for dimer selectivity by shifting αC helix inward

The monomer-selective inhibitors do not contact the center of the αC helix and their co-crystal structures only adopt αC-out state (Figure 2b). To test our hypothesis that the h-bond formation with Glu501 is critical for restricting the αC helix to the αC-in states, we examined the density plots of the αC position vs. the distance between the amide nitrogen of PHI1 or LY and the carboxylate of Glu501 in the holo dimer simulations (Figure 3g,h). In the PHI1-bound dimer simulations, the PHI1–Glu501 h-bond is stable with only occasional breakages, as seen from the density maximum centered at the N4–Glu501 distance of 2.9 Å and αC position ∼18 Å (Figure 3g). In the LY-bound dimer simulations, however, the LY–Glu501 h-bond is weaker and less stable than the counterpart of the PHI1-bound dimer, as seen from the local density maximum centered at ∼3.4 and the global maximum near sim4.5 Å (Figure 3g,h). The stronger h-bond between PHI and Glu501 may be attributed to the additional hydrophobic interaction PHI1 forms with Asn500, which is absent for LY (Figure 2b). It is also noteworthy that when the PHI1–Glu501 interaction switches from h-bonding to van der Waals interaction, the αC position is slightly shifted outward to ∼19 Å, which is similar to the position adopted in the LY-bound dimer simulations. This suggests that the stronger h-bond between PHI1 Glu501 may contribute to the inward αC position as compared to the LY-bound dimer.

To further dissect the mechanism of dimer selectivity, we examined the h-bond interaction between PHI1 or LY and Glu501 in inhibited monomer BRAF^V600E simulations. Strikingly, the PHI1– Glu501 interaction can become completely disrupted, with the distance moving beyond 6 Å to as high as 12 Å; correlated with the disruption of the PHI1–Glu501 interaction, the αC position is shifted out to the range of 21 Å–24 Å, similar to that sampled by the apo dimer (Figure 3i). In stark contrast, the LY–Glu501 interaction remains stable as in the holo dimer simulations (Figure 3j). These data are consistent with the previous simulations of the LY- and regorafenib (REG)-bound monomeric and dimeric BRAF^V600E based on different force fields, which showed that the root-mean-square deviation (RMSD) of the dimer-selective REG is increased in the monomer compared to dimer simulations, whereas the RMSD of the equipotent LY remains the same (Adamopoulos et al., 2021).

The correlation between the αC position and the LY–Glu501 interaction confirms our hypothesis that the h-bond interaction between the inhibitor and Glu501 is a key for restraining the αC helix and shifting it to the αC-in states. Since dimerization already restricts the motion of the αC helix and shifts it inward in the apo dimer, inhibitors capable of interacting with Glu501 can bind to the dimer via a conformational selection mechanism in addition to induced fit. On the other hand, conformational selection cannot be exploited for these inhibitors to bind the monomer, as the αC position in the apo monomer is outward. Compared to the equipotent inhibitors, the dimer-selective inhibitors such as PHI1 form much stronger h-bond with Glu501, which shifts the αC helix further inward. The latter may lead to a larger entropic penalty for the monomer binding as compared to the equipotent inhibitors.

PHI1 or LY binding has similar stabilizing effect on the dimer interface of BRAF^V600E

The aforementioned data demonstrates the importance of considering entropic penalty in monomer binding as a contributor to dimer selectivity. To rule out the possibility that the different degree of dimer (de)stabilization may also be a contributing factor for dimer selectivity, we turned to the difference contact network analysis (Yao et al., 2018). In this analysis, the BRAF^V600E dimer was first partitioned into four different communities based on the the residue-residue contacts, which resulted in each community largely corresponding to the N-lobe (blue or grey) and C-lobe (red or orange) of either protomer (Figure 4 left). Then, a community-community difference contact network between the apo and holo dimer simulation sets was calculated and mapped onto a graph, where the vertices represent the communities and blue and red edges represent the increased and decreased contact probabilities due to inhibitor binding (Figure 4 right). Since we are interested in testing the dimer stability in the presence of PHI1 or LY, the interprotomer contact probabilities(between N-lobe:A and N-lobe:B or C-lobe:B; between C-lobe:A and C-lobe:B or N-lobe:A) were calculated and summed up. Interestingly, for both PHI1 and LY, the total interprotomer contact probability is increased (by 0.3) in the holo relative to the apo simulations. This net increase is mainly due to the N-lobe:A to C-lobe:B interactions which compensates for the decrease in the C-lobe:A to C-lobe:B contacts. This analysis demonstrates that both the dimer-selective and equipotent inhibitors have the same slightly stabilizing effect on the BRAF^V600E dimer interface; this rules out the possibility that the dimer selectivity is due to the different degree of dimer stabilization between the dimer-selective and equipotent inhibitors.

Positive cooperativity of PHI1 is due to the allosteric modulation of the αC and DFG conformation in the opposite protomer

As previously mentioned, PHI1 was found to exert a more potent inhibition of the second protomers of the BRAF^V600E dimer whereas LY demonstrated similar potency in the inhibition of the two protomers (Cotto-Rios et al., 2020). To shed light on this cooperativity mechanism, we examined the simulations of the mixed BRAF^V600E dimers in which only one protomer is in complex with PHI1 or LY. We first compared the αC helix position of the apo protomer in the mixed dimers (Figure 5a). Surprisingly, the αC helix of the apo protomer in the PHI1-bound mixed dimer is restrained and shifted inward by 1 Å relative to the apo dimer; in contrast, the position of the corresponding αC helix in the LY-bound mixed dimer remains the same but becomes slightly more flexible (blue and grey, Figure 5a left). Because of the inward shift of the αC helix in the PHI1-bound mixed dimer, the Glu501 of the apo protomer has a 25% higher probability of forming a salt bridge with Lys483 as compared to the apo dimer; in contrast, the salt-bridge probability for the corresponding Glu501 in the LY-bound mixed dimer remains the same as in the apo dimer (blue and grey, Figure 5a middle). These data demonstrate that PHI1 binding in one protomer allosterically modulates the αC helix in the second apo protomer such that it moves inward and becoming more favorable for binding the second PHI1.

Figure 5.

From the distributions of the DFG pseudo dihedral, we can see a slight right shift in the peak position for the DFG in the apo protomer of the PHI1-bound mixed dimer relative to the apo dimer (blue and grey, Figure 5a right). Although the shift is small (the differences between means is approximately one standard deviation, see Supplementary Table 2), it suggests that PHI1 binding in one protomer can allosterically shift the DFG motif outward, making it favorable for binding a second inhibitor. In contrast, the DFG dihedral of the apo protomer in the LY-bound mixed dimer appears to be slightly smaller than the apo dimer with difference between means of approximately one standard deviation (Supplementary Table 2), which is unfavorable for binding the second inhibitor (orange and grey, Figure 5a right). The flexibility of the DFG motif in the apo protomer of the PHI1- or LY-bound mixed dimers is the same as in the apo dimer.

Next, we compared the αC helix position in the holo protomer of the mixed dimers (Figure 5b left). Remarkably, the αC helix in the PHI1-bound protomer of the mixed dimer (blue) is shifted outward by ∼1 Å relative to the holo dimer bound to two PHI1 (green, Figure 5b left), demonstrating that the second PHI1 binding allosterically shifts the αC helix further inward. Further analysis shows that in the holo dimer the αC helix of one protomer is on average 0.5 Å closer compared to the neighboring protomer (Supplementary Table 2). Nonetheless, the Lys483–Glu501 salt bridge is stable in both the mixed and holo dimers; this is because the αC helix predominantly samples the αC-in state in both cases (blue and green in Figure 5b middle). In contrast to PHI1, the αC position in the LY-bound protomer of the mixed dimer is similar to that in the LY-bound holo dimer (orange and red in Figure 5b left), although the probability of the Lys483–Glu501 salt-bridge in the LY-bound protomer in the mixed dimer is slightly lower than in the holo dimer.

Consistent with the effect of the second PHI1 on the αC position of the first PHI1 bound protomer, the second PHI1 allosterically shifts the peak of the DFG probability density further outward, as shown by the 30° larger DFG pseudo dihedral in the holo dimer relative to the mixed dimer (green and blue in Figure 5b right). In contrast, there is no significant difference in the DFG pseudo dihedral between the LY-bound mixed and holo dimers. This data demonstrates that the presence of PHI1 in one protomer modulates the αC and DFG conformation of the apo protomer such that the apo protomer becomes more favorable for binding.

Concluding Discussion

We explored the mechanism of dimer selectivity and cooperativity of BRAF^V600E inhibitors using MD simulations of the dimeric and monomeric BRAF^V600E, in the absence and presence of one or two dimer-selective (PHI1) or equipotent (LY) inhibitor(s). The simulations uncovered the atomic details of the remarkable allostery in BRAF^V600E dimerization and ligand binding (Figure 6), which offer explanation for why some BRAF inhibitors are monomer selective while others are dimer compatible, i.e., selective or equipotent. Specifically, our data showed that dimerization of BRAF^V600E leads to the restriction and an inward shift of the αC helix position relative to the monomer (Figure 6 top panel), which explains why inhibitors that can stabilize the αC-in states are dimer compatible whereas those that cannot are monomer selective. The fact that both dimerization and inhibitor binding induces αC to move inward contributes to the phenomenon of drug-induced RAF dimerization (Hatzivassiliou et al., 2010; Lavoie et al., 2013; Karoulia et al., 2016).

Figure 6.

The co-crystal structure analysis and MD simulations identified a h-bond donor (e.g., an amide linker in the dimer-selective PHI1 or the equipotent LY) as a key for dimer compatibility; the h-bond with the carboxylate of Glu501 stabilizes the αC helix in the αC-in states. Two factors make Glu501 a special and critical anchoring point for inducing the αC-in states. First, it is located at the center of the αC helix, which makes it easier (as opposed to the end of the helix) to induce a helix movement. Second, the h-bonding between the inhibitor and Glu501 is synergistic with the Lys483–Glu501 salt-bridge formation. In contrary, the lack of a h-bond with Glu501, e.g., in VEM, Debrafenib, or PLX7904, results in the monomer selectivity. Note, the DFG-in inhibitors can also donate a h-bond to αC-Glu501 and bind the BRAF^V600E dimer. An example is SB5909885, which donates a h-bond from the oxime group to αC-Glu501 and also forms a salt bridge with Lys483 (PDB ID: 2FB8) (King et al., 2006).

The difference between the dimer-selective and equipotent inhibitors is more subtle. The MD simulations revealed that PHI1 forms a more stable h-bond with Glu501 in the BRAF^V600E dimer as compared to LY, which is consistent with the ∼1 Å inward shift of the αC helix and more stable Lys483–Glu501 salt bridge. The latter differences are much smaller in the co-crystal structures; the αC positions and Lys483–Glu501 distances are only respectively 0.1 and 0.2 Å smaller in the PHI1-vs. LY-bound co-crystal structure. Since the monomeric BRAF^V600E has a flexible αC helix that predominantly samples the αC-out states, forming a tighter h-bond would incur a higher entropic penalty for monomer binding. This may explain why the PHI1–Glu501 interaction as well as the αC position are unstable in the monomer simulations but stable in the dimer simulations, in contrast to the LY-bound simulations. Therefore, the stability of the h-bonding with Glu501 may be a key for dimer selectivity.

Without the MD simulations, how would one determine if the h-bond between the inhibitor and Glu501 is stable? We found that the deviation between the αC position and/or K–E distance of the two protomers in the co-crystal structure offers some indication (Supplemental Table 1). With the exception of LY and Ponatinib, the αC position and/or K–E distance between the two protomers in the co-crystal structures of AZ628, TAK632, BGB283, SB5909885 deviate by 0.3 Å or higher (Supplemental Table 1). In contrast, the αC position and the K–E distance are (nearly) identical between the two protomers in the co-crystal structures of the dimer-selective inhibitors LXH254, RAF709, Sorafenib, and Belvarafenib (Supplemental Table 1). The identical αC position and K-E distance in the two protomers suggest that the αC helix is restrained by the inhibitor, i.e. it forms a stable h-bond with Glu501.

To additionally test this crystal structure-based hypothesis, we examined the co-crystal structures of GDC0879 and Tovorafenib, which were not analyzed in Ref (Adamopoulos et al., 2021). In the co-crystal structure of GDC0879 (PDB ID: 4MNF), the αC position deviates by 0.3 Å and the K–E distance deviates by 0.1 Å between the two protomers. In the co-crystal structure of Tovo-rafenib (PDB ID: 6V34), the the αC position deviates by 0.2 Å and the K–E distance deviates by 0.4 Å between the two protomers. These deviations suggest that the αC is not adequately restrained by the inhibitors and therefore we predicted GDC0879 and Tovorafenib to be equipotent. Note, GDC0879 is a DFG-in inhibitor, which is an additional indication for a equipotent inhibitor. Indeed, both GDC0879 and Tovorafenib were found as equipotent in experimental studies (Karoulia et al., 2016; Tkacik et al., 2023). These analyses led us to propose the following empirical assessment of a RAF inhibitor based on its co-crystal structure with BRAF^V600E: 1) lack of a h-bond with Glu501 indicates monomer selectivity; 2) presence of a h-bond with Glu501 but inconsistent αC position and/or K–E distance between the two protomers indicates equipotency; 3) presence of a h-bond with Glu501 and identical αC position and K–E distance between the two protomers indicates that the inhibitor is likely (but not necessarily) dimer selective. Given that the resolution of a resolved structure is often ∼2-3 Å, this proposed assessment is not intended to replace more rigorous tests, i.e. utilizing MD simulations.

Finally, the MD analysis uncovered a mechanism for positive cooperativity. Our findings are summarized in Figure 6; upon dimerization (top row) the αC-helix goes from αC-out and highly flexible to slightly restrained and inward shifted. The mixed simulations demonstrated that the first PHI1 binding in the BRAF^V600E dimer primes the second apo protomer by making the αC and DFG conformation more favorable for binding, i.e., shifting the αC inward and the DFG outward (Figure 6, bottom right panel). Importantly, without a second PHI1, the αC and DFG conformation in the first protomer is not fully shifted in or out, respectively, as compared to the two-inhibitor bound dimer (Figure 6, bottom left panel). These data suggest that the positive cooperativity of PHI1 is due to its ability to allosterically modulate the αC and DFG conformation in the second protomer. Taken together, our findings provide a mechanistic understanding for the remarkable allostery and conformational interplay between kinase dimerization and inhibitor binding. As we prepare the manuscript for submission, a biophysical experiment was published, which suggested that the first inhibitor binding dominates the allosteric coupling between type II inhibitor binding and BRAF dimerization (Rasmussen et al., 2023), consistent with our simulation data. The work presented here has implications for understanding the molecular mechanism of kinase signaling and contributes to the rational design of protomer-selective inhibitors.

Methods and Protocols

Intracellular homogeneous TR-FRET assay

SKMEL239 and SKMEL239-C4 cells were plated at 50000cells/well in white TC-treated 96-well plates in 100ul complete growth media (DMEM). Cells were incubated with the various RAF inhibitors for 1 hour at 37°C, 5% CO₂. ERK phosphorylation was measured using the THUNDER™ Extreme Phospho-ERK1/2 (T202/Y204) TR-FRET Cell Signaling Assay Kit (Bioauxillium) according to directions for the Standard 2-Plate Assay Protocol for Adherent Cells. Cells were lysed for 30 minutes at RT under shaking. Lysates were transferred to a white 384-well plate, sealed and incubated with the detection mix antibody at RT for 4 hours. TR-FRET signal was measured at 615 nm and 665 nm excitation using a TECAN SPARK plate reader.

System preparation for simulations

Simulations were prepared using a crystal structure of BRAF^V600E in complex with either PHI1 (PDB ID: 6P7G) (Cotto-Rios et al., 2020) or LY (PDB ID: 5C9C) (Peng et al., 2015). The initial structure of the apo simulations was taken from 6P7G. The a-loop is not resolved in either protomer in 5C9C, but is resolved for protomer B in 6P7G. Thus, the missing a-loop in the protomer B of 6P7G and in both protomers in 5C9C were built by rotating and translating the resolved a-loop from the first protomer using the alignment tool in PyMOL (Schrödinger, LLC, 2015). Chain B (which has the resolved a-loop) was first aligned to chain A using all residues except for the a-loop and the two end residues that connect it to the rest of the protein. Following the alignment, chain B except for the a-loop and its two end residues were deleted. The N-terminus was acetylated and the C-terminus was amidated. Hydrogen atoms were added using the HBUILD facility in the CHARMM package (version c37a2) (Brooks et al., 2009).

All-atom continuous constant pH molecular dynamics (CpHMD) simulations

The recently developed all-atom particle mesh Ewald CpHMD (PME-CpHMD) (Harris et al., 2022) with the asynchronous pH replica exchange sampling protocol (Wallace and Shen, 2011; Henderson et al., 2020) was used to determine the protonation and tautomer states of histidines. To prepare for the CpHMD simulations, the histidine residues were first set to HIP with the dummy hydrogens on the Nδ and Nϵ atoms. The system was solvated in a rectangular water box with at least 10 Å distance between the protein and the boundary (∼23,000 water molecules). The protein was represented by the AMBER ff14SB force field (Maier et al., 2015) and water by the TIP3P model (Jorgensen et al., 1983). The dimer structure was briefly minimized for 500 steps (first 200 were using steepest decent, following 300 used conjugate gradient) with a harmonic force constant of 100 kcal/mol/Ų applied on all heavy atoms of the protein. This was followed by 100 ps of heating to 300 K using the PME-CpHMD simulations at pH 7.0 with the restraints still applied. Once heated the restraints were gradually removed in six stages: in the first two stages the protein heavy atoms were restrained with a force constant of 100 and 10 kcal/mol/Ų; in the next four stages only the backbone heavy atoms were restrained with a force constant of 10, 1.0, 0.1, and 0.0 kcal/mol/Å. Each stage was simulated for 250 ps, for a total of 1.5 ns. A cutoff of 12 Å was used for the non-bonded interactions.

The equilibrated structure was then used to initiate the pH replica exchange PME-CpHMD simulations. The asynchronous pH replica exchange sampling protocol (Wallace and Shen, 2011; Henderson et al., 2020) was used to accelerate convergence of the coupled protonation and conformational states (Wallace and Shen, 2011). Five replicas were created at different pH conditions, from pH 6.5 to 8.5. Each replica was first equilibrated to its pH by repeating the final four stages of equilibration mentioned above. The pH replica exchange CpHMD was then conducted for 10 ns with attempted swaps of neighboring pH conditions occurring every 2 ps. All other settings are identical to Ref. Harris et al. (2022). For the calculation of protonation and tautomer state probabilities, the λ and x values above 0.8 or below 0.2 were used (default setting in the CpHMD analysis package (Henderson et al., 2022)). At pH 7.5 His477 was protonated at both Nϵ and Nδ while His510 was protonated at Nϵ only. These protonation/tautomeric states were used for all convention (fixed-protonation-state) simulations below.

Conventional fixed-protonation-state MD simulations

Eight BRAF^V600Esystems were simulated, consisting of monomeric and dimeric BRAF^V600E either in the presence or absence of PHI1 or LY (see Table 1). Monomer systems were prepared by eliminating one protomer from the prepared dimer structure. In the apo monomer and dimer systems, ligand(s) was removed. In the mixed or holo systems, one or both inhibitors from the co-crystal structure was kept. The protein was then placed in a rectangular water box with a minimum distance of 10 Å between the protein and edges of the water box using the LEaP program (Case et al., 2020). Based on the protonation states determined using CpHMD, sodium and chloride ions were added to neutralize the system and reach a physiological ionic strength of 0.15 M.

The conventional (fixed-protonation-state) MD simulations were carried out using the AMBER20 MD package (Case et al., 2020). The proteins was represented by the ff14SB force field (Maier et al., 2015) while inhibitors were parameterized by the general AMBER force field (GAFF) method (Wang et al., 2004). The TIP3P model (Jorgensen et al., 1983) was used to represent water. The Leapfrog integrator was used to propagate the coordinates. The SHAKE algorithm was applied to bonds involving hydrogen to allow for a 2-fs time step. Additionally, the hydrogen mass re-partitioning (Hopkins et al., 2015) was used to redistribute the mass between hydrogens and their bonded heavy atoms to allow for a 4-fs time step. A nonbonded cutoff of 8 Å was used as in the ff14SB validation study (Maier et al., 2015) while the electrostatic potentials were computed using the particle-mesh Ewald method (Darden et al., 1993) with a real-space cut-off of 12 Å and a sixth-order interpolation with approximately 1 Å grid spacing. Each system underwent minimization using 1000 steps of steepest descent followed by 19000 steps of conjugate gradient while the heavy atoms were harmonically restrained using a force constant of 100 kcal/mol/Ų. Following minimization, the system was heated to 300 K over 1 ns under an NVT ensemble using a Langevin thermostat (Feller et al., 1995) with collision frequency of 1 ps^-1 for temperature control. The systems then underwent a 6-stage equilibration in which the backbone restraints were gradually reduced to 10, 5, 2, 1, 0.1 and 0 kcal/mol/Ų over the course of 100 ns under a NPT ensemble. A Monte-Carlo barostat (Case et al., 2020) was used to control pressure at 1 bar using a relaxation time of 1.0 ps. Each system were run in three replicates each starting from different random velocity seeds and each run lasted 5 μs.

Simulation data analysis

CPPTraj (Roe and Cheatham, 2013) was used to analyze the protomer conformation (αC-helix position, DFG pseudo dihedral, etc.) and visualizations were produced using PyMOL (Schrödinger, LLC, 2015). The contact network analysis was conducted using the open source code developed by Yao and Hamelberg (https://github.com/The-Hamelberg-Group/dcna) (Yao et al., 2018). Unless otherwise noted, the last three μs trajectory frames were used for analysis. All probability distributions were created by combining the last three μs of each replica for each system, with each distribution consisting of 50 bins. Unless specified, distributions contain quantities from both protomers in dimeric simulations.

Data Availability

The MD simulation input files and analysis scripts are freely downloadable from https://github.com/JanaShenLab/RAF/. The raw MD trajectories are available upon request.

Acknowledgements

Funding support by the National Cancer Institute to J.S. (R01CA256557) and E.G. (R01CA238229) is acknowledged.