Modeling of RAS complexes supports roles in cancer for less studied partners

Background RAS protein interactions have predominantly been studied in the context of the RAF and PI3kinase oncogenic pathways. Structural modeling and X-ray crystallography have demonstrated that RAS isoforms bind to canonical downstream effector proteins in these pathways using the highly conserved switch I and II regions. Other non-canonical RAS protein interactions have been experimentally identified, however it is not clear whether these proteins also interact with RAS via the switch regions. Results To address this question we constructed a RAS isoform-specific protein-protein interaction network and predicted 3D complexes involving RAS isoforms and interaction partners to identify the most probable interaction interfaces. The resulting models correctly captured the binding interfaces for well-studied effectors, and additionally implicated residues in the allosteric and hyper-variable regions of RAS proteins as the predominant binding site for non-canonical effectors. Several partners binding to this new interface (SRC, LGALS1, RABGEF1, CALM and RARRES3) have been implicated as important regulators of oncogenic RAS signaling. We further used these models to investigate competitive binding and multi-protein complexes compatible with RAS surface occupancy and the putative effects of somatic mutations on RAS protein interactions. Conclusions We discuss our findings in the context of RAS localization to the plasma membrane versus within the cytoplasm and provide a list of RAS protein interactions with possible cancer-related consequences, which could help guide future therapeutic strategies to target RAS proteins. Electronic supplementary material The online version of this article (doi:10.1186/s13628-017-0037-6) contains supplementary material, which is available to authorized users.


Background
Tumor exome sequencing studies have uncovered RAS mutations in over 30% of solid tumors, yet to date no therapy has been discovered to effectively treat RAS driven cancers. RAS proteins themselves provide an attractive therapeutic target but efforts aimed at drugging RAS directly thus far have failed [1]. More complete characterization of RAS signaling may provide badly needed insights to support renewed efforts aimed at developing RAS-targeted therapies [1,2]. In particular, little is known about how RAS proteins interact with less-studied interaction partners, and how these proteins contribute to tumorigenesis.
RAS is a small GTPase that localizes to the plasma membrane (PM) and upon receiving extracellular stimuli, initiates multiple signaling cascades inside the cell [3]. There are 4 major isoforms of RAS, including KRAS4A, KRAS4B, HRAS and NRAS. All RAS isoforms transition between active and inactive states, a process that is regulated by guanine nucleotide exchange factors (GEFs) and GTPase activating proteins (GAPs). The isoforms are highly similar in amino acid sequence; all isoforms are 188-189 amino acids in length, with nearly 85% sequence identity [4]. The first 87 residues are 100% homologous and most of the sequence differences occur in the short C-terminal hypervariable region (HVR) (residues 166 to 189), which is only 10% conserved among RAS proteins (Fig. 1).
RAS signaling is deceptively complex. In contrast to their sequence similarity, different RAS isoforms have been associated with distinct functionality and show distinct patterns of expression in various cell lines and tissues [5,6]. Recent work has shown that RAS isoforms also form distinct homo- [7] and hetero-dimers [8]. In cancer, the clinical phenotype of patients with RAS mutations is highly context dependent. The prevalence of specific mutant RAS alleles differs among tumor types, as does the predominantly mutated RAS isoform. For example, head and neck [9] tumors carry mutant HRAS, melanomas [10] carry mutant NRAS and pancreatic tumors [11] carry mutant KRAS. For KRAS mutant driven tumors, non-small cell lung cancers tend to carry G12C mutations [12], while the G12D [13] substitution is predominantly found in pancreatic cancer. This complexity suggests that no single therapy will be able to universally target mutant RAS.
Recently, Chavan et al. [14] described a mechanism by which G12 mutations render RAS isoforms more active. According to this model, the G12 residue binds to the HVR, keeping RAS in a less active state by blocking interaction of RAS with effectors. Mutations at this site may impair HVR binding, increasing RAS interaction with downstream effectors. Nonetheless, differences in the prevalence of G12 amino acid mutants across tumor types cannot be explained by the weakening of the HVR binding by G12 mutations alone. If that were the case, the mutant amino acids should occur with the same frequency across tumors. One possible explanation for the difference in mutant prevalence is that G12 mutations also affect RAS binding with protein interaction partners, the expression of which may differ by tumor type [15]. Thus better understanding of RAS protein interactions will be essential for understanding its role in tumorigenesis.
From a therapeutic perspective, RAS proteins have been perceived as "undruggable" due to a lack of effective inhibitors, even after more than three decades of research [16]. However, very recently four groups have reported novel allosteric and orthosteric inhibitors targeting RAS protein protein interactions (PPI) by taking advantage of a pocket that becomes accessible upon compound binding [17][18][19][20]. In 2016, Athuluri-Divakar et al. [21] identified a high affinity small molecule inhibitor that binds to RAS effectors such as RAF, Ral-GDS and PI3Ks and disrupts their interactions with RAS. Subsequently, RAS inactivation through PPI disruption has become the leading approach for targeting RAS [22], further motivating studies of the RAS protein interaction network.
RAS PPIs have predominantly been studied in the context of the RAF and PI3kinase oncogenic pathways [23] and structural analyses have concluded that RAS isoforms bind to canonical downstream effectors and upstream regulators in these pathways using the highly conserved switch I (residues between 30 and 40) and II regions (residues between 60 and 76, [24]). Understanding of RAS interactions is largely limited to the interactors that have RAS binding related domains (RBrDs) (see Methods), here referred to as "canonical interactors". Other, non-canonical RAS interactors have been experimentally identified, including FNTA, FNTB and chaperone proteins PDE6D and LGALS1, but it remains poorly understood how these proteins interact with RAS.
To study RAS isoform specific PPIs, including partners with no well-defined RBrD, we predicted which residues The G12, G13, A59 and Q61 hotspot mutations are illustrated in red. The HVR region is partially (residues 166 to 171) present in the PDB and the figure on the surface of RAS interact with each binding partner. The resulting map of binding to RAS was used to model possible RAS complexes and identify possible competitive binding relationships. We then analyzed the putative effects of mutant alleles on RAS interactions with binding partners. Our analysis suggests a novel binding interface on the allosteric region/HVR of RAS proteins that is accessible when RAS is not localized to the plasma membrane. Interestingly, multiple interaction partners that bind via this new interface are putative cancer genes.

