The Structural Pathway of Interleukin 1 (IL-1) Initiated Signaling Reveals Mechanisms of Oncogenic Mutations and SNPs in Inflammation and Cancer

Interleukin-1 (IL-1) is a large cytokine family closely related to innate immunity and inflammation. IL-1 proteins are key players in signaling pathways such as apoptosis, TLR, MAPK, NLR and NF-κB. The IL-1 pathway is also associated with cancer, and chronic inflammation increases the risk of tumor development via oncogenic mutations. Here we illustrate that the structures of interfaces between proteins in this pathway bearing the mutations may reveal how. Proteins are frequently regulated via their interactions, which can turn them ON or OFF. We show that oncogenic mutations are significantly at or adjoining interface regions, and can abolish (or enhance) the protein-protein interaction, making the protein constitutively active (or inactive, if it is a repressor). We combine known structures of protein-protein complexes and those that we have predicted for the IL-1 pathway, and integrate them with literature information. In the reconstructed pathway there are 104 interactions between proteins whose three dimensional structures are experimentally identified; only 15 have experimentally-determined structures of the interacting complexes. By predicting the protein-protein complexes throughout the pathway via the PRISM algorithm, the structural coverage increases from 15% to 71%. In silico mutagenesis and comparison of the predicted binding energies reveal the mechanisms of how oncogenic and single nucleotide polymorphism (SNP) mutations can abrogate the interactions or increase the binding affinity of the mutant to the native partner. Computational mapping of mutations on the interface of the predicted complexes may constitute a powerful strategy to explain the mechanisms of activation/inhibition. It can also help explain how an oncogenic mutation or SNP works.


