Introduction

Genomes of many organisms encode large number of proteins with multiple domains1. Multi-domain proteins have an evolutionary advantage over large single domain proteins with respect to folding2. Further, multiple domains provide proteins with both structural3,4,5 (domain motion/interaction) and functional plasticity (involving new functional sites/binding/regulatory/allosteric sites). The evolution of various domain combinations in proteins is primarily governed by major recombination events (Duplication/ insertion/ deletion/ transposition)6,7,8,9. Further the evolutionary selection pressure for the newly generated domain combination is governed by functional advantage provided to the organism10. An ideal fitness function which could model the evolution of domains requires a biophysical term and a functional term10,11,12. The biophysical term serves to channel the sequence to remain in the stability zone (ΔGfolding) for the proteins to function normally under physiological conditions12. The functional term is a combination of the overall efficiency of the reactions co-ordinated by the given sequence13,14. The domain stabilities and the strength of inter-domain interactions in a multi-domain protein have an enormous influence on the functioning of a multi-domain protein15,16,17. The inter-domain communications facilitate proper functioning of most multi-domain proteins18,19. The most common way to understand the emergent effects of domain-domain interactions on stability and folding employs the comparative analysis of energetics of the constituent domains in isolation along with the native full length form of the protein. A very common situation of obtaining single domains of multi-domain proteins with stable (native) structure for X-ray or NMR analysis is often non trivial. No clear guidelines are currently available for designing constructs coding for single domain of a multi-domain protein that is supported by a thorough understanding of origin of destability of domain structures of multi-domain proteins when considered in isolation. Till date the intricate relationship between the domain stability, inter-domain interactions and independent stable existence remains obscure. In our current study, using biophysical analysis of the domain stability and their interactions in multi-domain proteins, we aim to understand the basis of inability of certain domains to exist independently with native fold. The reasons for not observing certain sequences to exist independently in the extant proteins are manifold. 1. Innate nature of the sequence which cannot fold into a stable unit/conformation. 2. Such sequences have not been explored by evolution and hence are not observed in extant proteins. 3. Such sequences can exist stably where these proteins have evolved to interact with other proteins in an obligate fashion (Homo/hetero oligomers). 4. Incompleteness of the data: Current sequenced genomes lack the domain of interest as a single domain protein. We explore the possibility of lack of foldability and/or stability of domains of multi-domain proteins to exist independently. We further quantify the contributions of domain-domain interactions and solvent-protein interactions in order to understand how far these factors affect the stability in isolation. We also address how lack of inter-domain interactions in stable single domain homologues has contributed to differential selection pressure of non-interacting surface residues and subsequent evolution of these domain folds. Information from stable single-domain homologues of multi-domain proteins has been tailored with energetics calculations in order to arrive at guidelines for safe residue substitutions by site-directed mutagenesis to artificially stabilize a single domain of a multi-domain protein with native fold, but, in isolation.

Results

Stability of domains in single and multi-domain proteins

We grouped all the protein structures into four main datasets, viz; MD_doms, MD_chains, 1D_homs and 1D_uniq (See Table 1) corresponding to domains of multi-domain proteins, full-length gene products of multi-domain proteins, single domain homologues of domains in multi-domain proteins and a control, independent data set of single domain proteins respectively. We computed the stabilities (ΔGFolding) of each of the proteins/domains in these data sets using the empirical effective energy function, FoldX20 (See Methods). We found that domains of multi-domain proteins considered in isolation (MD_doms; 7.70 ± 0.96 Kcal/mol) are significantly less stable in comparison to (1) full length polypeptide chains (MD_chains; −6.57 ± 2.37 Kcal/mol) (Unpaired t-test: t = 6.603; df = 1282; P = 5.89E-11), (2) single domain homologues of domains in multi-domain proteins (1D_homs; −2.15 ± 2.34 Kcal/mol) (Unpaired t-test: two-tailed; t = 4.534; df = 1115; P = 6.41E-06) and (3) the control dataset of single-domain proteins (1D_uniq; −0. 78 ± 2.77 Kcal/mol) (Unpaired t-test: two-tailed; t = 3.237; df = 996; P = 1.25E-03) (Figure 1). On the contrary the stabilities of full length multi-domain proteins (MD_chains) is not very different from that of independently stable single domain homologues (1D_homs) (Unpaired t-test: t = 1.244; df = 687; P = 0.2139) and control set of single domain proteins (1D_uniq; −0. 78 ± 2.77 Kcal/mol) (Unpaired t-test: t = 1.309; df = 568; P = 0.1910) highlighting the probable importance of inter-domain interaction energy involved in stabilizing otherwise unstable domains (MD_doms).