Results
In order to identify interaction interfaces on RAS proteins, we first curated an experimentally verified RAS isoform specific PPI network that consists of as many as 103 RAS protein interactions (see Methods). Of 77 RAS interacting proteins in our network, 45 proteins do not have an RBrD (Fig. 2). We first curated interface residues from seven RAS co-complex structures available in the Protein Data Bank [25] (PDB; 6 HRAS, 1 KRASs and no NRAS complexes), all of which capture canonical RAS interactors.
We next built RAS 3D complex models using PRISM [26], a template-based protein complex prediction server that has been extensively benchmarked and reports high accuracy for identifying interfaces [27]. PRISM predicted binding interfaces for 20 complexes involving the RAS proteins and 17 binding partners. PRISM was not able to identify binding interfaces for the other 37 partners with available structural data including ITSN1, RGS12, GRB1 and UBC. We used the 20 models of RAS interactions to determine which amino acids on the surface of RAS participate in each interaction (see Methods). This resulted in a map of interface residues on RAS proteins ( Fig. 3) enabling us to analyze patterns of interface residue usage across interaction partners (Fig. 4). This map revealed two categories of RAS interaction partners: proteins that bind via the switch regions, and proteins that bind via the allosteric region and/or the HVR. The HVR does not adopt a stable conformation and therefore is challenging to capture via X-ray crystallography. Given that the available PDB structures in general do not include the full HVR, PRISM predictions for binding in the allosteric/HVR region are best interpreted as "interactors that do not bind to the switch region".

