Exploring Molecular Mechanisms of Paradoxical Activation in the BRAF Kinase Dimers: Atomistic Simulations of Conformational Dynamics and Modeling of Allosteric Communication Networks and Signaling Pathways.

The recent studies have revealed that most BRAF inhibitors can paradoxically induce kinase activation by promoting dimerization and enzyme transactivation. Despite rapidly growing number of structural and functional studies about the BRAF dimer complexes, the molecular basis of paradoxical activation phenomenon is poorly understood and remains largely hypothetical. In this work, we have explored the relationships between inhibitor binding, protein dynamics and allosteric signaling in the BRAF dimers using a network-centric approach. Using this theoretical framework, we have combined molecular dynamics simulations with coevolutionary analysis and modeling of the residue interaction networks to determine molecular determinants of paradoxical activation. We have investigated functional effects produced by paradox inducer inhibitors PLX4720, Dabrafenib, Vemurafenib and a paradox breaker inhibitor PLX7904. Functional dynamics and binding free energy analyses of the BRAF dimer complexes have suggested that negative cooperativity effect and dimer-promoting potential of the inhibitors could be important drivers of paradoxical activation. We have introduced a protein structure network model in which coevolutionary residue dependencies and dynamic maps of residue correlations are integrated in the construction and analysis of the residue interaction networks. The results have shown that coevolutionary residues in the BRAF structures could assemble into independent structural modules and form a global interaction network that may promote dimerization. We have also found that BRAF inhibitors could modulate centrality and communication propensities of global mediating centers in the residue interaction networks. By simulating allosteric communication pathways in the BRAF structures, we have determined that paradox inducer and breaker inhibitors may activate specific signaling routes that correlate with the extent of paradoxical activation. While paradox inducer inhibitors may facilitate a rapid and efficient communication via an optimal single pathway, the paradox breaker may induce a broader ensemble of suboptimal and less efficient communication routes. The central finding of our study is that paradox breaker PLX7904 could mimic structural, dynamic and network features of the inactive BRAF-WT monomer that may be required for evading paradoxical activation. The results of this study rationalize the existing structure-functional experiments by offering a network-centric rationale of the paradoxical activation phenomenon. We argue that BRAF inhibitors that amplify dynamic features of the inactive BRAF-WT monomer and intervene with the allosteric interaction networks may serve as effective paradox breakers in cellular environment.

Introduction