Table 1 The Datasets used in the current study.
Full size table

Inter-domain interactions and effects on stability

We observed that domain-domain interactions (See methods) across the interface primarily comprised of short range (van der Waal's) contacts between the hydrophobic residues (Figure 2). Most of the contacts are made by Leu, Ile, Phe, Ala and Tyr (dark blue contours in bottom left zone in Figure 2). Few interfaces are observed to be studded with salt bridges (especially Arg-Asp/Glu and Lys-Asp/Glu) (Figure 2; top right quadrant). The strengths of these interactions are proportional to the observed interaction interface surface area (Linear regression: r2 = 0.9273; P < 0.0001; ΔGInt = −0.0157*SAInterface + 5.092). in silico Alanine scanning revealed that the residues (Primarily Hydrophobic) involved in most frequent interactions across the interface are also the most energetically favoured interactions in the domain-domain interfaces (Figure 3), thereby acting as hotspots in the domain-domain interface. Figure 3 also shows that change in interaction energies due to Ala mutations is correlated with the volume/size of the interface residues in the wild type.

Figure 2
figure 2

Interaction preferences in the domain-domain interface of multi-domain proteins.

Shown is the contour plot of frequencies of pair wise residue contacts between the two-interacting domains of a set of two-domain proteins (n = 428) . A distance based approach was used to determine the van der Waal's contacts between the two interacting domains (See methods). Note the preponderance of non-polar contacts and salt bridges (Bottom left quadrant and top right quadrant).

Full size image
Figure 3
figure 3

Effect of in silico Alanine Scanning mutagenesis of interfacial residues on the energy of association between the interacting domains.

Each of the interfacial residues of the 428 interfaces was mutated to Alanine in silico and the change in the interaction energies between the two domains was computed. All residues except Glycines and Alanines were mutated in this analysis. Most Ala-substitutions in the interface were found to be destabilizing. We find that the interface stability is majorly governed by large hydrophobic interactions (involving L, M, I, F, Y and W) and few polar contacts (involving R, D and E), as we observe the medians values of ΔGInteraction for these mutations are above the 1kcal/mol (Red line) threshold.

Full size image

Solvation effects and stability

To understand precisely the contribution of solvent effects on the stabilization of domains, we computed the solvation energies explicitly (See methods) and found that the solvation energies are significantly less favourable (Paired t-test: two-tailed; t = 2.262; df = 704; P = 2.41E-02) for domains of multi-domain proteins considered in isolation (−5819 ± 135.4 kT Mean ± s.e.m.) in comparison to single domain proteins (−6177 ± 165.7 kT) (Figure 4). We find that the sum of expected interaction energy (given the precise surface area of interaction) (ΔGint,(SA)) for a domain in multi-domain system (obtained from the linear fit in Supplementary Figure S2) and ΔGfolding of the same domain (−6.047 ± 1.05 kcal/mol) is comparable to the ΔGfolding for the single domain homologues (−2.156 ± 2.34 kcal/mol) (Unpaired t-test: two-tailed; t = 1.689; df = 1115; P = 9.15E-02) and the complete multi-domain chain (Unpaired t-test: two-tailed; t = 2.344; df = 1282; P = 8.15E-01) (Figure 1).

Figure 4
figure 4

Free energy of solvation of domains of multi-domain proteins and single domain homologues .

The solvation free energy computed as the difference in the reaction field energy (See methods) of the protein when the external dielectric is changed from ε = 1(air/vacuum) to ε = 80 (water) using the Delphi solvation model. The figure shows that the MD_doms have significantly higher solvation free energy in comparison to their single domain homologues (Paired t-test ; t = 2.262, df = 704, P = 0.0240).

Full size image

Residue changes in single domain proteins and their implications on domain stability

We quantified the sequence divergence among the corresponding protein pairs. We observed that there is considerable divergence (Mean sequence identity = 23.0% ±15.4) in the sequences of the domains of multi-domain proteins and their respective single domain homologues. To understand how sequence divergence of 1D_homs from the MD_doms, results in differential stability, we quantified residue changes specifically in the interface of domains of multi-domain proteins. We measured the propensities for each of the 400 possible mutation effects in these interfaces (See methods) (Figure 5a). By comparing domains of multi-domain proteins with their single domain homologues we note that most interfaces have been changed considerably with non-polar residues in the interface substituted by polar residues in the single domain homologues. We observe a clear electrostatic surface change in the domain pairs (Figure 5a: top left quadrant) for the regions corresponding to domain-domain interfaces. We further quantified the effect of these mutations by mutating in silico the non conserved interface sites of the multi-domain protein domain into the corresponding topologically equivalent residues of single domain homologues (Figure 5b). We performed swapping mutagenesis in silico by substituting residues in surface patches of the single domain homologues that are topologically equivalent to domain-domain interface of multi-domain proteins by the topologically equivalent residues in the domain-domain interface of multi-domain proteins and vice versa. We found that the swapped mutations showed no specific relationship among homologous domains (Pearson correlation: r2 = 0.02092; n = 5019; P < 0.0001); single domain homolog mutants are more destabilized (ΔΔG = 1.54 ± 3.49 Kcal/mol) than the multi-domain protein domain mutants (ΔΔG = 0.885 ± 1.76 Kcal/mol) taken out of multi-domain context (Paired t-test: two-tailed; t = 12.60; df = 5018; P = 7.34E-36) (Figure 5b).

Figure 5
figure 5

Substitution preferences in single-domain homologues of multi-domain protein domains and the effect of swapping mutagenesis on the domain stabilities of multi and single domain proteins.

(a) Structure based sequence alignment between single domain-homologues and multi-domain protein domains were generated and used to calculate the propensity of interface mutations (See methods). Observe the preferences for polar to apolar substitution frequencies in going from single-domain homologues to MD_doms. For e.g. the propensity of W to K substitutions in the interface. (b) The MD mutants represents the distribution of single point mutations (n = 5019), where each of the non-conserved interfacial residue was mutated to the corresponding residue found in the respective single domain homologue in the structure based sequence alignments. The 1D_hom mutants represents a set of single point mutations (n = 5019) on single domain proteins where the residues corresponding to the interface in the multi-domain protein homologue were mutated to those found in the MD_doms in the structure based sequence alignments. MD_mutants were significantly less destabilizing (Paired t-test; t = 12.60; df = 5018; P < 0.001) and show no relation (Spearman rank correlation test: r2 = −0.067; P = 1.50E-06) to stability effects in 1D_hom mutants.

Full size image

Selection of interfacial residues and guidelines for domain stabilization

In order to identify potential key sites in the exposed interface patch of the multi-domain protein domains that are mutable to enhance the stability of the domain, we classified the interfacial residues into stabilizing; neutral and destabilizing (See methods). This classification is based on the free energy contribution of the individual residue type towards the overall domain free energy for folding in isolation (See Supplementary Fig. S3 online). Stabilizing and destabilizing residues are those with energy contribution towards the overall domain stability are respectively higher than and lower than a 95% confidence limit of the contribution of a typical surface residue to overall stability. Contribution of neutral residues is not different from that of a typical surface residue. We computed residue conservation in single domain homologues in positions which are topologically equivalent to the interface positions of domains in multi-domain proteins. We contend that the residue conservation in the single domain homologues is dependent on the class of interfacial residues (Stabilizing/Neutral/Destabilizing) which they belong and not on the frequency of observation of these three classes. We observed a bias in the proportion of residues conserved for these classes (χ2 = 80.42, df = 2, P = 3.44E-18, n = 1574; the a priori expected proportions of stabilizing, neutral and destabilizing residues are 0.25, 0.95 and 0.25 respectively) in the corresponding single domain homologues (Figure 6a). Stabilizing residues are conserved better (29.84%; n = 94) than the neutral (14.19%; n = 1452) and destabilizing residues (12.44%; n = 225) (Figure 6a). Furthermore, binomial tests between proportion of residues conserved and their expected proportion are significant only for neutral and destabilizing residues (P < 0.0001, the expected conservation frequency for the binomial tests were computed as the 23% of total residues for each class, as the mean sequence identity between the MD_doms and 1D_homs is 23%. Figure 6a). To understand how this differential conservation has resulted we probed into the mechanisms by which a surface residue contributes to the overall stability. We measured three parameters, 1. Residue depth (Å) 2. Residue protrusion index and 3. Number of Cα—Cα contacts. We find that the stabilizing residues are located more deeply (Figure 6b) (on the average the residue depth is 2.13 Å ± 0.59) than the destabilizing residue positions (average residue depth: 1.51 Å ± 0.27). This resulted in having substantial intra-protein contacts (Figure 6d). Further the protrusion of these residues into the solvent is also considerably less, rendering locally compact structures (Figure 6c). We further tested the effect of point mutations on every interfacial residue site (See methods). We substituted each residue into 19 other residue types in silico using FoldX and computed the ΔΔG of the mutations. We found that the stabilizing residues in the locations of isolated domains equivalent to domain interface when mutated, lead to more destabilization of the domain than the neutral (Mann-Whitney U-test: U = 318000; rank sum : Stabilizing = 968900, Neutral = 22860000, P < 0.0001 ) and destabilizing residue sites (Mann-Whitney U-test: U = 3217; rank sum : Stabilizing = 40250, Destabilizing = 12400; P < 0.0001) (Figure 6e). Further the destabilizing residues when mutated have a larger chance of stabilization (% of stabilizing mutations i.e. with negative ΔΔG are 9.5, 27.28 and 55.55 when stabilizing, neutral and destabilizing residues were respectively mutated).

nsSNPs in the domain-domain interfaces

To understand the effect of single amino-acid substitutions in the interface, we analyzed the non-synonymous Single Nucleotide Polymorphism's (nsSNPs) in the interfaces of multi-domain proteins. We specifically studied human multi-domain proteins and identified nsSNP in the domain-domain interfaces (See Supplementary Table S1 online). We identified 33 nsSNPs in the domain-domain interface of the 14 human proteins from our dataset (See Supplementary Table S1 online) and studied their effect on the structure and stability of the interface (See methods). Most of the destabilizing changes (14/33; ΔΔG >1.0 Kcal) are observed where a wild type non-polar residue is mutated to a polar residue due to the base substitution (See Supplementary Table S1 online). This pattern has been observed in one third of the nsSNPs studied (See Supplementary Table S1 online). Also larger destabilizing effects are noted when those residue sites involved in ligand-binding (8/33 cases studied) are mutated.

Discussion

To be able to fold independently, the free energy for folding (ΔGfolding) for domains should be favourable. The FoldX20 computed free energies (See Methods and Supplementary fig. S4 online) showed that domains of multi-domain proteins in isolation (MD_doms) are significantly less stable in comparison to (1) full length polypeptide chains (MD_chains), (2) single domain homologues of domains in multi-domain proteins (1D_homs) and (3) the control dataset of single-domain proteins (1D_uniq) (Figure 1). This could be due lack of stabilization forces within the isolated domain, as we observe that the extent of stability of full length multi-domain proteins (MD_chains) is not very different from that of independently stable single domain homologues (1D_homs). This finding highlights the importance of inter-domain interaction energy involved in stabilizing otherwise unstable domains (MD_doms). The hydrophobic contacts predominate in the domain-domain interface (Figure 2). A few observed salt bridges are likely to provide specificity and shape complementarities required for stable interactions across the interface (Figure 2). As the predominant force driving the interface formation is hydrophobic, we expect a linear relationship between the interface surface area and the interaction energy (See Results). In silico Ala scanning can show the relative importance of a residue in forming an interface. The energy change greater than 1 kcal/mol upon Alanine mutation, indicates that the WT residue at that site contribute significantly to the overall interaction energy across the interface. The most interacting residues in (mostly hydrophobic and few polar) the domain-domain interface also show drastic destabilization (> 1kcal/mol; Figure 3) in terms of interaction free energies, when mutated to Alanine, thereby acting as hotspots. These hotspot residues stabilize the interaction and the full length proteins. In isolation the unstable nature of domains of multi-domain proteins (MD_doms), to exist independently (Figure 1), is believed to be due to lack of sufficient solvation energy21,22,23,24,25. The domains of multi-domain proteins as shown before (Figure 2) have a large number of hydrophobic residues at the domain-domain interface which are exposed in the absence of domain-domain interaction, resulting in poor solvation of these domains and high unfavourable ΔGfolding. Indeed it has been implicated in the hydrophobic collapse26 during the folding of multi-domain proteins. Supplementary Figure S1 shows an example of a multi-domain protein from our dataset highlighting the importance of hydrophobic character in the domain-domain interface. These extended flat hydrophobic surfaces between the domains, forces many multi-domain proteins to collapse and fold cooperatively. This process is driven by a quick desolvation of the interface during folding25. This is shown in our empirical free energy computations using FoldX (See Supplementary Fig. S1 online). We further corroborate our initial energy computations by computing explicit solvation energy (Figure 4). This clearly showed us that reduced solvation makes the isolated domains (MD_doms) unstable (Figure 1 and Figure 4). We hypothesize that multi-domain proteins compensate this loss of solvation energy by interacting favourably with their tethered domains, thus increasing the overall ΔGfolding. The interaction energy between two interfaces scales linearly with the interaction surface area (See Supplementary Fig. S2 online). In theory the interaction energy required to compensate the loss of solvation effects can be computed from this linear fit. This when added to the free energy of the isolated domain makes it favourable (Figure 1) and comparable to the free energy of folding of single domain proteins (Figure 1). This shows that the stabilization of an isolated unstable domain can be enhanced by incorporation of suitable residues to result in improved interaction energy. This also indicates that the interaction energy computed from the linear relationship with exposed surface area is approximately equal to the additional solvation energy of 1D_homs and the inter-domain interaction energy of the full-length multi-domain proteins (MD_chains). Moreover, this shows that the unstable nature of domains can be modulated by engineering interactions in the domain-domain interface patch.

The single domain homologues adopt the same fold and topology (i.e. belong to the same family in SCOP database) as the domains in multi-domain proteins, but have enhanced stability (Figure 1), in spite of high structural similarity (Average Dali Z-scores: 15.73±9.44; n = 705). We showed that there is considerable sequence divergence between the corresponding protein pairs and the single domain homologues have specific mutations at topologically equivalent sites corresponding to the domain-domain interface (Figure 5). These mutations make the surface more polar and amenable for solvation, thereby enhancing the overall free energy (Figure 4). Previous studies showed that polar-apolar residue substitutions in homologous proteins are less frequent than other patterns and are implicated in amyloidogenesis27. We have already showed that the single domain homologues have more favourable solvation free energy (∼ 350 kT difference) than the domains of multi-domain proteins (Figure 4). This result is a direct outcome of the observed residue changes in the interface (Figure 5a). The energetic effects of swapped mutations were more drastic on single domain homologues in comparison to the multi-domain protein mutants (Figure 5b). This shows that choice of residues (alterations and conservation) at interface sites is governed by stability contributions. Stabilizing mutations (non-polar to polar) are preferred at interface sites of domains of multi-domain proteins to aid in their independent existence (Figure 5b and Figure 6a). The single domain homologues are stable and belong to the same SCOP family as the isolated unstable multi-domain protein domain. These residue changes primarily modulate the surface electrostatics of the domains (Figure 5a) resulting in favourable solvation free energy (Figure 4). This pattern of specific residue changes observed in this analysis can aid in successful design of fragments and modulating domain stability in isolation. We propose a classification scheme to support this idea (See Supplementary Fig. S2 online). It is known from previous studies28,29 that surface residues are less conserved in proteins and conservation of non functional residues is generally attributed to constraints in the structural fold30. Stability contributes significantly to the overall fitness function and hence we would expect those sites with residue contribution to the overall stability low to be less conserved in homologous proteins9,31. The observed substitution frequencies and hence the nature of conservation at a given site would be proportional to the energetic effects of the residues' contribution towards the overall stability32,33. For a given domain-domain interfacial residue, we expect higher conservation of stabilizing residues in comparison to neutral and destabilizing residues. Here we showed that the energy based classification of interfacial residues clearly segregates the residues into those under differential selection pressure (Figure 6). The differential selection pressure is the direct outcome of the differential contribution towards stability. This resulted in different extent of conservation of the three classes. Further, effects on stability of the domains show a gradation when these three classes are mutated independently. The loss of free energy upon mutation thus is site dependent (specific class) and is more drastic for the stabilizing residues. This clearly shows that destabilization of a structure caused by a residue is not tolerated as it is mutated almost always. We also observed a trend in the residue depth, number of Cα—Cα contacts and the protrusion of a residue to be dependent on the three classes. Stabilisation is achieved by deeper residues having more contacts and less protrusion. This classification is thus shown to have a mechanistic basis. Using in silico mutagenesis we probed the plasticity of these sites to accommodate different residues and modulate stability. Most point mutations on proteins in such studies are destabilizing34. We propose that the site-specific nature of the classified residues determines the extent of stabilization/destabilization upon mutation. We found striking difference in the proportion of mutants gaining stability when destabilizing residues positions were mutated to the other 19 residues. Thus there are more chances of enhancing domain stability by selectively mutating destabilizing interfacial residues to naturally occurring substitutions found in their equivalent sites of single domain homologues.

Natural variations/mutations at the domain-domain interfaces can be studied by analysing the nsSNPs. Small variations in the amino acids at the domain-domain interface can have large functional consequences35,36. Most monogenic and several polygenic diseases have been known to be caused by nsSNPs37,38,39. The effect of nsSNP can be deleterious or neutral depending on the extent to which the variant amino acid affects the protein structure39,40 and function41,42. The severity of the nsSNP had been shown to be directly correlated with its effect on protein structure and stability43,44. Few of the observed nsSNPs (9/33) in the domain-domain interface have been studied by computational and experimental methods (See Supplementary Table S1 online)45,46,47,48. We believe that deleterious nsSNPs in the domain-domain interfaces act by breaking/altering the communication between the two interacting domains, required for proper folding/stability, thus leading to non-functional proteins followed by disease manifestation.

In conclusion, often, a domain of a multi-domain protein interacts with other domains in the same protein primarily by hydrophobic residues in the inter-domain interface, thereby stabilizing an otherwise unstable domain. Independently existing stable single domain homologues of these domains undergo selective non-polar to polar mutations which modulate the solvation and hence stabilize the fold. Here we also showed that the sub-set of residues in the inter-domain interface (self domain destabilizing residues) of multi-domain proteins can be selectively mutated to aid in the formation of stable single domain constructs with the ability to fold, which otherwise might form inclusion bodies inside cells. This is of potential importance to crystallography and NMR where routinely protein domains of large multi-domain proteins are generated to solve the domain structure in isolation. The present study suggests strategies for domain stabilization in such cases. Naturally occurring variants (nsSNPs) in the inter-domain interface often contribute to disease states by altering the extent and nature of interactions between the tethered domains.

Methods

Dataset generation: wild-type and mutant protein/domains

We obtained a set of single and two-domain monomeric proteins from PDB. Monomeric proteins have been identified by consulting PiQSi49 and PISA50. We used domain boundaries as proposed by SCOP database51. These structures were filtered using resolution cut-off of 2.5 Å and sequence identity cut-off of 40%. The SCOP families corresponding to multi-domain proteins have been used to identify single domain homologues. A control dataset of non-redundant single domain proteins from completely different families has also been formed. The summary of the datasets used in the current analysis is provided in Table 1. The lists of PDB and SCOP accession codes are provided in the Supplementary information (Datasets). All the structures have been optimized at 298K by using the optimize tool available in FoldX (version 3.0 beta)20 before free energy computations. The in silico modelling of mutations have been performed using mutate module of FoldX. Every interfacial residue has been mutated to 19 other residue types and the change in the free energy for folding between the mutant and wild-type (ΔΔGfolding) has been computed. This was then used to study the stabilizing and destabilizing mutations in the interface. In silico Alanine scanning mutagenesis has been performed for interfacial residues and the changes in interaction energy (ΔGInteraction = GInt,Ala − GInt,WT ) between the two interacting domains have been computed.

Energy computation and stabilities of wild type & mutant domains/proteins

We used an empirical effective energy function FoldX52 to compute the ΔGFolding for all the proteins. The FoldX energy function's stability module computes ΔGFolding and also provides a residue level contribution to the overall ΔGFolding . The FoldX energy function was originally fitted for the ΔΔG computations for mutations52 to be used as a modelling tool. In order to be able to compare ΔGFolding values of large datasets of proteins, we benchmarked FoldX calculations using experimental information. The free energies for folding computed using FoldX have been compared to experimentally determined free energies obtained from Protherm database53 for a set of test proteins (n = 35) (See Supplementary Information: Datasets). We filtered the Protherm database to obtain this test set of proteins, with known 3D-structure, which unfolded reversibly in a 2-state manner in the presence of chaotropic agents at 298K. Further, the FoldX calculations have been done at the same temperature and pH as the unfolding experimental conditions. Significant correlation (Spearman rank correlation: r = 0.5679; P = 4.00E-04) has been noted between the FoldX computed and experimentally determined free energy changes (See Supplementary Fig. S4 online). This encouraged us to use FoldX to compute free energy for folding for all the four datasets at 298K and pH 7.0 and also compare the distributions of ΔGFolding. Electrostatic free energy and consequently the solvation free energy for both multi-domain protein domains as well as their single-domain homologues have been computed independently (Figure 4). Solvation energy was computed by using the solution of non-linear Poisson–Boltzmann equation implemented in program Delphi54,55, as the difference between total grid energies in water and vacuum. The internal dielectric (εint ) of 4 has been used for the protein's interior and external dielectric (εext ) of 1 and 80 has been used in vacuum and solvent computations respectively.

Domain-domain interaction preferences & interaction hotspots at interfaces

Inter-domain interfaces have been recognized using the FoldX complex analysis module. To quantify interactions across the domain-domain interface, a distance based approach has been used. Residues across the interface with distances below the cut-off have been considered as interacting. The cut-off distance is computed as the sum of van der Waal's radii of interacting atoms plus 0.5Å56. The van der Waal's radii for the protein atoms were taken from Chothia (1975)57. Using this information a pair-wise interaction matrix has been computed to identify most frequent interacting pairs in the interface (Figure 2). In order to quantify the importance of these pair-wise interactions, we also computed the interaction energy between the two interacting domains for all the proteins. In silico Alanine scanning mutagenesis at all non-Ala and non-Gly interface sites has been performed to compute the strength of interaction hot spots. The interaction hot spots have been identified by an increased interaction energy ( >1 Kcal/mol destabilization) between two interacting domains on mutation into an Alanine (Figure 3).

Classification of interfacial residues

Inter-domain interfacial residues have been classified into stabilizing neutral and destabilizing. This classification is based on energy contribution towards the overall stability (ΔGfolding) of the domain taken in isolation for the free energy computations. All the interfacial residues energy contributions have been assessed using FoldX stability module. We used a residue-dependent cut-offs for classifying the interface residues into the above three classes. These cut-offs were obtained for all the twenty residue types using the control dataset (1D_uniq). The 1D_uniq (Table 1) was used to compute the free energy contribution of all the surface residues. These were then segregated into the 20 residue types. The energy values obtained have been tested for normality and fitted to normal distributions to compute mean (µ) and standard deviation (σ). This has been performed for all the 20 amino-acid types. The lower (µ − 1.98σ) and upper (µ + 1.98σ) limits of the 95% confidence interval of these distributions have been used as lower and upper cut-offs for the classification. This enabled us to overcome the bias in the FoldX computed energy-values for certain amino acids. If the interface residue energy is less than µ − 1.98σ, or more than µ + 1.98σ, then it is classified as a stabilizing residue and destabilizing residue respectively. All the residues with intermediate energy contribution are classified as neutral (See Supplementary Fig. S3 online).

Residue conservation and changes at interface

We quantified the preferred mutations occurring in the interface region by measuring the propensities of all 400 possible mutations. The structures of domains of multi-domain proteins and their respective single-domain homologues have been aligned using DALI58. The structure based sequence alignments thus generated have been used to compute the observed changes in the interface. The propensity of a particular mutation to occur in the interface of a single domain protein (Pi→j, where, i and j are topologically equivalent sites in multi-domain protein domain and single domain homologues respectively) has been measured as the ratio between the frequency of observed mutation in the interface () and the frequency of same mutation in the data set over the entire sequence ().

Pi→j > 1 indicates most preferred substitutions (Figure 5a) in the single domain homologues.

Residue micro-environment and stability

The residue depth; DPX index (defined as the distance of a non-hydrogen buried atom from its closest solvent accessible protein neighbor), residue protrusion; CX value (defined as the ratio between the external volume and the volume occupied by the protein within a sphere of 10 Å centered at every interfacial atom. This has been averaged over all the atoms for a given interfacial residue) and number of Cα—Cα contacts within a 12 Å sphere for every interfacial residue has been computed and used to quantify the local environment of any given residue. The DPX and CX indices59 have been computed using web-servers developed at ICGEBnet Protein tools.

Analysis of natural variations in the domain-domain interfaces

In-order to understand the stability effects of natural variations in domain-domain interfaces; we have analyzed non-synonymous single nucleotide polymorphisms (nsSNPs) in the human proteins (Supplementary information Table S1). We identified the nsSNPs in these proteins from the mapping database (SNPdb to PDB) LS-SNP/PDB60. We identified 33 interfacial nsSNPs in 14 human proteins. We studied the effect of nsSNP by modeling the structure of variant proteins and computing the ΔΔGfolding using the FoldX build module. We also obtained other parameters such as ligand binding, solvent accessibility and conservation of the nsSNP sites from LSSNP/PDB60.