RAS binding partners that do not bind to the switch region
Most studies of RAS signaling in cancer focus on effectors that bind the switch regions, with few studies addressing binding partners without RBrDs. Among these, one of the earlier studies addresses the involvement of the HVR in the KRAS -CALM interaction [6]. This result has only recently been revisited [2,14] after the discovery of RAS homodimerization [28]. It has previously been demonstrated that residues mediating RAS dimerization, R135, R161 and R164 located in the 4th and 5th helices, are not in the effector region. These helices are conserved across all three of the RAS proteins [29]. A different line of evidence for the non-switch regions' involvement in RAS interactions is present in the PDB complex between PDE6D and RHEB (PDB ID: 3T5G). RHEB is a GTPase with 80% sequence similarity to RAS proteins. RHEB primarily localizes to the ER and Golgi Apparatus [30] and according to a co-crystal structure (PDB ID: 3T5G) it interacts with chaperone protein PDE6D via its HVR. PDE6D is also reported to act as a chaperone for RAS isoforms, and most likely interacts with them via their HVR region as well.  The residues comprising the effector, allosteric and hypervariable regions are labeled. The residues in the switch regions are highlighted in yellow. The domain on the interaction partner that is predicted to interact with RAS is listed in parentheses next to the binding partner's name Fig. 4 Description of the pipeline for determining multi-protein RAS complexes. Protein interactions in STRING were manually evaluated and PDB structures were acquired. PRISM was used to predict complex formation. These complexes were combined with available RAS complexes in the PDB to gain a complete view of RAS binding surfaces. This was used to determine which binding partners could bind simultaneously. Simultaneous binding was then assessed in a 3D model to ensure that no steric clashes resulted Interaction partners that do not contain an RBrD [31] and are predicted to bind RAS proteins outside of the switch region include RABGEF1, LGALS1, RARRES3, SRC (HRAS) and CALM (KRAS). We further investigated published data on these RAS interactions to determine their effects on RAS signaling. Interestingly, many of these genes are determinants of RAS protein localization, thereby directing signaling toward or away from canonical RAS effector pathways.
Rab5 GDP/GTP exchange factor (RABGEF1) acts as an E3 ligase for HRAS, and HRAS monoubiquitination promotes its localization to the endosome, thereby decreasing ERK signaling [32]. RABGEF1 possesses a zinc finger (ZnF) domain similar to that of A20 with E3 ligase activity [33][34][35], and Xu et al. [32] suggested that the interaction between RABGEF1 and HRAS occurs via this domain. A mutation in the ZnF domain was reported to obstruct RABGEF1's ability to ubiquitinate RAS, providing support for ZnF-mediated RAS binding.
LGALS1 is a cytosolic protein that has a dual role in RAS biology: it acts both as a chaperone that transports depalmitoylated HRAS to the Golgi and a critical scaffolding protein [36] that is essential for RAS PM association [37]. It has been shown that the RAS inhibitor Farnesylthiosalicylic acid (FTS) simultaneously disrupts RAS PM localization [38] and HRAS interaction with LGALS1 [37]. FTS associates with RAS via the HVR region, and its perturbation of the HRAS-LGALS1 interaction is consistent with our prediction that LGALS1 interacts with the HVR.
The inhibitory interaction between retinoic acid receptor responder protein 3 (RARRES3) and HRAS occurs in the Golgi [39,40]. RARRES3 interaction is reported to suppress EGF-induced RAS activation by decreasing the plasma membrane localization of RAS [40]. The RARRES3 protein lacks well-defined domains, however there are three partitions of its sequence: cytoplasmic, transmembrane and lumenal. Our interface region predictions implicate the cytoplasmic region of RARRES3 as binding to RAS.
SRC is a proto-oncogene that acts upstream of RAS [41], and activates the oncogenic RAS pathway by phosphorylating SHC [42] in tumor cells. HRAS is also reported to directly interact with SRC and negatively regulate its kinase activity; when bound to HRAS, SRC auto-phosphorylation and phosphorylation of its target NR2A are inhibited in the brain [43]. Interestingly, in contrast to other kinases like RAFs, which have RBrDs, SRC binds to HRAS via its kinase domain. Both HRAS and SRC are known for their oncogenic activity, thus the inhibitory activity of HRAS toward SRC suggests tumor suppressive activity. Impaired HRAS signaling, whether due to mutation or therapeutic intervention, could reduce its inhibition of SRC, thereby promoting SRC's oncogenic activity.
CALM is thought to specifically bind KRAS-4B among RAS family proteins [6,8,44,45]. This interaction regulates KRAS membrane association [46], which is crucial for RAS activity in cancer. CALM translocates KRAS from the PM to the cytosol by sequestering its farnesyl group, thereby regulating the concentration of KRAS at the PM [47]. In a range of cell lines, inhibition of CALM induced KRAS and RAF-1 activity [48]. CALM is known to bind KRAS in a RAF/switch region independent manner [8] and the HVR is involved in this interaction [2,6]. Our model suggested that CALM binds HRAS in the allosteric region. Given the absence of the HVR from existing HRAS crystal structures, it is unclear whether CALM actually contacts the allosteric region as well as HVR, or whether building is mediated exclusively by the HVR. Nonetheless, our predictions are consistent with reports that CALM does not bind the switch region.
Based on this evidence and our structural predictions, we propose that when RAS is not localized to the PM, the allosteric and the HVR regions of RAS proteins are involved in interactions with proteins that regulate RAS localization and PTMs. Interestingly, some of these partners are reported to suppress tumorigenic RAS signaling, including RARRES3 [49] and RABGEF1 [50], CALM [51], while others such as LGALS1 [52] and SRC [53] are reported to promote it (Fig. 5). Thus non-switch mediated interactions of RAS may play a significant role in RAS's oncogenic potential.
We next compared PRISM's interface to interface predictions made by similar algorithms including pyDock [54], ZDOCK [55] and COTH [56]. For this comparison, we selected two predictions outside of the Switch regions (HRAS -RARRES3 and KRAS -CALM), and one interaction known to use the Switch region (HRAS -RAF1). While there was some disagreement on the specific residues involved, all methods agreed on whether or not Switch region residues were involved in the interaction (Additional files 1, 2 and 3: Figures S1, S2 and S3). For example, PRISM predicted that RARRES3 does not bind the HRAS switch regions, but instead binds predominantly to residues in the allosteric and hypervariable regions. The other three methods also predicted complexes that did not involve residues in the switch region (with the exception of pyDock which included 3 residues from one Switch motif out of 16 total residues at the predicted binding interface), and implicated multiple residues in the allosteric region (Additional file 1: Figure S1). All methods are capable of returning no result if there is no favorable binding interface. Altogether, this provides strong evidence that RARRES3 interacts with HRAS through an interface that does not involve the canonical switch region binding site.
Finally, we investigated the possibility that proteins binding to RAS isoforms outside of the switch region use a common structural motif. None of the proteins interacted with RAS using the same domain. We next used VAST to assess pairwise structural similarity of these interfaces using the PDB IDs listed in Table 1. This analysis did not identify common structure among interfaces binding outside of the switch region. This may suggest that the disordered HVR alone, or in combination with the allosteric region, allows interaction with multiple distinct structures. VAST did report similarity of binding interfaces on BLC2, NF1 and RALGDS, all of which bind at least partially to the switch region but none of which have an RBrD.