The human protein kinases are involved in regulation of many functional processes in signal transduction networks and represent one of the largest classes of clinically important therapeutic targets []. Protein kinases act as versatile activators and dynamic regulatory switches that are essential for regulation of cell cycle and organism development. A staggering amount of structural, genetic, and biochemical data on protein kinase genes has been accumulated in recent years, revealing a large variety of regulatory mechanisms, ranging from phosphorylation of kinase activation loops and autoinhibition to allosteric activation induced by dimerization or protein binding [. The steadily growing structural knowledge about conformational states of the kinase catalytic domain, regulatory assemblies, and kinase complexes with small molecule inhibitors has provided compelling evidence that conformational transformations between the inactive and active kinase states are central to the enzyme regulation and function []. Functional conformational changes in protein kinases are operated by several regulatory regions of the catalytic domain: the conserved catalytic triad His-Arg-Asp (HRD), the DFG-Asp motif, the regulatory αC-helix, and the activation loop (A-loop). The inactive kinase states are often characterized by the DFG-out and closed A-loop conformations, while the active kinase forms feature the DFG-in and open A-loop conformations []. These regions are also involved in the formation of the regulatory spine (R-spine) and catalytic spine (C-spine) networks that are assembled and stabilized during conformational transformations to the active kinase states []. Despite diversity of regulatory mechanisms, modulation of kinase activity through dimerization and conformational repositioning of the αC-helix emerged as a common mechanism shared by several important protein kinase families, including ErbB kinases [] and BRAF kinases []. Structural determinants of dimerization-induced regulation in the ErbB and BRAF kinases are rather similar, as the “off-state” of both enzymes is defined by a nonproductive αC-helix-out conformation supported by the presence of a short helical element in their A-loops that locks the enzyme in the inactive dormant form. Dimerization-induced allosteric regulation involves coordinated transitions of the kinase domain from the inactive monomer structure to a dimer configuration in which the αC-helix moves to an active ‘in’ conformation that ensures a productive alignment of the hydrophobic spines and catalytic residues required for activation. While a head-to-tail dimer assembly of the catalytic domains is characteristic of the ErbB kinases [], a symmetric side-to-side dimer arrangement represents structural ‘modus operandi’ of the BRAF kinase activation []. The crystal structure of the inactive BRAF kinase has revealed a non-productive monomeric state of the enzyme, in which the αC-helix-out conformation can disrupt structural environment of the catalytic and regulatory residues near the ATP-binding site that is required for activation [. Dimer-inducing BRAF inhibitors irrespective of their binding modes may restrict the inter-lobe dynamics of the catalytic domains and promote stabilization of the active kinase conformations that facilitate the effective side-to-side dimerization []. Curbing the initial enthusiasm of the BRAF drug discovery efforts, the recent breakthrough studies have revealed that most of the existing BRAF inhibitors can paradoxically activate the wild-type (WT) BRAF kinase by inducing conformational changes and promoting dimerization that triggers RAS-independent transactivation and leads to MEK/ERK phosphorylation [ The ability of dimer-inducing BRAF modulators to effectively inhibit the first monomer while conferring an active form of the second monomer may inadvertently stimulate kinase activation and promote ERK signaling []. As a result, therapeutic and clinical utility of BRAF inhibitors that exert paradoxical activation can be diminished.

Based on the crystallographic binding mode, kinase inhibitors have been classified in several distinct classes []. Type I kinase inhibitors target the catalytically competent, active (DFG-in/ αC-helix-in) conformation of the kinase domain. Type II kinase inhibitors recognize the inactive DFG-out/αC-helix-in kinase state, while type I½ inhibitors often bind to another inactive DFG-in/αC-helix-out conformation that precludes the formation of the active enzyme form. Type II and type I½ kinase inhibitors often exhibit a higher degree of selectivity as these molecules bind to specific inactive conformations that may be structurally unique for a given kinase target. Many type I BRAF inhibitors have been developed and the crystal structures of their complexes with the BRAF dimers have been solved (), including SB-590885 (pdb id 2FB8) [] and GDC-0879 (pdb id 4MNF) []. Although type I BRAF inhibitors are effective against the oncogenic BRAF-V600E mutant, these molecules are usually strong dimer-promoting agents causing paradoxical activation in the BRAF-WT cells. Type II BRAF inhibitors [] that stabilize the inactive DFG-out/αC-helix-in conformation can also stimulate enzyme transactivation by favoring an αC-helix-in position that promotes side-to-side dimerization (). The crystal structures of the BRAF-WT and BRAF-V600E complexes with the FDA-approved Sorafenib (pdb id 1UWH, 1UWJ) have demonstrated that type II inhibitors that strongly bind to both monomers and exhibit long residence time may be relatively weak dimer inducers []. Another dimer-inducing type II inhibitor LY3009120 (pdb id 5C9C) is a pan-RAF inhibitor of all three RAF isoforms that causes only minimal paradoxical activation due to structural occupation and effective inhibition of both monomers []. Other notable type II BRAF inhibitors include TAK-632 [], RAF-625 which is now advanced into clinical trials [], and AZ-628 []. Structural studies of the type II BRAF agents have demonstrated that their inhibitory effects may not directly correlate with a dimer-promoting potential and induction of paradoxical activation.

Type I½ BRAF inhibitors target specific DFG-in/αC-helix-out kinase conformations and have exploited structural features of the inactive BRAF-WT monomer with a nonproductive αC-helix-out position []. Several BRAF structures complexed with type I½ inhibitors feature a small helical segment in the A-loop, which is present in the crystal structure of the inactive BRAF-WT monomer []. This class of BRAF inhibitors has offered unique opportunities for targeted modulation of the dimerization mechanism by sequestering specific kinase conformations with αC-helix-out arrangements that compromise the fidelity of the dimerization interface and may therefore preclude effective signaling. The crystal structure of the BRAF complex with type I½ inhibitor PLX4720 (pdb id 3C4C) has revealed an asymmetric dimer, in which one of the two binding sites was preferentially occupied by the inhibitor targeting the DFG-in/αC-helix-out kinase conformation, while the inhibitor adopts a different binding mode and binds to a DFG-out/αC-helix-out conformation in the second binding site []. The crystal structure of a related sulfonamide-containing compound PLX3203 bound to BRAF-V600E mutant (pdb id 4FK3) [] has revealed a structural arrangement in which the inhibitor binds to the DFG-in/αC-helix-out conformation of the first monomer, while the second monomer is ligand-free and assumes an active DFG-in/αC-helix-in conformation (Figs ). Vemurafenib is a type I½ BRAF inhibitor that is FDA-approved for the treatment of highly prevalent BRAF-V600E dependent metastatic melanomas []. However, Vemurafenib can promote dimerization and transactivate drug-free protomers, leading to paradoxical pathway activation and enhanced signaling in the BRAF-WT cells. The crystal structure of BRAF-V600E oncogenic mutant in complex with Vemurafenib (PLX4032) showed that the inhibitor binding to a DFG-in/αC-helix-out kinase conformation of the first monomer induced a drug-free DFG-in/αC-in active conformation of the second monomer (Figs ) []. Dabrafenib is another type I½ inhibitor that binds to both monomers with similar DFG-in/αC-helix-out conformations (pdb id 4XV2, 5CSW) (Figs ) []. Dimer-promoting potential of these BRAF inhibitors may be associated with the extent of paradoxical activation. Indeed, while Dabrafenib is a relatively strong inducer of dimerization and paradoxical activation, Vemurafenib may cause a weaker dimer formation and a reduced paradoxical activation in biological assays [].

Chemical Structures and Structural Binding Modes of the BRAF Inhibitors.

Chemical structures of studied BRAF inhibitors: (A) PLX4720, (B) Vemurafenib (B), PLX7904 (C) and Dabrafenib (D). (E) The binding modes of the BRAF inhibitors in the first binding site (monomer A). PLX4720 (in green), Vemurafenib (in red), PLX7904 (in blue), and Dabrafenib (in cyan). (F) The binding modes of the BRAF inhibitors in the second binding sites (monomer B). Among studied inhibitors, only PLX4720 (in green) and Dabrafenib (in cyan) bind to the second monomer. Note structural similarity of the inhibitor binding mode and DFG-in/αC-out kinase conformation in the first monomer, while alternative inhibitor binding modes and conformational variability in the second monomer.

The recent discovery of BRAF inhibitors that can evade paradoxical activation has shown that small changes in the type I½ inhibitors may produce markedly different functional effects by disrupting kinase activation and inhibiting signaling in the BRAF-ERK-MAPK pathways []. The prime examples of these inhibitors are PLX7904 and its optimized analogue PLX8394 that are very similar to Vemurafenib in their chemical structure. These inhibitors selectively bind and inhibit the activity of both BRAF-WT and mutated BRAF forms, showing efficacy against tumors that are resistant to other BRAF agents. The crystal structures of the BRAF-V600E complex with PLX7904 (pdb id 4XV1), the BRAF-V600E mutant bound with PLX7922 (pdb id 4XV3), and the BRAF-PLX5568 complex (pdb id 4XV9) have revealed a binding mode and interaction profiles that are highly reminiscent of the BRAF-Vemurafenib complex. A weak paradox inducer Vemurafenib and a paradox breaker PLX7904 bind to a similar DFG-in/αC-helix-out kinase conformation in the first monomer, while inducing a drug-free DFG-in/αC-in active conformation of the second monomer (Figs ). By combining structure-functional analysis and computational modeling, a recent study has argued that small chemical and structural differences between paradox inducers and paradox breakers may produce appreciable changes in the collective dynamics of the BRAF dimer complexes []. According to this study, paradox inducers may enhance concerted movements of the αC-helix and A-loop regions, while paradox breakers could restrict the αC-helix dynamics and promote more flexible inter-lobe motions. Collectively, structural and functional studies have demonstrated that different classes of BRAF inhibitors can uniquely exploit the mechanism of BRAF activation and act as allosteric modulators of dimerization-induced regulation.

While many computational investigations have been conducted to elucidate molecular mechanisms of the ErbB kinases [], molecular modeling and simulation studies of BRAF kinase dynamics, binding and regulation have been fairly scarce []. Most of these studies have analyzed conformational dynamics of the BRAF kinase catalytic domain and have not thoroughly addressed the molecular basis of ligand-induced BRAF dimerization and paradoxical activation. Molecular mechanisms underlying the effects of inactivating BRAF mutations were explored in molecular dynamics (MD) simulations, showing that that the loss of BRAF function upon G596R alteration may depend on cooperative structural changes near the catalytic residue D594 of the DFG motif []. A combination of docking, MD simulations and QM/MM calculations was used to quantify binding and activity of type I BRAF inhibitors []. MD simulations and binding free energy calculations have been employed in computational studies of type II DFG-out inhibitors in complexes with the BRAF-V600E mutant []. Enhanced MD simulations have suggested a mechanism in which oncogenic BRAF mutations could kinetically trap the kinase in the active form by increasing the barriers of conformational transitions between the inactive and active kinase states []. By combining MD simulations and protein structure network analysis, we have recently examined binding of BRAF dimers with several type I and type II inhibitors []. Despite the growing effort to understand mechanisms of BRAF regulation, there is a significant gap between detailed structural information about the BRAF dimer complexes and limited biophysical and dynamic characterization of the dimer-promoting potential of BRAF inhibitors and its linkage to paradoxical activation.

In the current study, we expanded our previous work [] by systematically exploring molecular mechanisms of ligand-induced BRAF dimerization and paradoxical activation by a representative panel of type I½ inhibitors. The primary focus of this work was to identify molecular drivers that could differentiate paradox inducers PLX4720, Vemurafenib and Dabrafenib from a new generation of paradox breaker inhibitors exemplified by PLX7904. We characterized differences in the conformational dynamics and collective motions of the BRAF dimer complexes to show that paradox breaker may enforce structural and dynamic preferences of the inactive BRAF-WT monomer, which may be required to escape paradoxical activation. Allosteric regulation of multimeric protein complexes is often associated with ligand-induced cooperative effects. While positive cooperativity is a common mechanism for increasing the binding potential [], negative cooperativity occurs when ligand binding to the first binding site decreases the binding affinity in a second binding site, which is often accompanied by conformational entropy changes []. We reported the binding free energy analysis of the BRAF dimer complexes revealing that negative cooperativity and dimer-promoting potential of the type I½ inhibitors may be important molecular determinants of paradoxical activation. The relationships between inhibitor binding, protein dynamics and allosteric signaling in the BRAF dimers were also explored using a network-centric approach. A graph-based representation of protein structures was used as a framework for structural analysis of the residue interaction networks [], identification of functional residues [ and modeling of allosteric communications []. To characterize functional sites of allosteric regulation, we also explored coevolutionary relationships between BRAF residues. According to statistical coupling analysis (SCA) [] and mutual information (MI) model [] functional residues in protein systems are often coevolutionary connected. Networks of residues with high mutual information can characterize structural proximity of functionally important sites and this evolutionary signature can distinguish conserved functional residues []. Moreover, coevolving residues are often located close to each other in the protein structure [] and may form independent structural modules (or protein sectors) with distinct biochemical functions []. Functional conformational changes could be often enabled through networks of evolutionarily coupled residues [].

In the current work, we introduce a generalized network model of protein structure and dynamics in which coevolutionary residue dependencies and dynamic maps of residue correlations are integrated in the construction of the residue interaction networks. Different metrics that characterize topology of these networks were compared to understand mechanisms by which BRAF inhibitors selectively modulate the distribution of global mediating centers and allosteric signaling. We show that dimer-promoting inhibitor potential may be coupled with the induction of paradoxical activation through dynamic changes in the residue interaction networks controlled by the αC-helix and dimer interface sites. The results of this study suggest that molecular determinants underlying unique effects of paradox inducers and breakers may be associated with the differences in robustness and efficiency of allosteric interaction networks in the BRAF complexes. By simulating ensembles of allosteric communication pathways in the BRAF dimer complexes, we establish that paradox inducer and breaker inhibitors may activate specific signaling routes that correlate with the extent of paradoxical activation.

Results and Discussion

Atomistic and Coarse-Grained Simulations of the BRAF Dimers Reveal Ligand-Induced Modulation of Collective Dynamics

All-atom MD simulations of the BRAF dimers in complexes with a panel of type I½ inhibitors were carried out to determine whether conformational dynamics of the BRAF structures would differentiate between paradox inducer and breaker inhibitors. To characterize the inhibitor-induced changes, we compared conformational dynamics profiles of the BRAF dimer complexes with respect to the inactive BRAF-WT monomer (). We particularly examined conformational variations of the regulatory motifs that are responsible for large conformational transformations in the kinase structures: αC-helix (residues 491–505), αC-β3-loop (residues 485–490), αC-β4-loop (residues 506–516), and the A-loop (residues 593–616). Although the inactive kinase states are generally flexible and could exhibit dynamic inter-lobe movements coupled with positional changes of the αC-helix, the monomeric structure of the inactive BRAF-WT is stabilized in its dormant state through additional hydrophobic interactions between the αC-helix-out and a short helical segment of the A-loop []. We found that binding of paradox inducers may result in the increased mobility of the kinase domain in the BRAF dimer complexes (). Moreover, conformational variations in the second monomer were typically greater in the BRAF dimer structures. In particular, thermal fluctuations of the αC-helix and the adjacent αC-β3-loop were markedly larger in the BRAF dimer structures as compared to the inactive BRAF-WT monomer. As a result of these dynamic changes, paradox inducers may promote αC-helix 'in–out' oscillations in the second monomer, leading to the increased population of the closed states with a properly aligned dimer interface. By elevating flexibility of the kinase domain, type I½ paradox inducers may release the constraints of the inactive monomeric form and promote dimerization, which could offset their inhibitory effects. Instructively, binding of paradox breaker PLX7904 seemed to incur dynamic changes that enforce structural preferences of the inactive BRAF-WT monomer in the drug-bound monomer, while increasing conformational mobility of the drug-free monomer (). In particular, partly suppressed fluctuations of the αC-β3-loop and the αC-helix in the drug-bound conformation were comparable to the ones seen in the inactive BRAF-WT monomer. Accordingly, PLX7904 binding may restrain flexibility of the αC-β3-loop and encourage the αC-helix region to maintain its nonproductive ‘out’ position, therefore closely mimicking structural environment of the inactive BRAF-WT monomer. At the same time, the concomitantly increased mobility of the drug-free monomer could compromise cooperative interactions between the two monomers and weaken the dimerization interface. According to the recent studies, modifications in the αC-β3-loop that reduce its length could modulate conformational preferences of the αC-helix causing resistance to type I½ inhibitors []. Our findings corroborate with these experiments, suggesting that the dimer-promoting potential of paradox inducers may be linked with cooperative structural shifts of the αC-helix.

Conformational Dynamics Profiles of the BRAF Structures.

Conformational dynamics profiles of the BRAF dimer complexes with PLX4720 (pdb id 3C4C) (A), Dabrafenib (pdb id 4XV2) (B), Vemurafenib (pdb id 3OG7) (C), and PLX7904 (pdb id 4XV1) (D). Conformational mobility profiles of the BRAF dimer complexes are plotted with the reference to the dynamics profile of the BRAF-WT monomeric structure (pdb id 4WO5). The computed B-factors are annotated and colored as follows: the BRAF-WT monomer (in green), the first monomer of the BRAF dimer (in red) and the second monomer of the BRAF dimer (in blue). The positions of the dimer interface residues (D449, W450, K475, W476, H477, D479, R506, R509, H510, F516, Q562, D565, and Y566) are shown on the dynamics profiles filled maroon circles.

We also analyzed how inhibitor binding may affect conformational dynamics of the key dimer interface residues (D449, W450, K475, W476, H477, D479, R506, R509, H510, F516, Q562, D565, and Y566). The dimer interface is formed by several spatially separated interaction sections of the catalytic core. In our simulations, these regions showed moderate thermal fluctuations consistent with their role in modulating stability of the dimer interface. Of particular interest was the dynamics of the central interface section that includes residues D449, W450, R506, R509, and F516. The regulatory residues from the αC-β4 region (R509 and F516) are stable in both monomer and dimers structures and correspond to the hinge sites that coordinate positional shifts of the αC-helix (). However, the interface residues from the N-terminal tail region (D449, W450) could enjoy a considerable flexibility in the unliganded monomer but become rigidified upon dimer formation. Interestingly, some residues from the peripheral interface section (H577, R558, R562, D565, and Y566) may be also stable in the monomeric form and structurally preorganized for productive inter-monomer contacts in the dimers. Hence, spatially separated interface sections may be predisposed for dimerization, suggesting that allosteric interactions may be transmitted via of diverse communication routes. Based on these findings, we hypothesize that paradox inducer and breakers could selectively modulate stability of the dimerization interface regions and activate specific signaling pathways. To summarize, despite a similar pattern of the conformational fluctuation profiles in the BRAF dimer complexes, we detected specific dynamic signatures that could differentiate paradox inducer and paradox breaker inhibitors. The central finding of this analysis was that paradox breaker PLX7904 in its drug-bound conformation could strictly impose structural and dynamic features of the inactive BRAF-WT monomer that may be required for evading paradoxical activation.

Atomistic simulations suggested that the type I½ inhibitors may induce changes in the collective dynamics of the BRAF dimer structures that may be associated with their dimer-promoting abilities and paradoxical activation. To characterize collectivity of protein motions in the BRAF dimers, we complemented all-atom MD simulations with elastic network modeling (ENM) that can describe large cooperative transformations of protein regions in the space of slow low-frequency modes. According to this model, salient properties of protein dynamics can be uniquely defined by structural architecture of the protein fold and ligand-induced conformational changes may be embedded in the native fluctuations of the unbound proteins []. The essential mobility profiles of the BRAF structures were obtained using the Gaussian network model (GNM) []. In this model, the distribution peaks point to the protein regions involved in large structural movements, while local minima correspond to the functional residues that remain immobilized and serve as global hinge centers that coordinate collective dynamics. We first analyzed collective dynamics of the BRAF-WT monomeric state in the space of the slowest low-frequency modes (). The first principal mode of the BRAF-WT kinase domain corresponded to coordinated rotations of the N-terminal lobe and C-terminal lobe in opposite directions. In this mode, the distance fluctuations between residues within each lobe are small and residue movements are mostly positively correlated, while the inter-residue distance deviations between lobes are large as motions of these residues are negatively correlated (). In the second slow mode, the kinase lobes undergo coordinated opening and closing movements with respect to each other, which may be important for ATP and ligand binding. In the third principal mode, the kinase lobes undergo a shear motion sliding along their interface section which is accompanied by the forward move of one lobe and a concomitant inward displacement of the other lobe. We argue that the conservation of major slow modes in the kinase domain may be required for induction of functional motions during diverse regulatory processes. These principal motions of the catalytic domain are conserved across the protein kinase family as was affirmed in the NMR experiments [] and reproduced in computational studies [].

Functional motions along the slowest modes showed interesting differences between the BRAF dimer structures (). The principal movements of the first and second monomers were similar in the BRAF dimers complexed with paradox inducers PLX4720 and Dabrafenib. However, in the BRAF complexes with Vemurafenib and PLX7904, collective movements of the second monomer were reduced near the αC-β3-loop (residues 485–490) and secondary interface section (residues H510, K475, W476, H577, R558, R562, D565, and Y566) (). Global motions around this hinge site may compromise effective dimerization and reduce the extent of paradoxical activation. By examining collectivity of the slow modes in the BRAF dimer structures we found that up to five low-frequency modes may be required to adequately describe all global movements (). Subsequent slow modes have a markedly reduced level of collectivity and could only affect local dynamics changes (). The degree to which these slow GNM modes capture global structural changes was estimated by the cumulative overlap with the principal component modes derived from conformational variability in the BRAF dimer crystal structures []. The slowest five modes exhibited more than 90% overlap with the principal components of structural changes observed in the BRAF dimer complexes with different classes of inhibitors ( Based on the observed convergence of slow modes, we analyzed the essential mobility profiles from the cumulative contribution of the first five principal modes and computed the mean square displacements of protein residues driven by these slow modes. We observed that dimer formation may introduce an additional cooperative mode with an extended global hinge that runs across the entire inter-monomer interface (Figs ). This slow mode forms a dominant direction of collective dynamics in the BRAF dimers, where the rigid interfacial sites become colocalized with the αC-helix residues. The analysis of collective dynamics in the BRAF complexes with paradox inducers PLX4720 and Dabrafenib () showed a dominant peak corresponding to a highly flexible αC-β3-loop (residues 485–490) that controls structural preferences of the αC-helix. In these complexes, the central interface section forms a core of the extended hinge region that could promote functional movements of the αC-helix around the crystallographic inactive ‘out’ position. At the same time, the key residues in the secondary sections of the dimer interface (H510, K475, W476, H577, R558, R562, D565, and Y566) were more flexible. As a result, the collective movements of these BRAF dimers may involve concerted changes of the αC-helices around the central hinge and allow for proper alignment of the hydrophobic spines and the dimer interface. In the BRAF-PLX4720 complex, the second monomer adopts a more flexible DFG-out/αC-helix-out conformation, leading to larger movements in the αC-β3-loop (residues 485–490) and αC-helix (residues 491–505) (). A similar pattern of collective motions could be seen the BRAF-Dabrafenib complex (). Collective dynamics profiles in the BRAF complexes with Vemurafenib ( and PLX7904 () revealed that the paradox breaker could suppress movements of the αC-helix region, thereby compromising the proper alignment of the regulatory residues and interfacial residues in the global hinge. On the other hand, the ‘peripheral’ sections of the interface become more immobile (). Hence, though the overall shape of the extended hinge may be preserved in all BRAF dimer structures, paradox breaker could engineer changes in several key residues that would shift the major hinge center away from the central interface section. The central finding of this analysis was that paradox breaker may enforce global dynamic preferences of the inactive BRAF-WT monomer that would likely compromise coordinated allosteric changes in the αC-helix and dimer interface regions. This dynamic effect of PLX7904 may intervene with the effective side-to-side dimerization and therefore alleviate paradoxical activation.

The GNM-Based Essential Mobility Profiles of the BRAF Dimer Complexes with PLX4720 and Dabrafenib.

Collective dynamics of the BRAF dimers is analyzed using the cumulative contributions of the first five slowest GNM modes. Conformational mobility profiles and structural maps of collective motions are shown for the BRAF dimer complex with PLX4720 (pdb id 3C4C) (A, B) and Dabrafenib (pdb id 4XV2) (C, D). The slow mode shapes are annotated and colored in red (first monomer) and blue (second monomer). Structural maps of collective dynamics are based on fluctuations driven by the slowest five modes. The color gradient from red to blue indicates the increasing structural rigidity and refers to an average value over the backbone atoms in each residue.

The GNM-Based Essential Mobility Profiles of the BRAF Dimer Complexes with Vemurafenib and PLX7904.

Collective dynamics of the BRAF dimers is analyzed using the cumulative contributions of the first five slowest GNM modes. Conformational mobility profiles and structural maps of collective motions are shown for the BRAF dimer complex with Vemurafenib (pdb id 3OG7) (A,B) and PLX7904 (pdb id 4XV1) (C,D). The slow mode shapes are annotated and colored in red (first monomer) and blue (second monomer). Structural maps of collective dynamics are based on fluctuations driven by the slowest five modes.

Interestingly, a number of dimer interfacial residues (R509, H510, H577, R558, R562, D565, and Y566) are structurally rigid along slow modes in the BRAF-WT monomer. In particular, the side-to-side dimer architecture and central interface section are anchored by a critical residue R509 that is a rigid hinge site in the inactive monomeric structure (Figs ). Hence, structural stability of these sites may be determined by their global dynamics rather than local stabilization in a particular kinase conformation. These findings are consistent with the notion that global hinge centers that control collective dynamics may be intrinsically determined by the native protein topology and retain their mechanical role upon dimer formation and binding with inhibitors []. Differences in collective dynamics of the dimer structures may be associated with ligand-induced changes in the essential mobility of the αC-helix residues R506 and F516. These residues couple the αC-helix motions with the dimer interface interactions and are likely to be critical for the mechanism of paradoxical activation. Our results indicated that paradox inducer and paradox breaker inhibitors could selectively alter the essential mobility profile of several important sites in the global hinge region, while preserving the general shape of the slow modes and collective motions.

Negative Cooperativity and Dimer-Inducing Potential of the BRAF Inhibitors Are Drivers of Paradoxical Activation

Structural analysis of the BRAF dimer complexes showed a considerable similarity of the inhibitor binding modes in the first monomer, while binding of PLX4720 and Dabrafenib to the second binding site produced a different interaction pattern (). To identify the binding energy hotspots and quantify binding energetics in the first and secondary binding sites, we performed MM-GBSA calculations [] of the BRAF dimer complexes. In these simulations, we examined the molecular basis of negative cooperativity that may be induced by the BRAF inhibitors. MM-GBSA calculations of the BRAF dimer complexes were supplemented by a systematic alanine scanning of the binding site residues [] for each of the studied complexes. Our results reproduced the relative order of binding for the inhibitors (Figs ), particularly demonstrating a stronger binding affinity of Dabrafenib as compared to Vemurafenib and PLX4720. According to the experimental data, Dabrafenib binds to the BRAF kinase with IC50 ~6 nM [] and PLX4720 binds BRAF with IC50 = 160 nM [], while Vemurafenib showed weaker binding with IC50~260–360 nM [].

Binding Free Energy Calculations of the BRAF Dimer Complexes with PLX4720 and Dabrafenib.

Binding free energies and alanine scanning of the binding site residues in the first and second binding sites of the PLX4720-BRAF complex (A,B) and Dabrafenib-BRAF complex (C,D). The standard errors of binding free energy differences, which are the standard deviation of the mean values, were ~ 0.11–0.25 kcal/mol for the PLX4720-BRAF complex and 0.18–0.22 kcal/mol for the Dabrafenib-BRAF complex. The results of alanine scanning are shown in blue bars for the first binding site (A, C) and in red bars for the second binding site (B, D).

Binding Free Energy Calculations of the BRAF Dimer Complexes with Vemurafenib and PLX7904.

Binding free energies and alanine scanning of the binding site residues in the drug-bound monomer of the Vemurafenib-BRAF complex (A) and the PLX7904-BRAF complex (B). The standard errors of binding free energy differences, which are the standard deviation of the mean values, were ~ 0.15–0.25 kcal/mol for the Vemurafenib-BRAF complex and 0.2–0.3 kcal/mol for the PLX7904-BRAF complex. A close-up of the Vemurafenib binding mode (C) and PLX7904 binding mode (D). The inhibitor is shown in atom-colored sticks and annotated. The binding site residues are shown in orange sticks and annotated.

Computational alanine scanning identified common binding energy hotspots of all studied type I½ inhibitors: V471, V482, K483, L505, L514, F516, W531, F583, and F595 residues (Figs ). Some of these binding hotspots are located in the regulatory regions and selectivity pocket that are involved in collective dynamics of the BRAF dimers and allosterically coupled with the dimer interface residues. In our previous work we demonstrated that type I BRAF inhibitors may be less sensitive to conformational and variations of the R-spine residues, while more selective type II inhibitors could be strongly coupled with the dynamics of the hydrophobic spines []. According to alanine scanning performed in the current study, binding of type I½ inhibitors may be strongly influenced by the energetic hotspots from the catalytic salt bridge (K483), the αC-helix (L505), the αC-β4-loop (L514, 516), and the DFG motif (D594, F595). (Figs ). Some of these sites (L505, L514, F516) were especially critical for binding of paradox breaker PLX7904 as alanine modifications of these residues caused a substantial reduction in the binding affinity (). These residues link the αC-helix movements with the central section of the dimer interface (R506, R509, F516, and W450) and may be important for the inhibitor-induced dimerization and activation. Our results are consistent with mutagenesis studies that affirmed the critical role of these residues for the inhibitor binding []. According to these experiments, L505H substitution is a spontaneously occurring oncogenic BRAF mutation that can confer resistance to PLX4032 and Vemurafenib. Moreover, BRAF-L505H mutant can activate MAPK signaling, but has a smaller activating and oncogenic potential than the prevalent BRAF-V600E mutation. Paradox breaker PLX7904 is more sensitive to BRAF-L505H than PLX4720, but still insufficiently effective to overcome drug resistance conferred by this mutant []. Notably, only a new paradox breaker inhibitor PLX8394 can successfully combat mechanisms of drug resistance and inhibit the MAPK pathway in the BRAF mutant cells that are resistant to Vemurafenib and PLX7904 [].

We also computed binding affinities and carried out alanine scanning of PLX4720 and Dabrafenib interactions in the second binding site (). These results revealed a drastically reduced binding affinity of PLX4720 in the secondary site, where the inhibitor assumed a different binding mode and was bound to the DFG-out/αC-helix-out conformation (). As a result, even though PLX4720 occupies both monomers in the crystal structure, it may effectively inhibit only the first monomer. The observed negative cooperativity of PLX4720 may be associated with the increased conformational mobility of the second monomer, which is consistent with the reduced occupancy and faster off-rates of PLX4720 dissociation seen in the experiments []. The observed connection between the reduced inhibitor binding and the increased flexibility in the secondary monomer reflects the salient feature of negative cooperativity. Indeed, this effect is associated with allosteric mechanisms that occur without global structural changes but could trigger significant changes in conformational entropy []. Our results corroborated with the experimental studies showing that binding of some BRAF inhibitors may be often accompanied by negative cooperativity effects that promote paradoxical activation []. Dabrafenib binding induced similar structural environment of the binding sites and concerted motions of the two monomers, leading to small differences in the binding free energies between the first and second monomers (). Accordingly, this inhibitor may display an appreciable dimer-promoting potential but only a marginal negative cooperativity. Hence, dimerization and the induction of paradoxical activation by type I½ inhibitors may be associated with a different degree of negative cooperativity. Our results indicated that dimer-promoting potential and negative cooperativity of type I½ inhibitors may represent important molecular determinants of BRAF regulation and paradoxical activation.

Protein Sectors of Coevolutionary Residues Mediate Dimerization in the BRAF Complexes

Conformational dynamics and binding energetics of the BRAF dimer complexes suggested that paradoxical activation may be associated with ligand-induced modulation of allosteric interactions. To investigate this relationship, we proposed a model that incorporated both dynamic and coevolutionary residue correlations in the construction and analysis of the residue interaction networks. First, we explored coevolutionary residue dependencies to test whether key functional sites mediating allosteric interactions in the BRAF dimer complexes may display a strong coevolutionary signal. By using MI approach [] we evaluated the Kullback-Leibler conservation score () and the cumulative mutual information (cMI) score that estimates the degree of shared MI of a given residue with other protein residues (). Highly coevolving residues were primarily assembled in the ATP binding site, in the catalytic loop region and near the substrate binding motif (residues 619-WMAPE-623). Structural mapping of the high cMI residues revealed a well-connected network with a characteristic θ-like shape [] in which coevolutionary networks connected the ATP binding site with the substrate binding region through the hydrophobic spines (). In this network, the dimer interface residues R506 and R509 may be directly coupled with the R-spine residues L505 an F516. Hence, functional residues involved in collective motions and dimerization may be also coevolutionary coupled. We also found that coevolutionary residues may form three independent structural sectors: the binding site sector (Q530, W531, C532, F583, S536), the regulatory region sector (R506,F516, I572, H574, H568, F595) and the dimerization sector (D449, W450, E451, L505, R509, L514) ().

Coevolutionary Dependencies of the BRAF Kinase Residues.

Sequence conservation and coevolutionary propensities of the BRAF residues. (A) The KL conservation score. The residue numbering in the sequence conservation profile corresponds to the residue numbering in the BRAF crystal structures. (B) The cMI profile measures cumulative accumulation of mutual information per residue. The residue profiles are shown in blue bars. The inter-domain interface residues are shown in red circles. The R-spine residues are highlighted in green squares. (C) Structural mapping of high cMI residues (in blue spheres) onto the crystal structure of the Vemurafenib-BRAF complex (pdb id 3OG7). These residues included V471, A481, L515, F516, W531, L537, H574, L577, K578, S579, F583, L584, F595, W619, M620, D638, Y640, F642, and I644. The first monomer is shown in green ribbons and the second monomer is shown in cyan ribbons. (D) Structural mapping of the protein sectors of coevolving residues is shown in both monomers. The binding site sector residues Q530, W531, C532, F583, and S536 are shown in orange spheres; the regulatory sector residues R506, F516, I572, H574, H568, and F595 are shown in red spheres and the dimerization sector residues D449, W450, E451, L505, R509, L514, and F516 are shown in blue spheres.

These sectors of coevolving residues comprise contiguous clusters that only partially overlap and involve key residues from the ATP-binding site and regulatory αC-helix and αC-β4-loop regions. The topology and functional role of these modules are consistent with the main characteristics of protein sectors such as physical connectivity in the tertiary structure and biochemical independence in mediating protein function []. Strikingly, several key residues (L505, R509, L515, and F516) are strategically located at the cross-section of coevolving modules connecting the ATP binding site with the regulatory αC-helix and the dimer interface. These observations suggested that key regulatory residues connecting the αC-helix and the dimer interface may be subjected to coevolutionary constraints. Structural mapping of protein sectors also indicated that these modules of highly coevolving BRAF residues comprise a global interaction network that may promote allosteric interactions and facilitate conformational changes during dimerization.

We also computed proximity-based mutual information (pMI) score (), which is a structure-sensitive parameter defined as the average of cMI scores of all residues within a predefined distance from a given residue. Importantly, many regulatory sites displayed high pMI scores, suggesting that to maintain a particular function their local structural environment may be enriched by highly coevolving residues (). Interestingly, the dimer interface residues are located in structural proximity of high pMI sites. One group of the interfacial residues (D449, W450, R506, and R509) is linked to the high pMI residue F516. Another interfacial cluster (R558, R562, D565, and Y566) is located near the R-spine residues I572 and H574 that also showed high pMI score (). We argue that functional importance of the R-spine residues as drivers of regulatory changes may be determined by their coupling to clusters of flexible coevolving residues that undergo coordinated conformational changes upon dimerization and activation. This may highlight the functional role of the interfacial residues as transmitters of allosteric signals and conformational changes that are orchestrated by the regulatory spine residues.

Proximity-based Coevolutionary Profiles of the BRAF Dimer Complexes.

The ensemble-based pMI profiles of the BR AF complexes with PLX4720 (A), Dabrafenib (B), Vemurafenib (C), and PLX7904 (D). pMI values for each residue position are evaluated as the sum of cMI values of all residues within 5Å distance from a given residue. The distance between each pair of residues in the structure was calculated as the shortest distance between any two non-hydrogen atoms from respective two residues. pMI profiles are computed using average values obtained from MD trajectories and ensemble-based definition of the local residue environment. The residue profiles are shown in blue bars. The inter-domain interface residues are shown in red circles. The R-spine residues are highlighted in green squares.

BRAF Inhibitors Modulate Centrality and Communication Propensity of Allosteric Centers in the Residue Interaction Networks

A novel methodological aspect of this work was an integration of coevolutionary residue dependencies and dynamic maps of residue correlations in the construction and analysis of the residue interaction networks. In this model, residues with high level of centrality correspond to key mediating centers that connect coevolutionary and dynamically coupled residues in a global interaction network. By employing this protein structure network model, we determined how paradox inducer and paradox breaker inhibitors could modulate organization of the interaction networks and activation pathways. Different network parameters such as residue-based centrality and communication propensity were explored to characterize how BRAF inhibitors could affect mediating capabilities of kinase residues and the distribution of allosteric centers. First, we compared the residue centrality profiles of the BRAF-WT monomer structure and BRAF dimer complexes (). The residue centrality of the BRAF-WT monomer revealed several conserved peaks corresponding to the αC-helix region (residues 491–516), the 574-HRD-576 catalytic motif, the 593-DFG-595 motif, and the organizing αF-helix (residue 635–651). The integrating role of the αC-helix and the αF-helix regions is known to be critical for modulation of kinase dynamics and activity []. According to our results, paradox inducer inhibitors could alter centrality of the kinase residues by amplifying the peak near the αC-helix/αC-β4-loop region and central section of the dimer interface (R506, R509, F516, and W450) (). In some contrast, paradox breaker PLX7904 may induce important network changes by reducing centrality of the regulatory αC-helix region to the level observed in the BRAF-WT monomer (). The major peak of the distribution shifted towards the secondary interfacial cluster formed by C-terminal residues (R558, R562, D565, and Y566). As a result, network preferences of functional centers may be partly distorted in the PLX7904-BRAF complex due to stabilization of a nonproductive αC-helix-out conformation. The impaired allosteric coupling between the αC-helix and the central section of the dimer interface may force the system to “dispatch allosteric traffic” through longer and arguably less efficient paths proceeding through the secondary interfacial cluster. These findings are consistent with the analysis of collective dynamics, showing that paradox breaker may induce dynamic and network signatures that are characteristic of the inactive BRAF-WT monomer.

Residue-Based Centrality Profiles in the BRAF Structures.

The residue centrality profiles of the BRAF dimer complexes with PLX4720 (pdb id 3C4C) (A), Dabrafenib (pdb id 4XV2) (B), Vemurafenib (pdb id 3OG7) (C), and PLX7904 (pdb id 4XV1) (D). The distribution for the BRAF-WT (in green) is shown for comparison on all panels. The centrality profiles for the first monomer are shown in red lines with filled red squares and for the second monomer in blue lines with filled blue circles. The centralities of the R-spine residues in the drug-bound monomer are highlighted as filled maroon squares.

By expanding a force constant model of protein dynamics [], we computed communication propensities (CP) of protein residues defined as average variations of the effective “distance” metric that measures both distance fluctuations and variations in the pMI scores between a given residue and all other residues. In this model, coevolving and dynamically correlated residues whose effective distances fluctuate with low or moderate intensity are expected to communicate with the higher efficiency than the residues that experience large fluctuations. We found that paradox inducers showed a similar CP profile in which the largest dominant peak corresponded to the αC-helix/ αC-β4-loop region showing high values for the interfacial residues R506, R509, H510, and F516 (). In addition to the interfacial residues, the αC-helix and R-spine residues may be also characterized by high communication capabilities. Consequently, by enhancing global mediating capacity of these residues, paradox inducers would likely direct the shortest inter-residue pathways through the central interface section, which could promote the effective side-to-side dimerization and activation. In the BRAF-PLX7904 complex, the CP values of the αC-helix and central interfacial section were markedly lowered (). Similarly to the residue centrality profile, the major peak shifted to the secondary interfacial region. Hence, paradox breaker may trigger dislocation of mediating centers in the central hinge region that would weaken dimer formation and favor alternative activation routes. Structural mapping of the interfacial regions illustrated differences in the communication preferences of the interface regions that may affect allosteric signaling in the BRAF dimers (). To conclude, our analysis of the residue interactions networks demonstrated that paradox breakers may circumvent paradoxical activation by closely mimicking structural, dynamic and network properties of the inactive BRAF-WT monomer.

Allosteric Communication Propensities in the BRAF Structures.

The residue-based CP profiles of the BRAF complexes with PLX4720 (A), Dabrafenib (B), Vemurafenib (C) and PLX7904 (D). The profiles for the first monomer are shown in red lines with filled red squares and for the second monomer in blue lines with filled blue circles. The R-spine residues are highlighted as filled maroon squares.

Paradox Inducer and Paradox Breaker Inhibitors Induce Specific Allosteric Communication Pathways

Using matrix of communication distances between protein residues in the BRAF dimer complexes, we obtained ensembles of shortest paths connecting each pair of residue nodes in the BRAF structures. In this model, allosteric communication pathways proceed via a network of coevolving and dynamically coupled residues. By simulating ensembles of short paths that connect the ligand binding sites in the BRAF dimers, we explored whether paradox inducer and breaker inhibitors may activate specific communication pathways and exploit different sections of the dimer interface. We reported atomistic details of the most probable paths that connect the regulatory DFG motifs in the first and second monomers. First, we found that single, high occupancy communication routes may dominate the ensemble in the BRAF dimer complexes with paradox inducers PLX4720 () and Dabrafenib (). The most probable path (82% occupancy) in the PLX4720-BRAF dimer connected F595 (monomer A) with I513, V511, and N512 residues from the αC-β4-loop before reaching the key interfacial residue R509 of the monomer A. At this point, the path bridged two monomers by connecting R509 (monomer A) with F516 of the monomer B. The path then navigated through the monomer B, moving from F516 to the binding site residues L514, Q530, W531, C532, and F583 residues (). Interestingly, this path traversed through most of the binding energy hotspots in the second monomer. In this BRAF complex, a critical inter-monomer juncture linked R509 (monomer A) with F516 (monomer B) by exploiting the central section of the inter-monomer interface (R506, R509, F516, W450).

Structural Map of Allosteric Communication Pathways in the BRAF Dimer Complexes.

The most probable communication pathways in the PLX4720-BRAF complex (pdb id 3C4C) (A), Dabrafenib-BRAF complex (pdb id 4XV2) (B), Vemurafenib-BRAF complex (pdb id 3OG7) (C), and PLX7904-BRAF complex (pdb id 4XV1) (D). The BRAF structures are shown in ribbons, with the monomer A in orange and monomer B in cyan. The key residues along the communication routes are annotated and shown in green spheres. The error bars on the pathway occupancies are within 5%.

Notably, the observed communication pathways navigated through protein sectors formed by the coevolutionary residues. These findings suggested that coevolutionary networks may form direct pathways for transmitting allosteric signals []. The most probable path in the Dabrafenib-BRAF complex (85% occupancy) could proceed from F595 of monomer A to the R-spine residues of the αC-helix (V504, L505), subsequently reaching to the interfacial residue R509 of the monomer A (). At this point, the path bridged over the inter-domain interface by connecting to R506 and R509 of the monomer B. After crossing the dimer interface, the route proceeded to Y566 and H510 residues in the monomer B before reaching I513 of the αC-β4 region and arriving at the final destination F595 of the second monomer (). In this case, the most probable path utilized direct coupling of the R-spine residues (V504, L505) with the central interfacial sites R506 and R509. Hence, paradox inducers may facilitate rapid and efficient communication through a single pathway connecting the binding site with the central interfacial section (.

Of special interest was a comparison of communication pathways in the BRAF complexes with Vemurafenib and PLX7904, whose structural and dynamic attributes are very similar. In these structures, signaling pathways become more diverse with 45% occupancy for the most optimal path in the Vemurafenib-BRAF complex and 55% occupancy in the PLX7904-BRAF complex (). The favorable route in the Vemurafenib-BRAF complex proceeded from F595 of the drug-bound monomer to I592 and the αC-β4-loop residues (N512, I513) before reaching H510 residue. At this point, the path crossed the interface to the drug-free monomer, connecting with M517, F516 and the αC-helix L505 residue (). The path then linked with T508 before reaching F595 of the second monomer. Interestingly, this communication route proceeded through the interfacial residue H510, while the central section of the interface (R506, R509) was avoided. In the BRAF complex with a paradox breaker PLX7904, the ensemble of communication pathways was also fairly diverse and none of the dominant routes exploited the central interface section. The first optimal path (55% occupancy) was longer and less optimal in comparison with favorable routes found in other BRAF dimers (). Starting from F595 residue of the drug-bound monomer, the path preceded to the R-spine residues V504 and L505. The route then followed to K507 before making a turn and reaching out to Y566 and R562 residues. At this point, the path crossed the dimer interface and connected to D479 and K475 residues of the second monomer. Subsequently, this path linked K475 with I457, A526, I527 residues from the N-terminal lobe in the second monomer before reaching to E501 (from the catalytic salt bridge) and connecting to F595 in the second monomer. Hence, the most favorable signaling route in the PLX7904-BRAF complex may bypass the central interface section and navigate through other interfacial residues.

The central finding of this analysis was the emergence of suboptimal and less efficient communication pathways for the paradox breaker as opposed to direct and more efficient transmission routes detected in the BRAF complexes with the paradox inducers.

Rapid and efficient communication in the BRAF complexes with paradox inducers may come at the expense of high dependency on critical modes in the αC-helix, αC-β4-loop, and dimer interface region (R506, R509, and F516). Of particular significance are key residues of the inter-monomer hinge F516 and R509 that control proper alignment of hydrophobic spines and allosteric coupling between the αC-helix and dimer interface. These high centrality residues that minimize the short path length and ensure the efficiency of allosteric communications are experimentally known to be important for kinase activity and regulation. As a result, even minor modifications of these residues could alter activation signaling pathways and trigger drug resistance. A different organization of optimal pathways in the BRAF-PLX7904 complex may reduce the efficiency of allosteric communications and impair activation signaling, while making drug binding more resilient to mutations of allosteric hotspot nodes. In network terms, the increased diversity of communication pathways would likely diminish network efficiency, but may potentially increase level of tolerance against mutations of mediating residues []. By linking paradoxical activation with the ligand-mediated changes in allosteric networks and communication pathways, our results may provide a network-based rationale of the distinct functional effects exerted by the type I½ BRAF inhibitors.

Conclusions

In this work, we explored molecular mechanisms of dimerization-induced BRAF regulation by a panel of type I½ inhibitors. Our study was motivated by recent revelations that most of the BRAF inhibitors can paradoxically activate kinase by inducing dimerization and promoting dimer-dependent transactivation. These experimental discoveries have shown that therapeutic and clinical utility of BRAF inhibitors that exert paradoxical activation can be diminished, creating a need for a quantitative assessment of dimerization potential induced by the inhibitor. By combining atomistic and coarse-grained simulations with modeling of the residue interaction networks and coevolution, we have identified molecular determinants of differential functional effects produced by a series of paradox inducer and breaker inhibitors. According to our results, dimer-promoting potential and negative cooperativity of type I½ inhibitors may be associated with the induction of paradoxical activation. These functional effects may be controlled by several regulatory sites in the extended hinge formed by the αC-helix and dimer interface residues. The important result of this study is that paradox inducers and paradox breaker may have specific dynamic and network signatures that uniquely characterize their distinct functional effect on BRAF activation. The central finding of our investigation is that paradox breaker PLX7904 may evade activation by closely mimicking structural, dynamic and network properties of the inactive BRAF-WT monomer. We argue that BRAF inhibitors that amplify structural features of the inactive BRAF-WT monomer and suppress dimer-promoting function may alleviate the extent of paradoxical activation and enhance drug efficacy. By simulating allosteric interaction networks and communication pathways in the BRAF structures, we determined that dimerization may promote efficient and robust activation pathways that are enabled through coordinated structural changes of the αC-helix and dimer interface residues. Furthermore, paradox inducer and breaker inhibitors may activate specific signaling routes that correlate with the extent of paradoxical activation. The results suggested that paradox inducer inhibitors may facilitate a rapid and efficient communication by a single pathway, while a paradox breaker may induce a broader ensemble of suboptimal and less efficient communication routes.

Understanding of system-based relationships between network efficiency and drug binding may prove useful in the design of novel BRAF kinase inhibitors. The lessons from our study could also inform discovery efforts aiming in mitigating adverse biological consequences and tailoring existing BRAF inhibitors into paradox breakers. According to the current view, design of type I½ BRAF inhibitors capable of stabilizing a non-productive αC-helix-out conformation may help to block the dimerization interface. Our results suggested that these inhibitor-induced changes may result in the induction of alternative signaling routes and could only partially alleviate dimerization-induced transactivation. An alternative design strategy may involve direct targeting of several binding energy hotspots (L505, R509, L515, and F516) that also serve as critical hubs in the allosteric interaction networks, coevolutionary networks and signaling pathways. By exploiting the unique functional role of these residues and simultaneously enforcing structural features of the inactive BRAF monomer, next generation of sulfonamide inhibitors may have a potential to act as paradox breakers in crystal and cellular environments.

Materials and Methods

Atomistic and Coarse-Grained Simulations

MD simulations were performed for the monomeric structure of the BRAF-WT, and BRAF dimer-inhibitor complexes that included BRAF-WT bound with PLX4720 (pdb id 3C4C), BRAF-WT complex with Dabrafenib (pdb id 5CSW), BRAF-V600E bound with Dabrafenib (pdb id 4XV2), BRAF-V600E complexed with Vemurafenib (pdb id 3OG7), and BRAF-V600E complex with PLX7904 (pdb id 4XV1). A combination of several long MD simulations and multiple shorter MD runs was undertaken to ensure a broader and more unbiased sampling. We have carried out two independent 500 ns and five independent 200 ns MD for each of the studied BRAF structures. The equilibrium ensembles and average properties of the BRAF dimer complexes were obtained by aggregating simulations from all independent MD trajectories. All crystal structures were obtained from the Protein Data Bank (RCSB PDB www.rcsb.org) []. Structure preparation process included several important steps that were previously reported in detail in our studies of protein kinases [] and molecular chaperones []. In this protocol, hydrogen atoms and missing residues were assigned using the WHATIF program []. The missing segments of the A-loop in the BRAF crystal structures were modeled and reconstructed with the ArchPRED program []. The initial protonation states were assigned using the WHATIF program and optimized with the aid of the H++ web server []. The ligand charges and parameters were derived using a standard protocol that was described in detail in our recent studies [77, ]. According to this procedure, the crystallographic inhibitor conformations were initially optimized at the B3LYP/6-31G level using the Gaussian 09 package []. The charge parameters for the BRAF inhibitors were assigned initially by the ParamChem web server [ and subsequently refined through the restrained electrostatic potential charge fitting procedure (RESP) []. The force field parameters for the BRAF inhibitors were obtained using a VMD plugin ffTK []. MD simulations were carried out using NAMD 2.6 package [] with the CHARMM22 force field [] and the explicit TIP3P water model []. The employed MD protocol was previously reported in detail in our studies of protein kinases [] and molecular chaperones []. The initial structures were solvated in a water box with the buffering distance of 10 Å. Long-range nonbonded van der Waals interactions were computed using an atom-based cutoff of 12Å with switching van der Waals potential beginning at 10Å. The smooth particle mesh Ewald (PME) method [] was employed to treat the long-range electrostatics. The simulation steps preceding the production period are consistent with our recent studies and were detailed previously []. An NPT production simulation was run on the equilibrated structures for 500 ns (or 200 ns for shorter simulations) keeping the temperature at 300 K and constant pressure (1 atm) using Langevin piston coupling algorithm.

The elastic network-based GNM approach [] was used to describe collective motions and functional dynamics of the BRAF dimers. The formulation and implementation of this approach was previusly detailed in our recent studies []. In a nutshell, the protein structure is described as an isotropic fluctuating network of N residue nodes, each represented by the respective alpha-carbon atom. The protein topology is described by N × N Kirchhoff matrix of inter-residue contacts, Γ. The diagonal elements represent the coordination number of each residue. The pairs of residues located within an interaction cutoff distance r = 7.0 Å are connected by springs with a uniform spring constant γ. The GNM modes are determined by diagonalization of the Kirchhoff matrix, Γ = UΛU.

Here, U is a unitary matrix, U = U−1 of the eigenvectors u of Γ and Λ is the diagonal matrix of eigenvalues λ. The elements of the kth eigenvector u describe the displacements of the residues along the kth mode coordinate, and the kth eigenvalue, λ, scales with the frequency of the kth mode, where 1 ≤ k ≤ N– 1. KB is the Boltzmann constant, T is the absolute temperature. The mean square fluctuations of a residue are evaluated as a sum over the contributions of all modes. The fluctuation of the ith atomic degree of freedom along the eigenvector u reflects the mobility of residue i in the kth mode.

Conformational mobility profiles in the essential space of low frequency modes were obtained using the oGNM web server [].

Binding Free Energy Analysis

The binding free energy was calculated using the MM-GBSA approach []. In this model, the free energy of binding between an inhibitor and the BRAF kinase was evaluated as follows:

The gas-phase contribution <ΔG> to the binding free energy is the difference in the molecular mechanics energy of the complex and the isolated protein and ligand. This contribution is the sum of the differences in the intramolecular energies ΔEintra, the van der Waals interaction energy ΔE, and the electrostatic interaction energy ΔE: where E is the energy of the bonded terms (bonds, angles, dihedral angles, and improper angles) of a given molecule; E is the van der Waals energy of the molecule; and E is the electrostatic energy of the molecule. ΔG is the solvation free energy, and −TΔS is an entropy term.

The solvation free energy ΔG is the difference between the solvation energy of the complex and solvation free energies of the isolated protein and ligand. This term is computed as the sum of nonpolar and polar terms:

The nonpolar contribution is evaluated from the solvent accessible surface area (SASA) model as +β where σ = 0.00542 kcal/ (mol*Ǻ2) and β = 0.92 kcal/mol. The electrostatic contribution to the solvation free energy of the system was estimated using the analytical generalized Born (GB) model. The entropy term includes contributions of translational ΔS, rotational ΔS and vibrational ΔS components:

The vibrational entropy terms were computed using normal mode analysis in the VIBRAN module of the CHARMM program. All energy terms are calculated using MD trajectories of the complexes, which is followed by separation of the complexes into isolated protein and ligand structures and subsequent minimization of the isolated molecules. Computational alanine scanning [] was performed by mutating the binding site residues and computing the binding free energy for the mutated system using 10,000 snapshots from 500 ns MD trajectory of the original complex. A single trajectory variation of the MM-GBSA approach was employed in which MD simulations of the BRAF complexes were followed by separation of the individual snapshots into isolated protein and ligand conformations. This protocol has been proven to be robust and efficient in binding free energy computations of the ligand-protein complexes, avoiding uncertainties and reducing cancellation errors associated with the differences in the intramolecular energies [].

Mutual Information Model and Coevolutionary Analysis

Coevolutionary dependencies of BRAF residues were evaluated using MI model and MISTIC approach []. The details of the MI model employed in this study were reported in our recent study of protein kinases []. The following MI parameters were computed: the Kullback-Leibler conservation score, a cumulative mutual information score (cMI), and a proximity mutual information score (pMI). The Kullback-Leibler (KL) score KLConsScore measures sequence conservation in the protein kinase family. For each column of the MSA, the KL conservation is calculated according to the following formula:

Here, P(i) is the frequency of amino acid i in that position and Q(i) is the background frequency of the amino acid in nature calculated using an amino acids background frequency distribution obtained from the UniProt database. cMI is a sequence-based MI score per residue position that estimates the contribution of a given residue to the MI network. cMI is calculated as the sum of MI values above a threshold t = 6.5 for every pair in which a particular residue of interest appears [].

pMI score is a structure-based MI metric that is defined as the average of cMI scores of all the residues within a given distance from a given residue in the protein structure. The distance between each pair of residues was calculated as the shortest distance between any two heavy atoms that belong to each of these two positions. The threshold t = 5 Å is used to define structural proximity of a residue:

pMI scores of protein residues in the BRAF dime complexes were computed using ensemble-based definition of local residue proximity. In this dynamics-based model, the amount of mutual information shared by a given residue with the spatially close residues was averaged over 10,000 conformations from the MD trajectories of the BRAF dimers.

Modeling of the Residue Interaction Networks: Communication Residue Propensities and Allosteric Pathways

In the construction of a protein structure network, a graph-based model of protein topological connectivity is used where residues are considered to be network nodes and are connected by edges representing residue interactions. The details of the graph construction using a particular interaction cut-off strength (Imin) were extensively discussed in the initial reports [] and our previous studies []. We expanded the original model as the network edges that define residue connectivity are weighted based on both cross-correlation dynamic fluctuation and coevolutionary dependencies measured by the MI scores. The weight w of a particular edge is defined as the element of a matrix r(x,x) measuring the dynamic and coevolutionary correlations for each pair of residues []. The length of the edge that connects nodes i and j is calculated from the corresponding generalized correlation coefficient as w = −log[r(x,x)]. The weighted graph model defines a residue interaction network that favors a global flow of information through edges between residues associated with dynamics correlations and coevolutionary dependencies, which may be important for allosteric communication in the protein structure. The ensemble of shortest paths is determined from matrix of communication distances by the Floyd-Warshall algorithm [] that compares all possible paths between each pair of residue nodes.

Using the constructed protein structure network, we compute the residue-based centrality parameter. The centrality of residue i is determined as its network betweenness computed as a fraction of the shortest paths between all pairs of residues that pass through residue i: where g denotes the number of shortest paths connecting j and k, and g(i) is the number of shortest paths between residues j and k that navigate through the node n. Residues that populate a significant portion of the shortest paths connecting all residue pairs are characterized by high betweenness values (high residue centrality). For each node n, the betweenness value can be normalized by the number of node pairs that exclude node n which is given as (N - 1)(N - 2) / 2, where N is the total number of nodes in the connected component that node n belongs to. The normalized betweenness of residue i can be expressed as follows: g is the number of shortest paths between residues j and k; g(i) is the fraction of these shortest paths that pass through residue i.

We also computed communication propensities (CP) of protein residues defined as average variations of the effective “distance” metric that measures both distance fluctuations and variations in the pMI scores between a given residue and all other residues. Small fluctuations and respectively large CP values for a given residue are associated with high communication propensities and signal a high level of dynamic and coevolutionary cooperativity. For each residue, CP metric is evaluated as follows: where d is the instantaneous distance between residue i and residue j, ΔpM is the respective difference between pMI scores of residues i and j at the current 3D positions x and x; k is the Boltzmann constant, T = 300K. 〈 〉 denotes an average taken over the MD simulation trajectory and d = 〈d〉 is the average distance from residue i to all other residues in the protein (the sum over j* implies the exclusion of the atoms that belong to the residue i). ΔpM is the average difference in pM scores.

In the expression (13), w1 and w2 are weighting factors that are adjusted to optimize efficiency of the residue interaction network. The network efficiency is related to the average length of shortest paths among residue nodes. The length of a path L(n,n) between distant nodes n and n is the sum of the edge weights between the consecutive nodes (n,n) along the path:

Collective Dynamics of the BRAF-WT Monomer Structure.

(A). The GMM-based conformational mobility profiles of the BRAF-WT monomer structure (pdb id 4WO5). The essential mobility of kinase residues is shown in the space of first two slowest modes (in green), the first three slowest modes (in red) and the first ten slowest modes (in blue). The position of the dimer interface residues in these profiles are highlighted in filled maroon-colored circles. (B) The inter-residue distance fluctuations map in the BRAF-WT is based on residue movements along the first three slow modes and derived using the ANM web server []. The large inter-residue fluctuations are shown in blue and small fluctuations are shown in red. (C) Structural maps of collective dynamics are based on fluctuations driven by the slowest three modes. The color gradient from red to blue indicates the increasing structural rigidity and refers to an average value over the backbone atoms in each residue.

(TIF)

Click here for additional data file.

The Slow Modes Shapes of the BRAF Dimer Structures.

The GMM-based conformational mobility profiles along the three slowest modes are shown for the BRAF dimer complex with PLX4720 (pdb id 3C4C) (A), Dabrafenib (pdb id 4XV2) (B), Vemurafenib (pdb id 3OG7) (C), and PLX7904 (pdb id 4XV1) (D). The slow modes shapes are colored in green (first monomer) and red (second monomer).

(TIF)

Click here for additional data file.

Collectivity of the Slow Modes in the BRAF Dimer Structures.

(A) The degree of collectivity of first 20 slow modes in the BRAF dimer structures is shown for the PLX4720-BRAF complex (in blue), Dabrafenib-BRAF complex (in red), Vemurafenib-BRAF complex (in green) and PLX7904-BRAF complex (in orange). High collectivity of the first five modes indicates that conformational movements in the BRAF structures can be described using a cumulative contribution of these slow modes. (B) Cumulative overlap of the GNM slow modes with the principal components describing conformational changes in the ensemble of BRAF crystal structures. This ensemble included all known BRAF dimer complexes with different classes of inhibitors as described in Fig 1. The softest five modes described ~90% of structural variations observed in the BRAF crystal structures.

Fig 1

Structural Landscape of the BRAF Kinase Dimer Complexes with Small Molecule Inhibitors.

(A) Structural alignment of the crystal structures of the BRAF kinase dimers with type I inhibitors (DFG-in/αC-in kinase conformation) are shown in green (monomer A) and cyan (monomer B) and included pdb entries 2FB8, 3D4Q, 3PSB, 3PRF, 3PRI, 3PPK, 3PPJ, 3Q4C, 3PSD, 4E26, 4H58, 4MNF, and 3OMV. (B) The structures of BRAF complexes with type II inhibitors (DFG-out/αC-in kinase conformation) included pdb entries 1UWH, 1UWJ, 5V9C, 5CT7, 4KSP, 4KSQ, 4FC0, 4G9R, 4G9C, 4DBN, 3Q96, 3II5, 3IDP, and 4JVG. (C) The crystal structures of BRAF complexes with type I½ inhibitors (DFG-in/αC-out kinase conformation) included pdb entries 3C4C, 3OG7, 4FK3, 4EHG, 3SKC, 3TV4, 3TV6, 4E4X, 4EHE, 4MBJ, 4PP7, 4CQE, 4XV1, 4XV2, 4XV3, 4XV9. In this class of BRAF complexes, some crystal structures (pdb entries 3SKC, 3TV4, 3TV6, 4E4X, 4EHE, 4MBJ, 4PP7) have a small helical motif in the activation segment. (D) The superposition of the crystal structures of BRAF dime complexes with PLX4720 (pdb id 3C4C) (in green), Vemurafenib (pdb id 3OG7) (in red), PLX7904 (pdb id 4XV1) (in blue), and Dabrafenib (pdb id 4XV2, 5CSW) (in cyan). The regulatory regions are annotated and structural arrangements of the αC-helix and the DFG motif in BRAF complexes are highlighted.

(TIF)

Click here for additional data file.

Structural Mapping of the Dimer Interface Regions and R-spine Residues.

(A) Structural mapping of the central interface section. The first monomer is shown in green ribbons and the second monomer is shown in cyan ribbons. The residues from the central section of the dimer interface (R506, R509, and F516) are annotated and shown as green spheres (first monomer) and cyan spheres (second monomer). The R-spine residues (V504, L505, F516, I572, H574, and F595) are annotated and shown as red spheres (first monomer) and blue spheres (second monomer). (B) Structural mapping of the secondary interface regions. The interface residues in this region (K475, W476, H577, R558, R562, D565, and Y566) are annotated and shown as green spheres (first monomer) and cyan spheres (second monomer).

(TIF)

Click here for additional data file.