Introduction
Interleukin-1 (IL-1) is a large family of cytokines (small cell signaling proteins) that mediate innate immune responses to defend the host against pathogens. The IL-1 family has 11 member proteins (IL-1F1 to IL-1F11) and they are encoded by 11 distinct human and mouse genes [1][2][3]. The first discovered family members, IL-1a (newly named IL-1F1) and IL-1b (IL-1F2), are secreted by macrophages and epithelial cells in response to pathogens and have strong proinflammatory properties leading to fever (affecting hypothalamus) and activation of T cells and macrophages.
IL-1 family members have been intensely studied (especially IL-1a and IL-1b) unraveling their roles in a number of autoinflammatory diseases [4][5][6]. Signaling initiated by the IL-1 cytokines increases the expression of adhesion factors on endothelial cells resulting in immune cells (such as phagocytes and lymphocytes) migration to the site of infection. The autoinflammatory disease is a class of chronic inflammation with increased secretion of active IL-1b, thus blocking IL-1b is therapeutically beneficial [7].
IL-1a and IL-1b can induce mRNA expression of hundreds of genes, including themselves (a positive-feedback loop), and their gene regulatory actions are conducted via a conserved signaling pathway [8]. Signal propagation mainly depends on mitogenactivated protein kinases (MAPKs), MAPK kinases (MKK/ MAP2Ks), MKK kinases (MKKK/MAP3K/MEKKs) and the downstream proteins of MAPKs, finally leading to activation of transcription factors that regulate the expression of host defense proteins ( Figure 1). The signal initiates by binding of IL-1a or IL-1b ligands to type I receptor (IL-1R1) and propagates with the help of the co-receptor IL-1 receptor accessory protein (IL-1RAP), forming a trimeric complex [9]. In this trimeric complex, the Tolland IL-1R-like (TIR) domains on the cytoplasmic regions of IL-1R1 and IL-1RAP receptors get close to each other resulting in the recruitment of myeloid differentiation primary response gene 88 (MYD88), Toll-interacting protein (TOLLIP) [7] and IL-1 receptor-associated kinase 4 (IRAK4) [10,11]. A stable complex is formed between IL-1, IL-1R1, IL-1RAP, MYD88 and IRAK4 [10]. MYD88 binding triggers phosphorylation of IL-1 receptorassociated kinases IRAK4, IRAK2 and IRAK1, leading to the recruitment and oligomerization of tumor necrosis factor-associated factor 6 (TRAF6) [12][13][14]. TRAF6 and phosphorylated IRAK1 and IRAK2 dissociate and migrate to the membrane to associate with TGF-b-activated kinase 1 (TAK1) and TAK1binding proteins TAB1 and TAB2 [7]. The TAK1-TAB1-TAB2-TRAF6 complex migrates back to the cytosol, where TRAF6 is ubiquitinated and TAK1 is phosphorylated [7]. From this point, the signal can propagate via two main paths: IKK -IkB -NF-kB and/or MKK -MAPK/JNK/ERK ( Figure 2). In the first path, phosphorylated TAK1 activates the inhibitor of nuclear factor kappa-B kinase subunit beta (IKKb) and activated IKKb phosphorylates the nuclear factor kappa-B inhibitor (IkB) which gets degraded so that nuclear factor kappa-B kinase (NF-kB) is released and migrates to the nucleus [7]. TAK1 can also activate mitogenactivated kinases (MAPK) p38, c-Jun N-terminal kinases (JNK) and extracellular signal-regulated kinases (ERK) by interacting with MAP kinase kinase (MKK) proteins. Downstream in this path, are transcription factors such as c-Jun, c-Fos, c-Myc and ATF2.
MAP kinase signaling pathways, which are conserved among eukaryotes, mediate cellular events triggered by extracellular signals such as cytokine binding [15] and they are essential for IL-1 signaling (Figures 1 and 2). This pathway builds upon a triple kinase cascade consisting of a MAP kinase kinase kinase (MKKK/ MEKK), a MAP kinase kinase (MKK/MEK) and a MAP kinase (MAPK) and these kinases sequentially phosphorylate and activate each other [15]. The JNK and p38 MAP kinases, called stress activated MAP kinases, have roles in tumor suppression and can be both directly phosphorylated and activated by MKK4, which is also a tumor suppressor [16][17][18]. The successive activation mechanism takes place as follows: MEKK interacts with inactive MKK and phosphorylates it; the complex dissociates, releasing the free and active MKK, which is ready to interact with inactive JNK to activate it [15]. Activation of JNK leads to disruption of the MKK-JNK interaction, freeing the active JNK to phosphorylate its downstream targets. There are several mechanisms through which stress activated MAP kinases regulate tumor suppression, including promoting apoptosis (p53, Bax, Bim/Bmf), inhibiting proteins that inhibit apoptosis (Bcl2, Bcl-XL, 14-3-3, Mcl-1), inhibiting tumor development (TGF-b1) and tumor growth (CDC25, CyclinD1/CDK4) [16]. Somatic mutations were identified in the JNK pathway via large-scale sequencing analyses of human tumor cells [19,20], and JNK3 encoding gene (MAPK10) has been speculated to be a putative tumor suppressor gene as almost half of the brain tumors that were examined included mutations [21]. ERK1 and ERK2, the other members of MAPK family, are also upregulated in tumors [22].
Recently, inflammation has been related to cancer [23][24][25]. Cancers are mostly due to somatic mutations and environmental factors, and chronic inflammation is implicated by most of these risk factors [26]. Chronic inflammation, due to autoimmune diseases or infections, causes tumor development via several mechanisms, including oncogenic mutation induction [26]. Oncogenic mutations and single nucleotide polymorphisms (SNPs) are key players in inflammation-related cancers and it is crucial to map the mutations/SNPs on the corresponding 3D structures of the proteins to gain insight into how they affect protein function [27,28]. SNPs that cause diseases, if not in the core of the protein, are frequently located in protein-protein interface regions rather than elsewhere on the surface [26,28]. Structural knowledge can clarify the conformational and functional impact of the mutation/ SNP on the protein [27][28][29]. The effect of a functional mutation can be expressed by a change in the specificity of the interactions between a mutated protein and its partners [30]. Quantitatively, a mutation changes the binding free energies of the mutant's interactions with its partners with respect to the free energies of its interactions in the native form [30]. The functional impact of the mutation on the specificity differs. The mutation can destabilize the protein and/or its interaction, leading to 'loss-of-function'; or can lead to a change in the specificity of protein-partner interactions, resulting in a 'gain-of-function', or can gain new binding partners and hence a new biological function, i.e. result in a 'switch-of-function' [30].
Two recent studies used structural pathways for mapping mutations on protein-protein interfaces, one on smaller scale