Implications of somatic hotspot mutations for RAS interactions
RAS mutations are most commonly observed at residues 12, 13, 59 and 61 in tumors and these mutants have all been reported to favor RAS's GTP-bound state. For KRAS, 80% of mutations occur at the 12th amino acid, while for NRAS, 60% of mutations affect residue Q61. HRAS, on the other hand, displays a less biased distribution with 50% and 40% of mutations occurring at codons 12 and 61 respectively [57]. In vitro studies have associated G12C mutations with the formation of bulky DNA adducts caused by carcinogens in tobacco smoke [58] and lung cancers are enriched for this specific mutation [57]. Studies of carcinogenesis have also revealed that UV-radiation causes pyrimidine dimers and this process can generate RAS Q61 mutations, which are prevalent in melanoma [59]. Colorectal cancers and haematopoietic/lymphoid cancers harbor codon 13 mutations at high frequency. Thus mutational signatures associated with different exposures and etiologies contribute to tumor-type specific differences in RAS codon mutation prevalence. Different codons are also associated with clinical characteristics of patient tumors. For example, colorectal cancer patients with codon 13 mutations respond to treatment with anti-EGF receptor cetuximabbased therapy, however, patients with mutations at G12 do not [57]. We therefore sought to determine whether specific PPIs could be affected differently by different codon mutations.
We used FoldX [60] to determine the effects of RAS mutations on binding to interaction partners using the new protein interaction models. FoldX has previously been applied to analyze the consequences of cancer mutations affecting the RAS/MAPK pathway [61]. Cancer mutations that were assigned destabilizing energies by FoldX were found to occur in regions of RAS proteins associated with RAS activation, with few predicted to impair protein folding. We adopted a similar approach to examine common effector region mutations (G12, G13, A59 and Q61, 57] and two recently reported mutation hotspots (K117N and A146T) in the allosteric region [62] (Table 2). Our RAS complex models did not predict direct binding to RAS K117 or A146 residues, Fig. 5 HVR and sub-cellular localization of RAS. a Our structural models suggested the involvement of the RAS allosteric region and HVR in protein complex formation. According to the evidence in the literature we hypothesize that such interactions most likely occur when RAS is not anchored to the PM. Our structural models (b) recapitulate the well-established RAS interactors that host an RBrD, (c) and suggest additional interactions that regulate oncogenic RAS signaling however BCL2 was predicted to bind a neighboring surface region on HRAS that spans the residues 147-150.
Several mutants were predicted to stabilize complexes that promote tumorigenesis. A59G and G12D mutations on HRAS and G13D mutations on both HRAS and KRAS are predicted to stabilize interactions with RALGDS. RALGDS is an activator of HRAS and KRAS, thus stabilization of this interaction should increase the amount of RAS protein in the active state, consistent with RAS pathway activation during tumorigenesis. Additionally, we observed that A146T is predicted to stabilize the HRAS -BCL2 interaction that inhibits RAS-induced apoptosis [63].
Interestingly, interactions with FHOD1 were destabilized by some mutants, and stabilized by others. FoldX predicted that G12 and G13 mutations de-stabilize HRAS-FHOD1 binding, while Q61 mutations stabilize it. These conflicting predictions suggest that FHOD1 may promote oncogenic RAS signaling in some settings and interfere with it in others. To further investigate this possibility, we investigated FHOD1 expression in two settings: head and neck cancer and thyroid cancer and thyroid cancer. In head and neck cancer, HRAS mutations are biased toward G12 and G13, whereas in thyroid cancer, HRAS mutations are most commonly observed at Q61. We observed that expression of FHOD1 is generally higher in thyroid tumors than in head and neck cancers, suggesting that Q61 may be more advantageous in thyroid cancer in part because of tumor type specific FHOD1 activity (Fig. 6, Additional file 4: Figure S4). Alternatively, this pattern might suggest that a shift in HRAS binding preference from FHOD1 to another partner such as RALGDS is favorable for head and neck cancer but less favorable for thyroid cancer.

Inferring RAS co-complexes from predicted binding residues
Using the map of residues that mediate interactions with different RAS binding proteins, it is possible to predict  higher order protein complexes that could form. We first eliminated any possible complexes including RAS interaction partner pairs that competed for interface residues on RAS, then used a steric clash filter to determine whether the remaining partners could bind RAS at the same time without interfering with each other (see the Methods; Table 3). This resulted in a list of 38 putative multi-protein complexes composed of RAS first-degree interactors and RAS proteins (Table 3). We found support for several of these predicted complexes in the literature. These RAS multi-protein complexes may help explain how RAS binding partners interact to regulate RAS signaling. Nussinov et al. [2] recently proposed the existence of a multi-protein complex composed of KRAS, PI3K and CALM. According to their model, CALM both sequesters KRAS from the PM and activates PI3Kα. The experimentally validated RAS interactions we used to build our models did not include an interaction between KRAS and PI3Kα, however, one predicted multi-protein complex model included KRAS, CALM and RAF1. Both RAF1 and PI3K bind to RAS using an RBrD. In addition, KRAS has been shown to bind to both RAF and CALM proteins at the same time and the binding of KRAS to CALM is independent from RAF [8].
Galectin-1 (LGALS1) can directly bind HRAS [37] and diverts RAS signals to RAF1 at the expense of PI3K [64]. According to our co-complex models, LGALS1 and switch region-binding partners such as RAF and PIK3CG can bind HRAS at the same time without causing steric clashes. However, it is known that RAS mediated activation of PIK3CG is dependent on recruitment of the PIK3CG catalytic subunit to the PM by the regulatory subunit PIK3R5 [65,66]. Our models suggest that LGALS1 occludes the binding site for the larger PIK3CG-PIK3R5 complex, thereby favoring LGALS1-HRAS-RAF1 complex formation.