Author Summary
Structural pathways are important because they provide insight into signaling mechanisms; help understand the mechanism of disease-related mutations; and help in drug discovery. While extremely useful, common pathway diagrams lacking structural data are unable to provide mechanistic insight to explain oncogenic mutations or SNPs. Here we focus on the construction of the IL-1 structural pathway and map oncogenic mutations and SNPs to complexes in this pathway. Our results indicate that computational modeling of protein-protein interactions on a large scale can provide accurate, structural atom-level detail of signaling pathways in the human cell and help delineate the mechanism through which a mutation leads to disease. We show that the mutations either thwart the interactions, activating the proteins even in their absence or stabilize them, leading to the same uncontrolled outcome. Computational mapping of mutations on the interface of the predicted complexes may constitute an effective strategy to explain the mechanisms of mutations-constitutive activation or deactivation.
pathways [27] and the other on large scale [28]. Mosca et al. [27] mapped mutations onto proteins as an application of their useful computational modeling technique with the data limited to the Kyoto Encyclopedia of Genes and Genomes (KEGG) complement cascade pathway and the interactions of complement component 3 (C3) and complement factor H (CFH) which it includes. In a pioneering work, Wang et al. [28] explored genotype-phenotype relationships on a large scale. They systematically examined thousands of mutations and mapped them onto interaction interfaces and experimentally validating their predictions for the MLH1-PMS2, WASP-CDC42 and TP63-TP73 interactions. These and other studies emphasized the need for computational methods for large-scale interactome studies [31,32]. We propose a method similar to ones used in the works of Mosca et al. [27] and Wang et al. [28], while introducing the advantage of in silico mutagenesis to observe the effects of mutations on protein-protein interactions on a large scale.
Here, we construct the IL-1 signaling pathway by combining the related pathways and information from the literature [8,[33][34][35] (Figure 2). We observe that there are approximately 100 interactions between proteins that have experimentally identified 3D structures. However, only 15 of the interactions have structures of protein-protein complexes ( Figures 3A and S1) in the Protein Data Bank (PDB, http://www.rcsb.org/pdb/) [36]. The structural coverage of the pathway is under 20z%. Our major aim is to expand the structural apoptosis pathway [37], and map oncogenic mutations and SNPs to reveal their mechanism.

Results/Discussion
Reconstruction and structural analysis of IL-1 pathway We are able to increase the structural coverage from 15% up to 71% by predicting the protein-protein complexes throughout the pathway using the PRISM algorithm [38][39][40], such that the  The Structural IL-1 Signaling Pathway and Cancer PLOS Computational Biology | www.ploscompbiol.org number of predicted structures of protein-protein complexes that have binding energies lower than 210 energy units is 74 out of 104 interactions ( Figure 3). The predicted structures of complexes are also compared with the 15 available PDB structures of the complexes in the pathway. Except for three where the proteins interact with peptides, all PDB structures could be successfully reproduced via the predictions (Table S1).

Mapping and analysis of SNPs and oncogenic mutations in the IL-1 pathway
A large scale analysis of the distribution of oncogenic mutations and SNPs in the predicted and experimental structures of proteinprotein complexes in the IL-1 pathway indicated that both oncogenic mutations and SNPs correspond significantly to interface regions and their 'nearby' residues (Tables 1-4 and  Table S2) with p-values 0.00013 and 0.0009, respectively. Nearby residues are defined as the residues that are at most 6 Å away from interface residues in three dimensional space. Along similar lines, David et al. showed that SNPs, if not in the protein core, are more likely to occur at the protein-protein interfaces rather than the remaining surface [41].

Oncogenic mutations in the interfaces of MKK4 with JNK2 and JNK3
Three oncogenic mutations and a SNP are observed either directly on the predicted interfaces of MKK4 with JNK2 and JNK3 (the SNP is a computational hot spot residue in the MKK4-JNK3 interaction) or as nearby residues ( Figure 4). In Figure 4A, the predicted MKK4-JNK2 interaction has a binding energy of 212.29 energy units. Analysis of interface residues of the predicted complex reveals that MKK4 Ser251 and JNK2 Gly268 correspond to the interface region. These residues are important because the p.Ser251Asn mutation is involved in metastatic melanoma and p.Gly268Ala is a SNP. Figure 4B shows the predicted complex of MKK4-JNK3. On MKK4, Arg154 is mapped to the interface as a computational hot spot residue and its mutation to Tryptophan (p.Arg154Trp) is involved in colorectal adenocarcinoma. Moreover, Gln142 is a nearby residue of the interface and the MKK4 p.Gln142Leu mutation is involved in lung squamous cell carcinoma. Although it is known that MKK4 interacts with MAPKs (JNK1/2/3 and p38a/b) through its MAPK-docking site (D-site) [42], here we identify a complementary binding interface between MKK4 and JNK2/3. The D-site is located at the N-termini of MKK4 (residues 37-52), but the Ntermini is missing in all of the 3D structures (PDB codes: 3aln, 3alo and 3vut; including residues between 80 and 399) and hence could not be included in our predictions. However, we modeled a complementary interface, which is not at the D-site but in the kinase catalytic domain. This model can be explained by the finding that although the D-site facilitates the activation of MAPKs through MKK-MAPK signaling, the specificity is also affected by allosteric cooperativity among other binding motifs in the kinase catalytic domains [43]. The mutations involve important residues in the critical domains of MKK4 such as Ser251 in the activation loop and Gln142 in the aC-helix. Gln142 was also proposed by another study [30] to be in the interface of the MKK4-human kinase STK4 interaction and the p.Gln142Leu mutation is speculated to switch-off native partners, gathering nonnative interaction partners instead. Also, the MKK4 protein has a large interface for the dimer formed in the experimentally identified crystal structure (PDB code: 3alnAB) which includes Gln142. Similar to the dimer in the PDB, in our predictions both MKK4 and JNK3 interact with each other through their kinase domains, forming a large interface area including Gln142 as a nearby residue. Similarly, a putative interface of MKK-MAPK interaction was found to overlap with the ERK2 dimerization interface [44]. Wilsbacher et al. also listed the important regions on the MAPK for MKK-MAPK binding; one of them is the tip of the C helix [44], which is in accordance with our model as Gln142 in the aC-helix of MKK4 is one of the mutants and found to be important in MKK4-JNK3 interaction. In short, the predicted interaction does not show a simple recognition (D-site) binding but rather involves a specificity binding model.

The mechanisms of the MKK4 mutations
In order to better understand the mechanism through which the mutations may relate to cancer, these residues are mutated computationally and the resulting predicted binding energies of the complexes are compared (Table 5). After individual energy minimization of the target structures and re-running PRISM on the energetically minimized wild type and mutant targets, we observe models of protein-protein complexes with slightly different interface residues and conformations for the minimum energy solutions. The comparisons are based on these reference interactions. For the MKK4-JNK3 interaction, results are also obtained using a different template (1p4oAB) than the original (2ab0AB). None of the complexes predicted using the original template is favorable due to the predicted positive binding energy values after the minimization and we are unable to assess mutations based on unfavorable reference interaction. However, predictions of the MKK4-JNK3 interaction based on the 1p4oAB template clearly show the effect of p.Gln142Leu and p.Arg154Trp mutations as the wild type complex is still favorable after the minimization (212.66 energy units, Table 5 and Figure 5A). The mutation of Arg154 to Tryptophan abolishes the interaction as the predicted binding energy becomes positive (12.84 energy units), implying that the interaction is not favorable anymore ( Figure 5B). We can explain this change in terms of the specificity as an effect of a loss-of-function type mutation, as MKK4 cannot phosphorylate and activate JNK3 and its tumor-suppression function is lost. On the other hand, Gln142 is a nearby residue and its mutation to Leucine causes a dramatic decrease in the predicted binding energy (241.14 energy units) with respect to the newly predicted reference MKK4-JNK3 interaction (212.66 energy units) ( Table 5 and Figure 5C). We observed that the interaction not only still takes place after the mutation occurs but also it is predicted to have a tighter binding.
The contribution of this change to cancer development can be interpreted in two ways: a stronger complex between MKK4 and JNK3 takes place and the interaction, previously transient, gets so strong that it becomes obligate and cannot dissociate so that the activated JNK3 cannot activate its downstream targets, resulting in tumor in predisposed persons; or the essentially constitutively active MKK4 phosphorylates the JNK3 targets and the elevated, unregulated JNK activation leads to tumors as shown previously [45]. A similar case was also previously shown for EGFR, the activation loop of which could not switch to the inactive conformation due to a mutation. Hashimoto et al. showed that single cancer mutations in kinase domains destabilized inactive states [46]. Moreover, our structural analysis explains the experimental data in the study of Ahn et al. [18], who performed site-directed mutagenesis using human MAP2K4 (MKK4) cDNA in order to create somatic mutants which are subjected to a mutant JNK1 as a substrate and concluded that p.Arg154Trp is one of the loss-of-function mutations whereas p.Gln142Leu is a gain-offunction mutation which results in a highly active MKK4 similar to a synthetic constitutively active mutant. In the second case study, the analysis of the predicted MKK7-JNK3 complex reveals the effect of the p.Arg178Cys and p.Arg178His mutations in MKK7. Similarly, the individual targets are energetically minimized, in silico mutations are performed, and PRISM is re-run to obtain new values of predicted binding energies (Table 5). When the values are compared to that of the new reference MKK7-JNK3 interaction (211.66 energy units), a jump in both mutant cases (Table 5) and a significant change in the conformations are observed, making the energies positive so that these interactions cannot take place anymore. Thus, the mutation of Arg178 to Cys and His are both inferred as loss-offunction mutations abrogating the interaction with the JNK3 partner. Since MKK7 is a tumor suppressor, the inhibition mechanism of the MKK7-JNK3 interaction could presumably lead to tumor progression due to the loss-of-function of MKK7.

The interactions of IL-1 receptors IL1R1 and IL1RAP
Finally, we concentrate on two IL-1 pathway specific interactions, namely the interactions of IL-1 with receptors IL1R1 and IL1RAP. After mapping the oncogenic mutations and SNPs onto the interfaces of these receptors, we observe that the predicted structure for IL-1a-IL1R1 interaction is important due to containing 11 SNPs and 5 mutations on its interface and nearby residues (see the last row of Table 2). Comparison of the predicted binding interface to experimentally reported interactions reveals common residues, confirming the predicted complex structure of IL-1a-IL1R1. In the work of Labriolatompkins et al. [47], using oligonucleotide-directed mutagenesis, they determined seven important residues on the binding site of human IL-1a for IL1R1 binding. When we mapped these residues on the PDB structures and compared them to the interface of the predicted wild type IL-1a-IL1R1 interaction (with 271.92 units of predicted binding energy), we observed that 4/7 residues (Arg16, Ile18, Ile68, Trp113) are in common and the remaining three residues (Asp64, Asp65, Lys100) are nearby residues of the interface.
After confirming that the predicted structure is in accordance with related data in the literature, the binding energies of the mutants are compared to the energetically minimized reference interaction (219.5 energy units). Being a ligand-receptor interaction, this complex has a large interface containing 75 residues in addition to 104 nearby residues. Intuitively, we expect to observe insignificant effects of SNPs or single point mutations on the interaction. The results are in parallel with our expectation and there are no significant changes in the predicted binding energies for the interactions with mutant proteins. However, we notice an exception for the p.Ile68Asn SNP on IL-1a ligand (Table 5). This SNP is observed to be very important for binding as it blocks the interaction of the mutant ligand not only with the wild type receptor but also with most of the mutant receptors (Table 5). This observation may be supported by the information that Ile68 is both a computational hot spot residue and experimentally determined as critical in binding [47]. Thus, importantly, this SNP may affect the innate immune system and inflammatory response as this interaction is critical for the initiation of IL-1 signaling, with the only alternative pathway being the IL1B(IL-1b)-IL1R1 interaction. However, interestingly this is not always the case. A ligand bearing the p.Ile68Asn substitution (a computational hot spot) appears to interact with a p.Ile276Thr mutant receptor (a nearby residue of the interface) with the energy of 240.64, suggesting that this mutation on the receptor is in fact useful for canceling out the effect of the SNP on the ligand ( Figure 6). Although these two residues are not in contact in three dimensional space (with a distance of 16.1 Å between them), the change in the binding energy is significant implying that the difference stems from slight conformational changes. This observation reveals that compensatory changes do not necessarily involve residues in contact. The p.Ile276Thr mutation on IL1R1 is important: according to the COSMIC database [48,49] it is related to endometrioid carcinoma and IL1R1 expression was found to be increased in active endometriotic lesions [50]. This fact supports our finding that IL-1a-IL1R1 interaction is still favorable between this mutant receptor and IL-1a with the p.Ile68Asn SNP, and emphasizes the usefulness of our strategy.
To conclude, in the absence of experimental data identifying the location of the mutations and SNPs on the protein structures, we are able to computationally clarify the mechanism of inhibition/activation by mapping these mutations on the interface of the predicted MKK4-JNK3, MKK7-JNK3 and IL-1a-IL1R1 complexes; and our account is corroborated by available experimental data.

Conclusions
Here, we reconstructed the IL-1 signaling pathway by combining related pathways and information in the literature, expanding the previously constructed structural apoptosis pathway. By predicting protein-protein complexes throughout the pathway using the PRISM algorithm, the structural coverage of the pathway increased up to 71%. The distributions of oncogenic mutations and SNPs in the predicted structures of protein-protein complexes indicated that they significantly correspond to interface and adjoining residues, and more importantly, in some cases to computational hot spot residues. While oncogenic mutation and SNP data are reported for single proteins, by mapping them onto interfaces we are able to determine the critical binding partners, interactions with whom are affected by the mutations. Additionally, in silico mutagenesis of the corresponding residues and comparison of the change in the binding energies between the wild type and mutant shed light on the mechanism of cancer

Materials and Methods
The reconstructed IL-1 signaling pathway (Figure 2) is used as the target pathway for predicting protein-protein interactions.

Protein-protein interaction prediction algorithm -PRISM
The 3D structures of protein-protein complexes in the reconstructed IL-1 signaling pathway are predicted by using the PRISM (PRotein Interactions by Structural Matching) [38][39][40] tool (Figure 7). PRISM is a large-scale protein-protein interaction prediction, modeling and structure assembly tool [40] that previously was successfully applied to several pathways [37,[52][53][54]. Two input sets, the template and the target sets, are required by the PRISM algorithm to obtain protein-protein interaction predictions (Figure 7). The structurally nonredundant (nonhomologous) unique interface dataset described by Tuncbag et al. [55] is used as the template set in this work. The dataset was generated by hierarchical clustering of 49512 two-chain interfaces (extracted from all of the available PDB complexes in the version of February 2006) into 8205 clusters. In this work, one representative interface for each of the clusters is combined into the template set. However, the number decreases from 8205 to 7922 due to the update in the PDB (since 2006) and replacement of previous structures with the better ones under different PDB codes. PRISM is a prediction algorithm based on threedimensional protein structures and therefore it can only be applied to target proteins with known (experimental or high quality modeled) structures. The target set may contain a minimum of two proteins and the number of proteins in the set can increase up to any desired number, that is, all the proteins in a given pathway. Here, the target set is composed of IL-1 related PDB chains, the interactions and complex structures of which we want to predict. We focused on every single one of these and applied the prediction algorithm for each of the binary interactions (104 edges) between the targets that have 3D structures (50 nodes) (Table S3). However, for one interaction there might be more than two structures in the target set as each of the proteins may have more than one experimentally identified structure. We included all available PDB codes for the target proteins in each run to have a complete picture of possible protein-protein complexes.
The prediction algorithm has four main steps: target surface extraction, structural alignment of the template interface and the targets, transformation of the targets onto the template and eliminating collisions, and flexible refinement of the resulting complexes ( Figure 7). Firstly, the surface residues of target chains are extracted via the Naccess program [56] which calculates the accessible surface area of residues. The criterion for a residue to be accepted as a surface residue is that its relative surface accessibility should be greater than 5%. In the second step, the structural comparison of the template interface chains with the target chain surfaces are carried out using the MultiProt structural alignment tool [57] based on some filters. For example, if the template chain has less than or equal to 50 residues, then 50% of the template residues should match the target surface residues; if larger than 50, a 30% match of template to target residues is required. In addition, at least one 'hot spot' residue on the template interface should match one of the hot spots on the target surface. In the third step, the resulting set of target surfaces from the previous set are transformed onto the corresponding template interfaces to form a complex so that these two targets are potential partners interacting with each other through an interface similar to the template interface architecture. Following transformation, PRISM eliminates the residues of the target chains that collide with each other. Finally, the FiberDock algorithm [58,59] is used to refine the interactions allowing some flexibility, to resolve steric clashes of side chains, to compute the global energy of the complex and to rank the solutions based on the calculated binding energies. The calculated binding score is correlated with the experimental binding free energy and taken as the approximation of the binding free energy function [58,60]. We accept a threshold of 210 energy units for a favorable interaction [54]. The detection of computational hot spot residues that are used in the analyses is done by using HotRegion Server [61].

Computational analysis of the SNPs and mutations on the structures of the protein-protein complexes
Once all possible structures of protein-protein complexes in the IL-1 signaling pathway are predicted, the analysis of the interfaces  in terms of SNP and mutation distribution is done. In addition to mapping the available SNPs and mutations onto predicted complexes and checking whether they correspond to the interface, in silico mutagenesis is also performed to observe the change in the predicted binding energy of the wild type and mutant complexes. The list of SNPs and mutations related to the target proteins is obtained from LS-SNP/PDB which is a web tool for genome-wide annotations of human non-synonymous SNPs mapped to Protein Data Bank structures [62] and COSMIC (Catalogue of Somatic Mutations in Cancer) database [48,49], respectively. Based on the list, 394 SNPs and 535 oncogenic mutations are mapped onto the target PDB structures of the proteins in IL-1 pathway and these numbers decrease to 371 and 302, respectively when only the SNPs and mutations on the PDB surfaces are considered (Table 1). Then, these SNPs and mutations are mapped onto the interface residues and their nearby residues of both predicted and  experimental interactions in the IL-1 pathway ( Table 2). The statistical significance of these mappings are calculated by applying t-test to the following two sets: the ratio of mutations and SNPs mapped on interface and nearby residues to total interface size including nearby residues versus the ratio of mutations and SNPs on a noninterface surface region to total noninterface surface region (Tables 3 and 4). These ratios are calculated for each of the 77 interactions listed in Table 2. Note that the surface residues of target chains are extracted via Naccess program [56] as described in the previous subsection.
To understand the effect of oncogenic mutations and SNPs, those predicted to be at the interface (or nearby residues) are mutated in case studies using the FoldX plugin [63] for the YASARA molecular viewer [64]. There are two steps for computationally mutating the residues: energy minimization and residue change. The ''RepairPDB'' module of YASARA is used for energy minimization both before and after the residue change and the ''BuildModel'' module is used for mutation. For setting the FoldX options, default values are used except temperature, which is taken the same as the experimental temperature of the PDB The wild type target structures are energetically minimized, related residues are mutated and PRISM is re-run to predict the protein-protein interactions. *The new reference interactions for comparisons are the predicted interactions between energetically minimized wild type targets and for the three cases they are: MKK4-JNK3 interaction predicted using 1p4oAB template with a predicted binding energy of 212.66 energy units; MKK7-JNK3 predicted using 1ftaAC template with a predicted binding energy of 211.66 energy units; and IL1A-IL1R1 predicted using 1itbAB template with a predicted binding energy of 219.5 energy units.   structure. After obtaining the mutant structures, an energy minimization process is performed for the second time and the resulting structures are saved as PDB files and used as new targets for re-running the interaction prediction algorithm PRISM. The predicted binding energy values of the interactions are then compared to the reference values (Table 5). Figure S1 Experimentally determined structures of protein-protein complexes mapped to the IL-1 signaling pathway. The PDB codes used in the figure are: IL1b-IL1R1-IL1RAP, 4dep; MYD88-IRAK2-IRAK4, 3mop; IKKa/IKKb-IKKc, 3brt; MKK6-MAPKp38a, 2y8o; MAPKp38a-MAPKAPK2, 2onl; JNK3-ATF2, 4h36; ERK2-Mnk1, 2y9q; c-Jun-c-Fos, 1fos; NF-kBp105-NF-kBp65, 3gut; NF-kBp50-NF-kBp65-IkBa, 1nfi. The blue color represents the proteins that precede its partners in the information flow.