Farnesyltransferase and HRAS in complex
Finally, we investigated the RAS CAAX motif modifier Farnesyltransferase (FNT), which is responsible for the attachment of RAS proteins to the PM [67]. The FNT complex, which consists of FNTA and FNTB proteins, binds directly to RAS and modifies the CAAX Cys residue with a 15-carbon farnesyl lipid [67]. Although the structure of the FNT complex is available in the PDB (PDBID: 2H6F), the 3D structure of its interaction with RAS is not available. We therefore modeled the interaction between FNT and HRAS using PRISM [26] (binding energy = −11.0) and observed that FNT most likely binds via the allosteric/HVR regions on the RAS surface, using amino acids 4, 44-48, 157, 158, 161, 164, 165 and 168 (Fig. 7). Our model further predicts that the HRAS-FNT complex is capable of simultaneously binding SOS1, but not ERBB2IP, LGALS1, RABGEF1, RARRES3, or SRC. LGALS1 interaction with HRAS is known to be mutually exclusive to FNT interaction. These proteins are sequential regulators of HRAS localization. LGALS1 recognizes the farnesyl group on HRAS [68] after it has been modified by FNT [67].

Discussion
RAS mutations are found in almost 30% of tumors, making these proteins an attractive target for therapy. However, targeting RAS therapeutically has proven challenging, due in part to incomplete knowledge of RAS protein interactions. Most therapies have focused on targeting effector interactions with the switch regions [18,19,69,70], however new evidence suggests that  Understanding of RAS protein interaction interfaces is still largely restricted to common effectors and proteins with RBrDs, which comprise only~40% of reported RAS binding proteins. Using structural modeling techniques, we were able to investigate protein complex formation between RAS proteins and the experimentally implicated RAS isoform-specific interaction network (Fig. 4). Of these, 45 proteins did not have a RAS binding domain, including post-translational modifiers such as FNTA and FNTB, scaffolding proteins such as PDE6D and LGALS1, ubiquitination pathway proteins such as UBC and BRAP and proteins that are reported to impair oncogenic RAS such as RARRES3 and RABGEF1. Our binding site predictions suggest that multiple RAS binding partners without RBrDs interact with RAS proteins via the allosteric and HVR regions. In particular we found 7 interaction partners (1 KRAS and 6 HRAS) that are predicted to bind the allosteric regions in the same region as NS1.
Interestingly, for 37 reported RAS interaction partners with available structural data, PRISM could not find a viable binding site. These could represent false interactions, or could be false negative predictions resulting from crystal structures that are incomplete or fail to capture a binding-essential conformational state.
The importance of cellular localization as a determinant of RAS signaling is becoming increasingly apparent. It has now been demonstrated that RAS signaling takes place in distinct cytosolic regions including the PM [72], the ER [73] and the Golgi complex [74,75]. Switchregion based interactions with common oncogenic effectors are reported to occur predominantly when RAS is located at the PM [72], whereas many of the interactions that we mapped to allosteric and/or HVR regions have been reported in the cytoplasm [36,40,76]. RAS dimerization is mediated at least in part through the allosteric domain [28,71]. One possible explanation for the apparent correlation between RAS localization and binding region could lie in a propensity for PM localized RAS to dimerize, thereby occluding the allosteric region. Notably, several supposedly cytosolic interactions, such as BCL2 [76], RARRES3 [40] and LGALS1 [36] have clear significance for oncogenic RAS signaling. Schmick et al. [77] speculated that changes in the concentration of RAS at the PM versus at organelle endomembranes could act as a switch for RAS signaling. Interestingly, we find that multiple proteins predicted to bind the RAS allosteric/HVR region are determinants of RAS localization.
In addition to localization, the proportion of RAS bound in particular complexes could have implications for oncogenic signaling and could help explain the activities of particular therapeutic inhibitors. Using our model of RAS isoform interaction interfaces, we were able to predict competitive and synergistic binding of RAS interaction partners. This allowed us to narrow a large combinatorial space of interactions to 38 likely RAS multi-protein complexes that included HRAS-RAF-LGALS1, HRAS-RAF-RARRES3 and KRAS-RAF-CALM. Interestingly, both LGALS1 and Farnesyltransferase use overlapping binding regions on RAS, suggesting mutual exclusivity of these two entities. Given that they have distinct roles in RAS membrane association, this finding may guide studies of the dynamics of RAS membrane association.
Of 6 proteins in the RAS interaction network that are targeted by FDA approved drugs (RAF1, SRC, BRAF, AGTR1, PSMB2, BCL2) [22], two were implicated as binding RAS residues outside of the switch regions, BCL2 and SRC. PRISM did not predict interfaces for AGTR1 and PSMB2, suggesting that these proteins do indeed bind RAS, but they are unlikely to do so via the switch region. Interestingly, the small molecule inhibitor, ABT-263, that targets BCL2, is reported to cause synthetic lethality in KRAS mutant tumors when combined with Trametinib, a MEK inhibitor [78]. Although we did not predict an interaction for BCL2 with KRAS, BCL2 was predicted to bind to both NRAS and HRAS. Our predictions further suggest that BCL2 and RAS effectors involved in MEK signaling, BRAF, ARAF, RAF1 are unlikely to bind to RAS isoforms simultaneously. This knowledge may be helpful for designing follow up experimental studies to further investigate the mechanisms resulting in the synthetic lethality of the ABT-263 Trametinib combination in mutant RAS cells.
Observed patterns of somatic mutations in RAS isoforms across tumor types suggest that biological context determines the advantage of particular mutations. A compelling possible explanation for the specificity of somatic RAS mutations is that tissue level patterns of gene expression results in different ratios of RAS interaction partners, and these differences modify the functionality to RAS mutations. Studying the implications of 5 common RAS mutations for interactions, we found evidence that different RAS binding partners could be affected in different ways by the same mutations. For example, G12 and G13 mutations favored RGALS1 binding, but impaired FHOD1 binding. Interestingly, the HRAS-FHOD1 interaction is subject to both stabilizing (Q61) and destabilizing (G12 and G13) mutations, however these events were associated with distinct tumor types. This suggests that the prevalence of particular RAS mutations in a tumor type may reflect the protein interactions that best contribute to oncogenic RAS signaling in that tumor.

Conclusions
We used structural modeling to map RAS binding interfaces and infer possible multi-protein complexes. Analysis of the resulting binding site map suggests that the allosteric region of RAS may play a more important role in RAS signaling than previously thought, and this role may go beyond RAS dimerization. Protein interactions utilizing the RAS allosteric region appear to be important for RAS localization, which is an important determinant of oncogenic RAS signaling. Further studies are needed to investigate the implications of this new aspect of RAS signaling in cancer.

RAS protein interaction network
We gathered experimentally validated HRAS, KRAS and NRAS protein interactions available in the STRING and PDB databases as of October 2015. We considered PPIs from STRING with confidence score greater than or equal to 0.4 (medium confidence) and manually checked each article associated with the given interaction. We discarded PPIs that had vague implications of physical interaction (such as the experiments that only used HRAS but claimed that the interactions were present for all 3 isoforms) from our dataset. Our network consists of 103 RAS protein interactions with 77 unique proteins; including 68 HRAS, 19 KRAS and 16 NRAS interactors (Fig. 2, Table 4). There is also evidence describing KRAS, NRAS and HRAS dimerization in the literature [28,29]. Edges in the network represent the presence of strong published evidence supporting a physical interaction. We note that there is a strong bias toward probing RAS interactions using HRAS likely due to experimental convenience, and this is reflected in our network.

RAS binding related domains (RBrDs)
Protein domain data was downloaded from the PFAM database in February 2015. We annotated any protein hosting a RAS binding domain (such as RBD, RA, PI3K_RBD, etc.), GAP domain (such as RAS_GAP) or GEF domain (such as RasGEF, RasGEF_N, RhoGEF, etc.) as having a RAS binding related domain (RBrD).  2MSE)). One PDB structure displayed a biological homodimer of HRAS (PDB ID: 3lo5, chains: A and C) but no equivalent structure was found for KRAS or NRAS. We used PISA [79] to differentiate biological homodimers from crystallographic artifacts.

Building structural models of RAS interactions
We aligned PDB structures to discard redundant PDB chains but retained distinct conformations/domains of a given protein. We used TM-align [80] for the structural comparison and considered two PDB chains to be redundant if they received TM-scores larger than 0.5 and had an RMSD under 2.5 Å. We clustered redundant PDB chains together and selected the chain with largest resolution and longest sequence as the representative of the group. On average we retained 2.8 PDB chains per protein (Table 5). We used PRISM [26] to identify physical binding configurations (complexes) between proteins. Previously, PRISM has been successfully applied for modeling PPIs from various pathways including apoptosis [81], ubiquitination [82], MAPK [83], the Toll-like receptor pathways [84], and for identifying possible drug interactions causing off-target effects [85]. In order to assess the performance of rigid-body prediction algorithms, docking benchmarks are widely used. PRISM was capable of building near native models for 87 out of 88 cases in the docking benchmark [27].
Protein complexes predicted by PRISM becomes more energetically favorable as the energy score decreases. The PRISM server reports only the models with energy score smaller than 0, and a threshold of −10 has been used previously in the PRISM literature [85,86]. PRISM v2.0 complexes receiving an energy score of −10 or smaller correspond to the 65th percentile of all models generated by the webserver as of 01/05/2017. For this work, we considered complexes receiving a score of −20 or smaller, corresponding to the 40th percentile (Table  1). At this threshold, we obtained 20 models covering 17 PPIs, with energy scores ranging from −52.35 to −20 (Fig. 3). PRISM has previously reported predictions as extreme as −829. If multiple configurations were implicated for the same RAS partner, the predicted complex with the lowest binding energy was used.
We ran PRISM with multiple distinct KRAS and HRAS chains, but only two chains generated complex predictions: 6q21 (A chain for HRAS) and 4dso (A chain for KRAS). There was only one PDB chain available for NRAS (3con:A) and it was missing the residues between 60 and 72. We therefore generated a complete model for NRAS via Swiss-Model [87] server.
For selected interfaces, PRISM predictions were compared to those of three other tools designed to predict protein complexes pyDock [54], ZDOCK [55] and COTH [56]. Predictions by these tools were obtained from their respective web servers by entering PDB IDs, and the top scoring predicted complex was selected.

Interface region extraction
We used the consensus set of KFC2 [88] and HotPoint [89] server predictions to extract the specific amino acid positions of the interface region from PRISM predicted complexes. The structures of the domains containing interface regions for proteins predicted to bind outside of the Switch domain were compared using VAST [90] via the NCBI structure database (https://www.ncbi.nlm.nih.gov/structure).

Steric clash filter
To determine whether RAS binding partners could bind RAS simultaneously, we used structural data to investigate whether simultaneous binding would create steric clashes. To do this we aligned 11 HRAS/4 KRAS complexes (both experimentally derived (Additional file 5: Figure S5) and predicted RAS complexes) to a reference HRAS/KRAS protein (HRAS PDB ID: 6Q21:A; KRAS PDB ID: 4DSO:A) via the "Profit" [91] protein aligner. Then, we calculated the Euclidean distance between the atoms of each protein in the complex. If the distance between any atom pair was less than 1 Å, we assumed the RAS interaction partners would result in steric clashes and therefore could not simultaneously bind RAS. We note that we used the HRAS-RAF1 complex available in PDB instead of the predicted complex generated by PRISM to evaluate steric clashes between RAF1 and other HRAS binding partners.
Simulating the effects of mutations on protein-protein interactions RAS hotspot mutations and tissue specificities were obtained from the literature [57,62]. We used FoldX [60] to simulate the effects of mutations on PRISM-predicted protein complexes. For this purpose we calculated the binding affinity change after a mutation was introduced. As suggested in the FoldX manual, we assumed that free energy changes with a magnitude smaller than 0.5 were neutral. Larger positive free energy changes correspond to a destabilizing effect on complex formation, whereas larger negative free energy changes correspond to a stabilizing effect on complex formation.