the ability of proteins 
 to interact with 
 different partners is 
 essential to the life of the cell .
indeed , the withdrawal or damage 
 of proteins with a high number of connections in protein  protein 
 interaction ( ppi ) networks ( hubs ) is often lethal or associated with disease .
binding to multiple partners can be promoted by different 
 factors , including multiple domains , domain repeats , multiple interaction 
 sites , intrinsically disordered regions , conformational changes , highly 
 charged surfaces , post - translational modifications , and alternative 
 splicing . even if using different definitions , hub proteins have been traditionally classified in 
 two groups , according to their tendency to use the same (  date  , 
  singlish - interface ,  or  sociable  hubs ) 
 or distinct (  party  ,  multiple interface ,  
 or  stable  hubs ) interfaces to interact with different 
 partners .
promiscuity in binding sites that are reused in different 
 complexes has been often found to rely on the ability of the protein 
 to adopt different conformations according to the partner , with changes ranging from local side chain reorientations to global 
 rearrangements and to disorder - to - order transitions .
the role 
 of intrinsic disorder in hub proteins has been thoroughly 
 investigated in the past . however , there are many examples of multipartner proteins with an 
 ordered structure , such as calmodulin , ubiquitin , or ras proteins . 
 in these cases ,
a possibility is that , in the absence of intrinsic 
 disorder , promiscuity is promoted by specific dynamical properties 
 pre - existing in the protein isolated state .
indeed , in the last decades 
 an increasing number of experimental and theoretical studies have shown that the conformational changes related with a change 
 in the protein state ( binding to ligands , post - translational modifications , 
 interaction with light ) can be sampled by the protein also in the 
 initial unperturbed state .
some recent works highlighted that in specific 
 cases the formation of multiple arrangements at a protein  protein 
 interface can be assisted by conformational fluctuations in the unbound 
 state in one of the partners .
however , only a few large - scale studies of structure - based ppi networks have explicitly considered the inherent 
 flexibility of proteins and evidence of a general relationship between 
 intrinsic dynamics and binding promiscuity in ordered proteins is 
 still missing .
a first problem in the determination of a link 
 between intrinsic 
 dynamics and binding promiscuity is that , while classifications in 
 ppi networks are usually made at the protein level ( e.g. , hub and 
 nonhub ) , the regions in a protein that are involved in interactions 
 with partners ( semi - interfaces or simply interfaces in the following ) 
 can have a largely heterogeneous composition in terms of different 
 properties , including mobility , evolutionary conservation , binding promiscuity , and binding affinity .
thus , a promiscuous protein can 
 be composed of both nonpromiscuous and promiscuous residues , each 
 with a different degree of intrinsic flexibility .
moreover , the flexibility 
 of a protein can be measured in terms of different space and time 
 scales , with previous works usually focusing on either backbone motions , slow and highly collective , or side chain motions , faster and more localized . here ,
we present the first 
 large scale study of the intrinsic dynamical 
 properties of promiscuous residues , where we mapped residue - based 
 measures of conformational flexibility and binding promiscuity on 
 interfaces of a large data set of proteins .
both backbone and side 
 chain intrinsic motions were described with in silico generated conformational ensembles of the isolated proteins .
the 
 analysis revealed higher flexibility levels for promiscuous residues 
 compared to nonpromiscuous ones .
this additional flexibility was found 
 to be highly organized in correlated motions both on a local and on 
 a global scale , indicating that , when intrinsic disorder is not a 
 major factor , an ordered form of flexibility can take over to promote 
 binding promiscuity .
remarkably , the magnitude of intrinsic motions 
 in promiscuous residues showed a reduced dependence from evolutionary 
 conservation when compared with nonpromiscuous ones , providing an 
 unprecedented indication that binding affinity and promiscuity can 
 have opposite effects on residue dynamics .
the functional importance 
 of promiscuous residues was also confirmed with an investigation of 
 the distribution of single nucleotide polymorphisms ( snps ) across 
 interface regions , which suggests that residues in promiscuous positions 
 have a reduced tolerance to genetic variations , related to the necessity 
 to preserve their binding polyvalence .
a data set of proteins was generated 
 using the pisite database ( figure 1 ) .
a nonredundant list of proteins from pisite was 
 used , composed of 7739 proteins with sequence identity < 30% identified 
 by clustering the sequences of 110325 proteins from 51482 pdb entries . for a given protein
, pisite collects all the 
 pdb entries containing the protein itself or very close homologues 
 ( > 90% sequence identity ) . the original protein from the nonredundant 
 pisite list and its homologues
, 
 we will refer to sequence families as simply families and to the original 
 proteins in the nonredundant pisite list as family representatives 
 ( frs ) .
in pisite , all the pdb complexes containing a member of the 
 family are used to define the family binding properties ( number of 
 partners and binding states ) , which are then mapped onto the fr .
the 
 basic assumption is that , because of their high sequence similarity , 
 members of the same family have a very similar behavior in terms of 
 binding interactions , so that structural information on partners can 
 be transferred from the family members to the fr . in this way ,
the 
 partner annotation of fr is enriched compared to when only fr occurrences 
 in the pdb are considered , and the bias due to the incompleteness 
 of the pdb is reduced .
the 
 flowchart shows the steps followed in the generation of the s and s data sets from the nonredundant pisite 
 database .
the operations indicated in blue on the right were performed 
 on all the members of each family .
see methods in the main text for 
 a full description of the procedures involved . for each family in nonredundant pisite ( figure 1 )
, we selected the members with known uniprotkb and scop ids 
 using the pdbsws pdb / uniprot mapping and 
 the scop ids using the astral scop database ( v 1.75 ) . only the families with unambiguous scop domain classification 
 across all the members and with at least one member with known uniprotkb 
 i d were retained .
we then selected 251 families that satisfied the 
 following requirements:(1)their members have only one scop domain.(2)the sequences of the resolved structures 
 in the family cover at least 75% of the corresponding uniprotkb sequences.(3)they have at least one 
 partner with 
 known structure.(4)there 
 is at least one structure in 
 the family with no gaps in the resolved main chain .
the ungapped x - ray 
 structure with the best resolution in the family was selected as the 
 structural representative ( sr ) to be used in the simulations and structural 
 analyses .
when no crystallographic structure with a complete main 
 chain was found , an ungapped nmr structure was selected as sr if available .
their members have only one scop domain . the sequences of the resolved structures 
 in the family cover at least 75% of the corresponding uniprotkb sequences .
there 
 is at least one structure in 
 the family with no gaps in the resolved main chain .
the ungapped x - ray 
 structure with the best resolution in the family was selected as the 
 structural representative ( sr ) to be used in the simulations and structural 
 analyses .
when no crystallographic structure with a complete main 
 chain was found , an ungapped nmr structure was selected as sr if available .
the 251 families define our full data set 
 s ( supporting information ( si ) 
 table s1 ) .
each family 
 includes on average 20 members , for a total of 4917 pdb chains .
the pdb and chain ids of partners from structural complexes ( or 
  structural partners ) were extracted from pisite for 
 each family in s and mapped to uniprotkb with the pdbsws 
 mapping . for consistency with the pisite database ,
redundancies were 
 removed by clustering the partner sequences with blastclust ( sequence identity 30% , sequence coverage 
 50% ) to define the total number of structural partners ( nps ) of the family .
this was used to classify the families of 
 s as monopartner ( nps = 1 , for a total of 151 families ) and 
 multipartner ( nps  2 , for a total of 100 families ) .
each of the 
 main scop classes (  ,  , / , + ) 
 was represented in both groups with at least 14 entries .
the proteins 
 in the families have an average size of 217 ( monopartner ) and 167 
 ( multipartner ) residues , as measured from the srs .
family partners 
 identified by nonstructural detection methods were 
 extracted from the intact database using 
 the uniprotkb ids of each family member .
the total number of unique partners ( np ) of each family 
 was obtained using blastclust ( see above ) on the sequences of all 
 the partners extracted from either intact or pisite . in the 
 following , we will often refer to a family as simply a protein , 
 with the fr as its sequence representative .
unless otherwise stated , 
 structural properties and simulated conformational ensembles were 
 calculated from the sr of the family , while binding multiplicity profiles 
 and ensembles of experimental structures were built - up using all the 
 family members ( see the following and inset in figure 1 ) .
md simulations ( see below ) were performed on a subset 
 ( s ) of s composed of 6 monopartner and 
 6 multipartner 
 proteins .
the s proteins were selected using stricter 
 criteria to increase the confidence of their classification , so that 
 multipartner proteins have at least 4 structural partners , the monopartner 
 proteins have both np and nps = 1 , scop classes are equally 
 sampled by the two groups and the average size of the proteins in 
 the two groups is comparable ( si tables s1 and 
 s2 ) .
a second data set ( s ) was considered 
 to test and 
 compare the findings observed on s. s is composed of sociable proteins ( si table s3 ) , defined as having at least 3 structural 
 partners and 3 different binding states .
starting from the original list of 102 sociable proteins , 69 were selected so that the fr has a complete 
 main chain structure .
the fr was used as the structural representative 
 sr for the simulations and structural analyses . since the original 
 data set of sociable proteins was not filtered for monodomain proteins , 
 s
only a limited number of 
 proteins in the data sets ( 4 in s and 5 in s ) contained post - translational 
 modifications ( ptms ) in the pdb structure , amounting to less than 
 3% of the total .
thus , ptms were not explicitly taken into account and residues with 
 ptms were modeled using the corresponding unmodified residue .
a systematic 
 study of the effects of ptms on interface dynamics would require ad 
 hoc data sets and a standardized framework including parameters for 
 a large range of ptms .
a complete set of tools for the inclusion of 
 ptms in md simulations has been recently proposed . for each s family f
, interface residues were mapped on the sr by analyzing 
 all the pdb binary complexes a : b ( where 
 a is a member of the family and b is 
 one of its partners ) listed in the pisite database for the selected 
 family members . for consistency with pisite , the binary complexes 
 were extracted from the pdb biological units .
the popscomp method was used to determine the profiles of normalized 
 solvent accessible surface area ( sasa ) that is buried upon complex formation for each chain a sasa(i , a : b ) = ( sasa(i , a ) 
sasa(i , a : b))/ sasa(i , a ) , where sasa(i , a ) and sasa(i , a : b ) 
 are the sasas of residue i in the isolated and bound 
 a molecule , respectively .
all the sasa profiles of the family were then mapped 
 onto the sr sequence using a  sr 
 sequence alignments given by the program profit .
the resulting sasa(i , sr : b ) values were used to define the sr interfaces .
a 
 sr residue i was considered to be part of a sr : b 
 interface if upon complex formation it buries at least 10% of its 
 sasa , that is , if sasa(i , sr : b ) > 0.1 .
if i is found in interfaces 
 involving n different partners , its binding multiplicity b(i ) is set to n. the 
 binding multiplicity was used to partition the interface residues 
 into three different binding classes .
residues with b  2 were classified as multipartner ( cmulti ) , while residues with b = 1 were grouped 
 in the cmono class if belonging to monopartner 
 srs and in the cmono_in_multi class if 
 belonging to multipartner srs ( si table s4 ) .
because of differences in the definition of the interface ( distance - based 
 in pisite and sasa - based here ) and of the selection process applied 
 in the data set generation ( see data set preparation ) , the binding multiplicity used for the s data set 
 may differ in isolated positions from the one recorded in pisite . 
 to ensure full consistency in the analysis of the s data 
 set , original b profiles from pisite
there was some redundancy in the pdb biological 
 units of each family , so that very similar a : b complexes were observed .
interface redundancy was eliminated 
 by clustering the sr - mapped sasa profiles 
 according to their overlap matrix oij = ( sasai(ij ) + sasaj(ij))/(sasai + sasaj ) , where sasai is the buried sasa 
 of sr in interface i and sasai(ij ) is 
 the surface area of the sr residues involved both in interface i and j , calculated using the coordinates 
 of interface i. the clustering was performed with 
 the complete linkage method and the optimal number of clusters was 
 chosen by maximizing the average silhouette value s. an average smax of 0.84 ( 0.19 ) was obtained for the s families , 
 indicating a good partitioning of the interfaces into the clusters . 

the interfaces with the largest sasa within a given cluster were chosen 
 as representatives to generate the final nonredundant set of 695 interfaces 
 ( si table s4 ) .
the interfaces were 
 partitioned into three classes according to 
 the maximum binding multiplicity of their residues and the binding 
 class of the corresponding proteins ( si table 
 s4 ) .
consistently with the residue classification , interfaces 
 containing at least one multipartner residue where defined as multipartner 
 ( cmulti ) , while the remaining interfaces 
 were classified as cmono and cmono_in_multi if belonging to mono- or multipartner proteins , 
 respectively .
the distributions of the total and relative hydrophobic 
 interface sizes for cmulti ( average sasa 
 = 1583  , average sasaphob = 62% ) 
 and cmono ( average sasa = 1346 
  , average sasaphob = 61% ) were 
 not statistically different , while the cmono_in_multi interfaces were significantly smaller ( average sasa = 803 
  , wilcoxon p - value < 2  10 ) and 
 less hydrophobic ( average sasaphob = 55% , wilcoxon p - value < 2  10 ) .
the values of the relative hydrophobic interface size ( sasaphob ) were calculated as the relative contribution 
 to the total sasa of the interface residues classified as hydrophobic 
 in pops .
the physicochemical properties of 
 interface residues in the different binding classes were analyzed 
 in terms of propensities relative to all solvent exposed residues 
 in the corresponding protein class .
the propensity px(ci ) of property x in the binding class ci was calculated as px(ci ) = px(ci)/px(surf ) , where px(ci ) and px(surf ) are the fraction of residues with property x in ci and at the surface , 
 respectively .
a px(ci ) value > 1 indicates a 
 higher 
 abundance of residues with a given property in ci than in the rest of the surface .
solvent - exposed 
 or surface residues were defined as having a relative accessibility 
 sasa  15% in the 
 sr of the family .
the sasa values were calculated normalizing 
 the sasa of each residue by that of the corresponding amino acid type 
 x in the axa tripeptide , where x is in an extended side chain conformation selected from the dunbrack and cohen rotamer 
 library as implemented in pymol .
propensities were calculated for the amino acid 
 identity , the evolutionary conservation grade as determined by consurf , the dssp secondary 
 structure definition , the relative accessibility of the isolated protein 
 sasa , and the extent of surface burial in the complex 
 sasa . to provide 
 a measure of the uncertainty associated with this calculation and 
 to rule out the possibility that the observed propensities were biased 
 by a few cases , confidence intervals ( ci ) at 95% were calculated by 
 bootstrap resampling with 1000 replicates .
the statistical significance 
 of differences between the propensities was estimated with a student s t test on the ci .
resampling and statistical 
 analyses were performed with r. interaction hot spots were predicted using different methods , namely 
 anchor , robetta , pisa , hotpoint , and kfc2a / b .
the 
 criteria used for the hot spot definition are summarized in si table s5 .
robetta , hotpoint , and the kfc2a / b 
 methods have been specifically parametrized for hot spot prediction 
 against experimental alanine scanning data and the default criteria 
 for residue classification defined in the original publications were 
 used .
the anchor server has instead been 
 designed to scan for residues that , besides contributing significantly 
 to the binding energy , are also 
anchors ( i.e. , they 
 have a large exposed surface in the monomeric state ) .
in addition , 
 we used pisa , an established method to determine the thermodynamic 
 stability of protein assemblies and to distinguish biological interfaces 
 from artifacts of crystal packing . here , 
 the single contributions of the residues to the interface stability 
 were used to estimate their importance in the complex formation
. thresholds 
 on anchor and pisa binding free energies were selected to obtain an 
 overall hot spot fraction comparable to the other methods .
hydration 
 scores shyd were evaluated 
 for each protein in the s data set from the spatial distribution 
 of water molecules observed in md trajectories ( see below ) .
water 
 density maps g(r ) were calculated at discrete points r defined by a 0.5- spaced rectangular grid around the 
 solute .
frames saved every 0.1 ps from the last 10 ns of the trajectory 
 were superimposed to a reference using c atom positions 
 to remove the overall roto - translational motion of the protein .
the 
 number density of the water oxygen atoms was then averaged at each 
 grid point over the md frames and normalized by the bulk density evaluated 
 in the 68  shell around the solute .
the hydration sites 
 were then identified as local maxima of the density map with g(r ) > 1 and used to define the hydration score shyd as previously described .
residues with a high shyd score were 
 surrounded either by a large number of maxima or by maxima with a 
 high density .
the shyd value of each residue i of type aa was standardized by calculating 
 the ratio ( shyd  hyd)/hyd , where hyd and hyd indicate 
 the average and standard deviation of shyd calculated over residues of type aa .
ensembles of 
 conformations were generated for each sr in the s and 
 s data sets with the tconcoord method . given a starting structure , tconcoord samples alternative 
 conformations by fulfilling a set of geometrical constraints as determined 
 from the initial coordinates and interaction types ( e.g. , covalent 
 bonds , hydrogen bonds , salt bridges , or hydrophobic interactions ) . 
 under - wrapped hydrogen bonds
each sr in the data sets was first energy - minimized 
 using the opls - aa force field with 250 
 steps of the steepest descent algorithm .
molecular dynamics ( md ) simulations of the 12 srs of the s data set ( si table s2 ) were performed 
 with gromacs 4.0 .
the initial structures 
 were solvated with a cubic box of spc water molecules , setting the 
 minimal distance between the protein and the box boundaries to 12 
  .
the berendsen algorithm was employed for temperature and pressure regulation , 
 with coupling constants of 0.2 and 1 ps , respectively .
the electrostatic 
 interactions were calculated with the particle mesh ewald method , with a 14- cutoff for the direct space 
 sums , a 1.2- fft grid spacing , and a 4-order interpolation polynomial 
 for the reciprocal space sums . for van der waals interactions ,
harmonic positional restraints with a force 
 constant of 4.8 kcalmol were imposed onto 
 the protein heavy atoms and gradually reduced to 1.2 kcalmol in 80 ps , while the temperature was increased from 200 to 
 300 k at constant volume .
the system was then simulated at constant 
 temperature ( 300 k ) and pressure ( 1 bar ) for 100 ps .
after the removal 
 of harmonic restraints , 2 ns of equilibration were run in npt conditions . 

gaussian network 
 model ( gnm ) calculations were performed by representing 
 each protein in s as a network of c atoms linked by harmonic springs .
each node was assumed to fluctuate 
 according to a gaussian distribution around its equilibrium position , 
 defined by the coordinates of the starting structure .
root mean square 
 fluctuation ( rmsf ) values were derived for each c atom by inverting the kirchhoff matrix built - up using a unitary force constant for the springs and a 7 
  cutoff on the distance between c atoms .
the conformational variability within each family of s was evaluated on the ensemble of pdb structures composed by all 
 the family members as defined in data set preparation .
the ensembles were built - up by collecting either x - ray or nmr occurrences 
 as listed in pisite .
ensembles with more than one structure were found 
 for 241 proteins in s and 11 proteins in s. in the generation of the ensembles , the pdb structures were first 
 aligned with profit . only the residues occurring in all the structures 
 were retained , so that the resulting ensembles could be analyzed as 
 pseudotrajectories with standard md analysis tools .
each structure 
 was also labeled as bound or unbound according to its binding state 
 as reported in pisite .
c and side 
 chain root - mean - square fluctuations ( rmsf ) of the residues from their 
 average position were calculated for the generated tconcoord , md , 
 and pdb ensembles .
c rmsf profiles were determined 
 after removal of the overall roto - translational motion by best - fit 
 superposition of the structures to all the c atoms 
 of the starting experimental structure ( sr ) used as a reference . for 
 side chains
, the rmsf values were calculated after removal of the 
 main chain roto - translation of the single residues , using the experimental structure as a reference for the 
 best - fit superposition .
the comparison of flexibilities between 
 different proteins and types of ensemble was performed using standardized 
 c and side chain rmsf profiles z - score(i ) = ( rmsfj(i ) 
j)/j , where rmsfj(i ) is the rmsf of residue i in 
 protein j , while j and j are the average and standard 
 deviation of the rmsfj distribution .
relevant rearrangements on protein interfaces could also arise 
 from subtle conformational changes of local structures .
these changes are often difficult to detect by traditional flexibility 
 analysis and require the isolation of local dynamics from the global 
 motion . to this end , the dynamics of 
 local structures in the tconcoord and md ensembles was analyzed with 
 a fragment - based approach .
it was previously shown that local conformational 
 changes and their correlation are effectively described by means of 
 a structural alphabet ( sa ) including prototypical geometries of backbone 
 fragments . in the present study ,
this sa comprises 25 representative 
 fragments of 4 consecutive c atoms , and it was specifically 
 designed to include the most typical local structures , as well as 
 to correctly describe conformational transitions sampled by molecular 
 simulations .
each fragment in the sa is labeled with a letter representing 
 a prototypical conformational state .
therefore , any 4-residue - long 
 segment in a protein structure can be annotated with a letter .
the 
 labeling is performed by identification of the most similar sa fragment 
 in terms of root - mean - square deviation ( rmsd ) between the protein segment and the letter . 

the conformation of a protein of n residues can be 
 encoded into a structural string of length n  

a column of this alignment summarizes the conformational states sampled 
 by a protein segment within the simulation .
the correlation 
 of conformational changes in a pair of protein 
 segments ( i , j ) can be calculated 
 as normalized mutual information ( mi ) between the associated columns 
 in the alignment ( eq 1):1where ci and cj are the 
 relevant columns in the structural string alignment , i(ci;cj ) is the mi , h(ci , cj ) is the joint entropy , and (ci;cj ) is the expected finite size 
 error .
it was previously demonstrated 
 that local correlated motions are instrumental for allosteric signal 
 transmission and may be involved in conformational changes of interacting 
 interfaces . to this end ,
the distribution 
 of statistically significant local correlations was compared for surface 
 and interface fragments in the three binding classes .
statistically 
 significant correlations were identified as previously reported . to analyze the flexibility of single 4-residue 
 fragments , a fragment 
 rmsf was calculated by defining n3 sliding 
 windows or fragments of 4 adjacent c atoms . for 
 all the structures in the ensembles ,
each fragment was superimposed 
 onto the reference experimental structure sr independently from the 
 rest of the protein , to remove local roto - translational motions .
the 
 fragment rmsf was then calculated as the quadratic mean of the rmsf 
 values of each c within the window .
the correlation between single interfaces and global 
 motions was 
 investigated with the functional mode analysis ( fma ) . given a functionally relevant property , fma
can be used 
 to identify the protein collective motion that is maximally correlated 
 with it ( maximally correlated motion or mcm ) .
this is usually expressed 
 as a linear combination of principal components ( pcs ) derived from 
 a principal component analysis ( pca ) of 
 the system trajectory .
we selected the fma approach where the coefficients 
 of the linear expansion i are determined 
 by maximizing the pearson correlation coefficient between the time 
 evolution of the functional property f(t ) and the projection of the trajectory onto the mcm . in this case , 
 the variance of f during the dynamics var(f ) can be approximated as var(f )  i=1,l i var(pci ) , where pci is the projection of the trajectory onto the ith pc and l is the number of pcs considered 
 in the mcm expansion .
the relative contribution of pci to var(f ) was thus evaluated as 
 pvari = i var(pci)/var(f ) .
the interface radius of gyration rg was considered 
 in 
 this study as a functional quantity related to the overall shape of 
 the interface .
rg was calculated over the c atoms of the interface residues for each nonredundant interface 
 in s. using an alternative property to describe the 
 interface geometry , namely the distance rmsd calculated over all the 
 possible pairs of c atoms in the interface , did not affect the conclusions described in results .
the mcm was expanded in the essential space 
 of all c atoms , composed of the first l pcs accounting for the 90% of the total c fluctuation . 

an average essential space size of 15 ( 7 ) pcs was observed in 
 the whole data set .
the first 450 structures of the tconcoord ensemble 
 were used for the pca analysis and as training set for the construction 
 of the linear model .
the remaining 50 structures were used for cross - validation . on average , the optimized pearson correlation 
 coefficient between the mcm and the rg was 0.81 ( 0.20 ) 
 and 0.79 ( 0.23 ) for the training and cross - correlation sets , 
 respectively .
this supports the validity of the linear model for the 
 mcm and rules out the possibility of overfitting in the determination 
 of the mcm .
since the motion along the mcm can be restricted in the 
 underlying energy landscape , while the 
 pcs in the essential space represent the directions along which the 
 protein motion has the largest amplitude , we also analyzed the single 
 pc contributions pvari ( see above ) to 
 the overall rg variance .
the number npc20 of distinct pcs with pvar  20% ( deviating from the average 
 value by 1.6  ) was determined .
the 
 comparison of the conformational spaces sampled by the md , 
 tconcoord and pdb ensembles was performed by calculating the normalized 
 overlap of the c covariance 
 matrices as2where 
 a and b are the matrices to compare , 
 tr is the trace operator and d(a , b ) is the matrix difference : d(a , b ) 
 = ( tr[(a  b ) ] ) .
the overlap ranges from 0 ( no overlap ) to 1 ( identical 
 matrices ) .
the cumulative overlap between two sets of pcs was calculated 
 as the average of the squared inner products between all the possible 
 pairs of pcs from the two sets .
the cumulative 
 overlap is 1 if the space spanned by the two sets is identical .
estimates of per  residue configurational entropy values 
 in the md ensembles were obtained using the heuristic schlitter s 

formula3where kb is the 
 boltzmann s constant , t is the temperature , e is the euler s number ,  is the plank s 
 constant divided by 2 , m is a diagonal matrix with the atomic 
 masses , and  is the covariance matrix .
entropy contributions 
 from the overall protein roto - translation were removed by superimposing 
 each frame in the md trajectories to the reference experimental structure 
 using all the c atoms . for each residue ,
the covariance 
 matrix  was then calculated considering all the atoms in the 
 residue or only the main chain atoms . in the all - atom case , to take 
 into account differences in the side chain length , entropy values 
 were normalized by the total number of atoms in the residue and multiplied 
 by 10 ( average number of atoms in a residue considering heavy atoms 
 and polar hydrogen atoms ) .
it is to be noted that the schlitter s 
 formula is an approximation to the upper bound of the entropy .
moreover , 
 the decomposition into residue contributions is not exact . indeed , 
 when applying eq 3 to single residues the correlation 
 between the residue and the rest of the protein is not explicitly 
 taken into account .
flexibility 
 analyses were performed on all frames of the tconcoord 
 and pdb ensembles and on md frames saved every 1 ps . most of the 
 findings presented in results are based on 
 the analysis of tconcoord conformational ensembles .
the tconcoord 
 method can be considered as a relatively fast method to explore the 
 conformational landscape of a protein . as explained above
, tconcoord 
 generates alternative conformations of a protein by satisfying the 
 constraints derived from a starting structure .
an all atom representation 
 of the protein is used , and even if energetic terms are not directly 
 involved in the generation of structures , the parameters defining 
 the upper and lower bounds of the constraints depend on the type of 
 interaction represented by each constraint .
subsequently , they were modified and combined with an estimate 
 of hydrogen bond stability to increase the portion of conformational 
 space explored by the ensemble . compared 
 to md simulations , tconcoord lacks an explicit representation of solvent 
 and long - range interactions , and it does not provide direct time or 
 energetic information on the simulated processes .
however , it is not 
 affected by convergence issues , and it can rapidly cover a large part 
 of the accessible conformational space . indeed , tconcoord ensembles 
 have been shown to be in general representative of the structural 
 variability of a protein on the basis both of md simulations and experimental 
 data .
in particular , they 
 can be effectively used as predictors of the relative flexibility 
 of protein residues , to distinguish rigid regions from flexible ones . 

moreover , large conformational transitions between different states 
 of a protein ( e.g. , from the open unbound to a close bound - like conformation 
 of calmodulin ) can be successfully reproduced 
 starting from a single state . to test the dependence of the 
 main results of the paper on the specific method used for the generation 
 of the ensembles
, we considered three more sources of information 
 on the conformational variability of the proteins .
first , we performed 
 short ( 40 ns ) equilibrium md simulations in explicit water for 12 
 proteins ( s data set ) selected from s ( see data set preparation ) .
these proteins were selected 
 so that the main conclusions can 
 be tested on them ; thus , ( 1 ) they are representative of the variability 
 of the original data set in terms of binding properties and fold , 
 ( 2 ) their partner annotation is the most reliable among the s proteins , and ( 3 ) they do not have a bias in terms of 
 size or secondary structure composition . with respect to tconcoord , 
 ensembles from short md simulations can be expected to provide a more 
 accurate representation of the thermal fluctuations around the starting 
 state , but they are less likely to sample the conformational space 
 farther from it , particularly when high - energy barriers are involved . as a third method , we used the gnm method to calculate equilibrium fluctuations around the sr structures from 
 the complete s data set .
elastic network models are 
 often used when large data sets are involved because of their reduced 
 computational cost . differently from tconcoord and md
, the gnm adopts a reduced , coarse - grained 
 representation of the protein , where usually only c atoms are considered .
their interactions are described by harmonic 
 potentials , so that a gaussian distribution of the fluctuations around 
 the equilibrium positions is assumed . despite its simplicity , the 
 gnm has been found to accurately describe the protein global motions 
 and in particular the equilibrium fluctuations in the neighborhood 
 of the starting structure . at last
, the conformational 
 variability actually observed in experimental 
 structures was analyzed by collecting all the pdb occurrences of a 
 given protein or of close homologues ( see generation 
 of conformational ensembles ) .
these collections of pdb structures 
 ( pdb ensembles ) can contain useful information on the protein conformational 
 variability .
indeed , flexibility indices 
 from pdb ensembles have been compared in the past with rmsf values 
 from simulations .
it is to be noted that pdb ensembles 
 can contain structures solved with different techniques ( x - ray and 
 nmr ) , in different experimental conditions and in different binding 
 states .
thus , in addition to conformational changes related to the 
 protein intrinsic flexibility , pdb ensembles can include changes induced 
 by external factors ( e.g. , binding to other proteins or ligands ) or 
 due to the different experimental conditions .
moreover , they can be biased by experimental 
 errors and by the fact that the pdb covers only a fraction of the 
 known states of a protein . for each of the described ensembles , 
 it is possible to calculate 
 indices measuring the extent of the structural change sampled by each 
 residue within the ensemble ( analysis of flexibility ) . as explained above ,
flexibility indices from different ensembles 
 include different types of contributions and derive from sampling 
 of conformations on different space and time scales . in particular , 
 the absolute magnitude of the structural changes observed in the selected 
 ensembles can be expected to be different .
for example , the comparison 
 of the pairwise rmsd distributions from tconcoord and md ensembles 
 of s proteins ( si table s7 ) suggests that on average larger portions of the conformational 
 space are sampled by tconcoord than by md simulations .
the relationship 
 between binding promiscuity and flexibility indices 
 derived from the different ensembles will be discussed in detail in 
 the results . however , the findings presented 
 in this paper depend mostly on the relative order of flexibility of 
 the residues in a protein and on the shape of the collective modes 
 of motion along which the protein is most free to move . hence , it 
 is useful to briefly report the direct comparison of the order of 
 flexibility predicted by the different ensembles and the calculation 
 of the overlap between the pcs extracted from them .
the calculation 
 of the correlation coefficients between the rmsf 
 profiles shows a generally good agreement among the tconcoord , md , 
 and gnm ensembles in the s database ( si table s8 ) , with average correlations ranging from 0.69 
 ( md / gnm ) to 0.75 ( c tcon / md ) .
these values are in 
 line with previous comparisons performed on different proteins . a good agreement between tconcoord and md simulations
was also found 
 for side chains fluctuations , with an average correlation coefficient 
 of 0.71 .
moreover , when the comparison between gnm and tconcoord c rmsf profiles is extended from the s to 
 the s data set , similar values of average correlation 
 coefficients are obtained ( 0.70 for s and 0.75 for s ) , indicating that s is not biased toward 
 cases with high correlation values .
as expected considering the higher 
 heterogeneity of pdb ensembles , smaller values were obtained for the 
 average correlation values between simulated and pdb ensembles , ranging 
 from 0.47 for gnm / pdb in s to 0.63 for md / pdb side 
 chain profiles in s. the conformational spaces 
 sampled by the tconcoord and md ensembles 
 of s were also compared ( si table 
 s9 ) by calculating:(a)the normalized overlap between the 
 overall covariance matrices of c atoms ( covmatover , 
 see eq 2 ) .
this value depends on the similarity 
 both of the shape of the collective modes of motion ( eigenvectors 
 of the covariance matrix ) and of the amplitude of the fluctuations 
 along these motions ( eigenvalues).(b)the cumulative overlap between the 
 essential subspaces spanned by the first 15 pcs ( cumover15 ) .
a value 
 of 15 was chosen for cumover15 since it is the average size of the 
 essential space used for the fma analysis of tconcoord ensembles ( analysis of flexibility section ) .
this overlap 
 index depends uniquely on the covariance matrix eigenvectors.(c)the maximum inner product 
 between 
 selected pairs of pcs from the compared ensembles ( maxinpr ) .
maxinpr 
 was calculated to measure how well a pc from the tconcoord essential 
 space can reproduce any of the first 3 pcs from md . similarly to cumover15 , 
 this index does not depend on the covariance matrix eigenvalues .
the normalized overlap between the 
 overall covariance matrices of c atoms ( covmatover , 
 see eq 2 ) .
this value depends on the similarity 
 both of the shape of the collective modes of motion ( eigenvectors 
 of the covariance matrix ) and of the amplitude of the fluctuations 
 along these motions ( eigenvalues ) .
the cumulative overlap between the 
 essential subspaces spanned by the first 15 pcs ( cumover15 ) .
a value 
 of 15 was chosen for cumover15 since it is the average size of the 
 essential space used for the fma analysis of tconcoord ensembles ( analysis of flexibility section ) .
this overlap 
 index depends uniquely on the covariance matrix eigenvectors . the maximum inner product 
 between 
 selected pairs of pcs from the compared ensembles ( maxinpr ) .
maxinpr 
 was calculated to measure how well a pc from the tconcoord essential 
 space can reproduce any of the first 3 pcs from md . similarly to cumover15 , 
 this index does not depend on the covariance matrix eigenvalues .
the comparison of tconcoord and md covariance 
 matrices ( covmatover ) 
 indicates a partial overlap between the two types of ensembles , with 
 an average value of 0.339 ( si table s9 ) . 

the observation of higher values for cumover15 ( average = 0.474 ) and 
 maxinpr ( average = 0.566 ) indicates a better agreement in the shape 
 of the collective motions ( defined by the relative amplitude and direction 
 of c motions ) than in the absolute amplitude of 
 the fluctuations along them .
the maxinpr index 
 was used also to compare the first 3 pcs from 
 the pdb ensembles with the pcs from the essential space of tconcoord 
 and md ensembles .
smaller values were obtained for both the pdb / tcon 
 ( average = 0.421 ) and the pdb / md ( average = 0.373 ) pairs compared 
 to the md / pdb ones ( average = 0.566 ) , indicating that collective motions 
 in the simulated ensembles are more similar to each other than to 
 the principal modes of structural change observed in the pdb ensembles . 

the larger values found in general for pdb / tcon pairs of pcs than 
 for pdb / md ones indicate that pdb pcs are better reproduced by tconcoord 
 ensembles for the set of proteins studied here .
this might be related 
 to a more complete coverage of large conformational changes by tconcoord 
 than by md simulations ( si table s7 ) .
human homologues of proteins in the 
 s and s data set were identified running 
 ncbi - blast v.2.2.26 + ( cutoff e - value 
 = 1  10 ) against a self - compiled library 
 composed of 2583 human proteins with solved or homologous structures , 
 annotated partner interactions and snp mapping .
information on nonsynonymous 
 snps was retrieved from the dbsnp database ( build 315 , ftp://ftp.ncbi.nlm.nih.gov/snp/organisms/human_9606/database/b135_archive/organism_data/b135_snpcontiglocusid_37_3.bcp ) , while disease - related snps ( dissnps ) 
 were extracted from the online mendelian inheritance in man ( omim ) 
 database ( www.omim.org/downloads ) ; therefore , only genetically inherited diseases variations 
 are considered here .
snps and dissnps were mapped back on the original 
 sequences of the s and s proteins using 
 the blast sequence alignment .
this procedure produced snp annotation 
 for 38 proteins of s and 25 proteins of s. per - protein propensity values were calculated for each binding 
 class and for each protein with the formulas described in interface analysis , using as reference the snps 
 and dissnps abundance found at the surface of the single proteins .
a data set of 251 monodomain proteins ( s ) was extracted 
 from the pdb and partitioned into 151 monopartner and 100 multipartner 
 proteins using the pisite database ( methods ) .
the composition of the data set in terms 
 of 7 general protein function categories ( si figure 
 s1 ) was obtained using a functional annotation of scop superfamilies .
as expected , monopartner proteins ( cyan ) showed an enrichment in the metabolism 
 and general categories , while multipartner proteins ( magenta ) were 
 particularly rich in the categories related to extra- and intracellular 
 processes , information , and regulation .
interface residues for 
 each protein in s were extracted 
 from its pdb complexes and partitioned into binding classes according 
 to their binding multiplicity b ( methods ) . in particular , two classes cmono and
cmono_in_multi were defined 
 for monopartner residues ( b = 1 ) belonging to monopartner 
 and multipartner proteins , respectively .
examples of both types of 
 residues are given in figure 2a , left ( cmono ) and middle ( cmono_in_multi ) panels .
residues with b  2 were assigned 
 to the cmulti class ( figure 2a , right panel ) .
a total of 12  622 interfaces residues 
 were found , with cmono , cmono_in_multi , and cmulti accounting 
 for 54 , 36 , and 10% , respectively , of the overall population ( si table s4 ) .
the results obtained on the s data set were compared with those from a second smaller 
 data set s , a collection of highly promiscuous proteins 
 ( methods ) .
as for s , the 4690 
 interface residues of s were partitioned into cmono_in_multi ( 59% ) and cmulti ( 41% ) residues ( si table s4 ) .
physicochemical properties of interface residues in the cmono , cmono_in_multi , and cmulti binding classes . (
a ) examples 
 of cmono ( left ) , cmono_in_multi ( middle ) , and cmulti ( right ) interfaces .
the residues involved in each interface are 
 color - mapped onto the surface in cyan ( cmono ) , blue ( cmono_in_multi ) , and orange 
 ( cmulti ) .
selected partners are represented 
 as cartoon . left : pdb i d 1x8d , chains
( b  e ) different types of calculated propensities 
 relative to the surface of cmono ( cyan ) , cmono_in_multi ( blue ) , and cmulti ( orange ) residues of s. error bars 
 and significance levels were estimated with bootstrap resampling . 

stars are drawn above cmono_in_multi and cmulti distributions , indicating the significance 
 levels ( * * * p - value < 0.001 , * * 0.001  p - value < 0.01 , 
 * 0.01  p - value < 0.05 ) of the comparison with cmono ( cyan ) and cmono_in_multi ( blue ) . ( b )
the conservation is expressed as consurf conservation grade ( gcons , ranging from 1 , less conserved , to 9 , 
 most conserved ) .
the three ss groups collect positions annotated 
 as helix ( h =  h  +  g  +  i  
 in the dssp dictionary ) , strand ( e =  e  +  b  ) 
 and loop ( l = blank +  s  +  t  ) .
three levels of sasa are considered , with sasa < 20% , 20%  sasa < 50% , and sasa  50% .
( f ) fraction 
 of s interface residues with a small ( < 20% ) , medium 
 ( 2050% ) or large ( 50% ) normalized buried sasa ( sasa ) .
the maximum sasa value is used for residues 
 involved in more than one complex .
( g ) standardized hydration score shyd in the s binding classes . simulated ensembles of accessible 
 conformations were generated 
 for each isolated protein in s and s using 
 the tconcoord method ( methods ) . to test the 
 dependence of results from the method used for the generation of the 
 ensembles , we performed also gnm calculations on all s proteins and molecular dynamics ( md ) simulations on selected cases 
 ( s data set , methods ) .
the intrinsic 
 flexibility of the residues was described measuring the rmsf from 
 average positions either of the c or of the side 
 chain atoms .
correlations in residue motions were calculated using 
 either c cartesian coordinates , providing collective 
 motions , or the fragment encoding from 
 a structural alphabet ( sa ) , providing local correlated motions ( methods ) . the same 
 analyses were performed on the ensemble of experimental structures 
 generated by collecting all the occurrences of a given protein in 
 the pdb ( methods ) .
a first characterization 
 of the physicochemical properties of interface residues will be followed 
 by the presentation of the central results on intrinsic flexibility . 
 in particular , the distribution of flexibility among the different 
 binding classes , its dependence from evolutionary conservation , and 
 the relationship between correlated motions and interface shape modulation 
 will be discussed .
the results obtained from the simulated ensembles 
 will then be compared with the conformational variability found in 
 the experimental structures .
the findings from the investigation of 
 the relationship between binding promiscuity and snps occurrence will 
 be then introduced .
the section will be closed by a case study exemplifying 
 the general conclusions . in this section
, we will characterize the residues 
 in the three 
 binding classes according to different physicochemical properties . 

in particular , amino acid identity , evolutionary conservation , solvent 
 accessibility , and secondary structure were determined and compared 
 with results from different data sets .
the overall interface composition 
 showed a previously documented enrichment with respect to the surface in large aromatic ( f , w , 
 y ) , hydrophobic ( m , a , i , l , v ) and specific polar / charged ( c , q , 
 r ) amino acids ( figure 2b ) . when compared to 
 the cmono and cmono_in_multi classes ( cyan and blue ) , cmulti residues 
 ( orange ) were found to be richer in the polar / charged q , d , and r 
 and , less significantly , in the aromatic f and y amino acids . on the 
 other side
, cmulti turned out to be particularly 
 poor in hydrophobic amino acids ( m , a , i , l , v ) .
these findings are 
 consistent with the more pronounced polar character generally found 
 in hub interfaces . moreover , 
 the polyvalence of the q amino acid ( which can be either an hydrogen 
 bond donor or acceptor ) and the long or bulky side chains of r , y , 
 and f are particularly suitable to adapt to different interfaces with 
 different local arrangements .
the analysis on s ( si figure s2a ) confirmed the cmulti enrichment in 
 q , r , f , and y , together with a preference , specific to the sociable cmulti residues , for the i and m amino acids . as expected , the conservation propensity pcons indicated that interface regions have a 
 higher abundance of conserved residues than the overall surface ( figure 2c ) .
indeed , pcons was 
 > 0 for consurf conservation grades 
 ( gcons ) > 5 in all the binding classes .
interestingly , 
 the cmulti group , while poorer in residues 
 with intermediate 
 conservation ( gcons = 7 , 8) was found 
 to be richer in highly conserved residues ( gcons = 9 ) than the cmono and cmono_in_multi groups , both for the s ( figure 2c ) and the s ( si figure s2b ) data set .
no large deviations 
 from the average surface values were observed for the abundance of 
 secondary structure elements in the three interface residue classes 
 of s ( pss  1 in 
 figure 2d ) .
however , when compared with the 
 two cmono ( cyan ) and cmono_in_multi ( blue ) classes , cmulti residues showed a tendency to be richer in loops and poorer in strands . 
 a similar preference for loops has been found for 
these regions are conceptually close , even if differently defined , 
 to our multipartner residues .
no significant differences in the secondary 
 structure composition were found between cmono_in_multi ( green ) and cmulti ( red ) residues in 
 the s data set ( si figure s2c ) , which both showed a general enrichment in -helices with 
 respect to the surface .
psasa ( figure 2e ) gives the propensity of an interface residue 
 to expose a small 
 ( < 20% ) , medium ( 2050% ) , or large ( > 50% ) fraction of 
 its 
 solvent accessible surface area ( sasa , see methods ) , as measured on the isolated protein .
interestingly , very exposed 
 residues ( > 50% ) were found to be particularly frequent in the cmulti group , as compared with both the cmono and cmono_in_multi groups , either for the s ( figure 2e ) and the s ( si figure 
 s2d ) data set .
similarly , if the fraction of sasa buried upon 
 complexation is considered ( figure 2f ) , a higher 
 proportion of residues with a large relative sasa ( > 50% ) was found for cmulti residues 
 than cmono and cmono_in_multi ones .
an analogous propensity for burying large 
 portions of sasa has been found in overlapping regions of date hub 
 proteins .
for the proteins in the 
 s subset of s , it was possible to evaluate 
 the hydrophilicity ( shyd , methods ) of the residues 
 from the distribution of water molecules around the solute in the 
 md simulations . to eliminate the dependence from the amino acid identity , 

each shyd value was standardized against 
 the shyd distributions of the corresponding 
 amino acid type . even if a reduced statistical significance was observed , 
 probably due to the smaller size of the s data set , the 
 comparison of the shyd z - scores showed 
 that monopartner residues in multipartner proteins tend to be less 
 hydrated ( figure 2 g ) than the other two classes . 

indeed , the presence of water molecules has been observed in pdb complexes 
 involving multipartner interfaces .
our 
 findings would indicate that the tendency to strongly coordinate water 
 molecules is mainly due to the promiscuous part of the interface . to summarize these results
, we found that multipartner residues 
 have characteristic physicochemical features that distinguish them 
 from monopartner residues belonging either to mono- or multipartner 
 proteins .
in particular , they tend to be richer in specific charged / polar 
 ( r , q , d ) and aromatic ( y , f ) amino acids and in loops , and to be 
 more conserved , more solvent exposed in the isolated protein , and 
 more buried in the complex . moreover , within multipartner proteins , 
 multipartner residues seem to be preferentially hydrated compared 
 to monopartner ones . in the following sections , we will analyze the 
 correlation between intrinsic conformational flexibility and binding 
 multiplicity .
we measured the intrinsic flexibility of s residues in terms of rmsf z - score values of either c or side chain atoms .
the distribution of flexibility values observed 
 for cmulti residues in tconcoord ensembles 
 was compared with that of cmono and cmono_in_multi residues ( averages are reported 
 in figure 3a and b , light colors ) .
all the 
 pairwise comparisons indicated a significant difference between the 
 distributions ( p - values < 0.01 ) , with an average flexibility order cmulti > cmono > cmono_in_multi for both c ( figure 3a ) and side chain ( figure 3b ) rmsf values .
a similar analysis on the s data set 
 ( si figure s3a and b ) confirmed that multipartner 
 residues ( red ) in sociable proteins are on average more flexible than 
 the monopartner ones ( green ) .
conformational flexibility of interface residues 
 in s and s. ( a / b ) average c ( a ) and side 
 chain ( b ) rmsf z - scores ( tconcoord ) calculated over cmono ( cyan ) , cmono_in_multi ( blue ) , and cmulti ( orange ) residues 
 in the s data set .
the significance levels from pairwise 
 wilcoxon comparison tests are reported with a star code ( see figure 2 legend ) .
stars are drawn above cmono_in_multi and cmulti bars 
 indicating the significance levels of the comparison with cmono ( cyan ) and cmono_in_multi ( blue ) .
averages and significance levels calculated considering 
 only the residues with the highest consurf conservation grade ( gcons = 9 ) are also reported in dark colors . 

( c ) dependence of tconcoord average c rmsf z - scores 
 ( dots ) from evolutionary conservation for cmono ( cyan ) , cmono_in_multi ( blue ) , cmulti ( orange ) residues .
residues are partitioned 
 into 9 groups according to their consurf conservation grade ( gcons ) .
( d / e ) average c rmsf z - scores calculated from tconcoord ( d ) and md ( e ) ensembles 
 over highly conserved ( gcons = 9 ) cmono ( cyan ) , cmono_in_multi ( blue ) , and cmulti ( orange ) residues 
 in the s data set .
the significance levels from pairwise 
 wilcoxon comparison tests are reported with a star code ( see figure 2 legend ) .
the intrinsic flexibility of interface residues in the different 
 classes was compared also using alternative conformational ensembles 
 ( methods ) . in particular ,
equilibrium fluctuations 
 around the native structure were evaluated using full - atom md simulations 
 and coarse - grained gnm calculations .
remarkably , multipartner residues 
 showed the largest average mobility among the binding classes for 
 all the considered ensembles and flexibility indices ( c , side chain and fragment rmsf ) , in both the complete s data set and the s subset ( si figure 
 s4 ) .
also , the relative flexibility order of the cmono and cmono_in_multi residue 
 groups was generally preserved , except for c flexibility 
 in s. in many cases , similar significance levels in the 
 distribution comparison were obtained from the different ensembles . 

this indicates that , despite the differences in the methods used for 
 the generation of the ensembles , in all of them , multipartner residues 
 have flexibility properties distinct from the monopartner ones . to assess the possible impact of these flexibility differences 
 on binding energetics , the configurational entropy of single residues 
 was estimated with the schlitter s formula ( methods ) applied 
 to the md ensembles of multipartner proteins in s ( si table s10 ) . the per - residue entropy term ts 
 ( where t = 300 k and s is calculated according to 
 eq 3 ) of multipartner residues was on average 
 higher than monopartner ones by 1.04 kcal / mol for the whole residue 
 ( considering an average residue size of 10 atoms ) and 0.44 
 kcal / mol for the residue main chain , with maximum differences of 1.52 
 and 0.65 kcal / mol , respectively ( si table s10 ) .
it has to be noted that the schlitter s formula gives only 
 an approximate estimate of the entropic term ( methods ) .
moreover , accurate measures of the impact of flexibility on the 
 binding free energy would require that protein partners are explicitly 
 taken into account in the calculation .
a further measure of 
 the relationship between flexibility and binding 
 was obtained by calculating the correlation coefficient between profiles 
 of tconcoord rmsf and binding multiplicity b for 
 each protein in s ( si figure s5 ) .
the observation of large correlation coefficients in this analysis 
 would require not only a difference in flexibility between mono- and 
 multipartner residues but also that , within multipartner residues , 
 higher b values correspond to higher rmsf values . 

indeed , 
 small average correlations were found for both mono- and multipartner 
 proteins , using either c ( si 
 figure s5a ) or side chain ( si figure s5b ) rmsf values .
however , the distribution of multipartner proteins 
 ( magenta ) was shifted toward significantly higher values than monopartner 
 ones ( cyan ) in both cases , with some multipartner proteins showing 
 correlation coefficients as high as 0.6 .
the presented results 
 show that , compared to monopartner residues , 
 multipartner residues have an average 
excess of intrinsic 
 flexibility , which could be used by the residues to adapt to the different 
 environments provided by the different partners when binding occurs . 
 in the subsequent section
, we will further refine this analysis by 
 focusing on the interface residues that are most important for the 
 interaction with the partner .
previous studies highlighted an 
 increased rigidity for hot spots and in general for evolutionary conserved interface 
 residues , which have been found to prefer preorganized bound - like 
 conformation even in the unbound state .
thus , we analyzed the relationship between evolutionary conservation 
 and conformational flexibility in our data set , to highlight possible 
 differences in the behavior of the three residue binding classes .
a small but significant anticorrelation between flexibility from 
 tconcoord ensembles and conservation was found by plotting the average 
 c rmsf z - scores of interface residues in s against the consurf conservation grade value gcons ( figure 3c ) .
a best - fit linear 
 regression showed that this dependence is slightly more pronounced 
 for the cmono class ( cyan ) , producing 
 an increased difference between cmono and cmulti c rmsf distributions 
 when only the most conserved residues are compared ( dark cyan and 
 dark orange bars in figure 3a ) .
the same slope 
 was found for the cmono_in_multi ( blue ) 
 and cmulti ( orange ) linear models ( figure 3c ) , so that cmono_in_multi average c flexibility was consistently smaller 
 than cmulti for all the gcons values .
similar findings , but with reduced differences 
 between the residue binding classes , were observed when considering 
 side chain rmsf values ( figure 3b and si s3f ) . to check
if these results are 
 affected by the method used to generate 
 the conformational ensembles , we compared tconcoord and md rmsf distributions 
 for the most conserved residues in the three binding classes of the 
 s subset ( figure 3d and e ) .
a 
 fully consistent picture was obtained from the two methods , confirming 
 the higher average flexibility of the most conserved multipartner 
 residues ( orange ) with respect to the monopartner ones ( cyan and blue ) . 
 moreover ,
the generality of the results was checked on the s data set ( si figure s3c / e ) , where cmulti residues ( red ) showed an even smaller 
 dependence of c rmsf from the conservation grade 
 than all the other classes .
we further restricted our analysis 
 to hot spot interface residues , 
 that is , the residues that contribute most to the binding energy upon 
 complex formation .
owing to the limited availability 
 of experimental information on hot spots , especially for multipartner 
 proteins , we identified candidate hot 
 spots using different prediction methods on the nonredundant set of 
 interfaces of s. for each protein and each method , 
 a given residue was classified as hot spot if it satisfies the method 
 criteria ( methods and si table s5 ) in at least one of the interfaces in which it 
 is involved . in agreement with previous findings , the majority of multipartner hot spots ( from 79 to 88% depending 
 on the specific method )
was predicted as such only in one interface 
 ( si table s11 ) , indicating that hot spots 
 are partner - specific .
hot spot populations 
 generated by different predictors had only a partial overlap ( si figure s6 ) , as expected from the differences 
 in the strategies adopted by the different methods ( methods ) . in spite of this , the comparison of the conformational 
 flexibility of mono- and multipartner hotspots produced surprisingly 
 consistent results ( figure 4 ) .
indeed , the 
 highest average tconcoord flexibility was observed for cmulti hot spots ( orange ) in almost all the cases for both 
 c atoms ( figure 4a ) and side 
 chains ( figure 4b ) .
average c ( a ) and side chain ( b ) rmsf 
 z - scores ( tconcoord ) of predicted hotspots 
 in the three binding classes cmono ( cyan ) , cmono_in_multi ( blue ) , and cmulti ( orange ) of the s data set .
the significance 
 levels from pairwise wilcoxon comparison tests between the rmsf distributions 
 are reported with a star code ( see figure 2 legend ) . for a given hotspot prediction method
, stars are drawn 
 above cmono_in_multi and cmulti bars indicating the significance levels of the comparison 
 with cmono ( cyan ) and cmono_in_multi ( blue ) .
the results presented in this section show that binding promiscuity 
 has a counteracting effect on the loss of flexibility expected for 
 more conserved positions . as a consequence ,
the preference of promiscuous 
 residues for a higher mobility is even more pronounced when considering 
 only the residues that are more likely to be determinant for the binding 
 to the partner . in the previous sections , we analyzed 
 the flexibility of residues in terms of their equilibrium fluctuations 
 from the average position in the simulated ensembles . here , we extend 
 this investigation to the correlation of these motions , to highlight 
 possible differences between monopartner and multipartner interfaces . 
 in particular , we aimed at detecting to which extent global motions 
 modulate the interface dynamics in the different classes . for 
 each protein in s , a set of nonredundant interfaces 

each interface was then 
 assigned to a binding class cmono , cmono_in_multi , or cmulti according to the maximum binding multiplicity of its residues ( methods and si table s4 ) .
the collective motions of each s protein were extracted 
 by a principal component analysis ( pca ) of its tconcoord conformational 
 ensemble . the first l principal components ( pcs ) 
 accounting for the 90% of the overall fluctuation of c atoms were selected . while other choices would be possible to represent 
 collective motions , using pcs ensured that the motions considered 
 in the analysis were allowed by the underlying energy landscape .
the 
 identification of the pcs mostly correlated with each interface was 
 performed through a functional mode analysis ( fma ) , using as functional property the radius of gyration of 
 the interface c atoms ( rg ) . the linear combination 
 of pcs that best correlates with the rg of each interface 
 ( maximally correlated motion or mcm ) was then calculated .
the contribution 
 of single pcs to the mcm , together with the value of the mcm - rg correlation , determines the fraction of the rg variance that 
 is due to the motion along a given pc ( pvar ) .
a pc was considered 
 to be correlated with an interface if pvar  20% .
the 
 number of rg - correlated pcs ( npc20 ) was calculated for each interface 
 ( figure 5a ) .
all the interfaces turned out 
 to be correlated on average with no more than 1 pc .
the behavior of cmono interfaces ( cyan ) was similar to that of 
 the multipartner ones ( orange ) , while the cmono_in_multi interfaces ( blue ) seemed to be slightly less affected by the protein 
 global motions .
if the number of unique rg - correlated pcs is 
 summed over all the interfaces of a given protein ( npc20 in figure 5b ) , a significant difference is found between monopartner 
 ( cyan ) and multipartner ( magenta ) proteins , indicating that the latter 
 use a larger number of independent global motions to modulate the 
 shape of their interfaces . since for multipartner proteins
each interface 
 is on average correlated with 0.9 pc and each protein has , 
 on average , four nonredundant interfaces , the peak at 3 for multipartner npc20 implies that 
 different interfaces tend to be correlated with different pcs .
qualitatively 
 similar results were obtained when performing the fma analysis on 
 md ensembles of the s data set ( si 
 figure s7 ) .
analysis of correlated motions in s. ( a / b ) 
 correlation 
 between global and interface motions by functional mode analysis ( fma ) 
 on tconcoord ensembles from the s data set .
( a ) number 
 of pcs accounting for at least 20% of the rg variance of a given 
 interface ( npc20 ) .
average values calculated over interfaces of the three binding 
 classes cmono ( cyan ) , cmono_in_multi ( blue ) , and cmulti ( orange ) are reported .
( b ) distributions of npc20 per protein ( npc20 ) , calculated as the number of unique pcs accounting 
 for at least 20% of the rg variance of any of the protein interfaces . 

values for monopartner ( cyan ) and multipartner ( magenta ) proteins 
 are reported .
( c ) distributions of normalized mi ( i ) calculated over pairs of surface residues 
 in monopartner ( light gray , smono ) and 
 multipartner ( dark gray , smulti ) proteins , 
 and over pairs of interface cmono ( cyan ) , cmono_in_multi ( blue ) , and cmulti ( orange ) residues . in a and c
the significance levels 
 from pairwise wilcoxon comparison tests are reported with a star code 
 ( see figure 2 legend ) .
the pcs used in the fma represent the main global 
 or 
 collective motions of a protein .
correlations between residues can 
 also be analyzed by considering their local dynamics , that is , by 
 removing the overall roto - translation of the protein from their motion . 

indeed , it has been shown that analyzing local dynamics can highlight 
 communication pathways within the protein that are difficult to identify 
 solely from the collective motions involving the entire structure . to this end
, local motions were analyzed with 
 a fragment - based approach and the extent of 
 correlation was estimated by normalized mutual information i ( eq 1 in methods ) . correlations within interface 
 residues of each binding class ( cyan , blue and orange in figure 5c ) calculated from the tconcoord ensembles were 
 found to be significantly higher than those within mono- and multipartner 
 surface residues ( light and dark gray ) , suggesting a higher level 
 of communication between interface residues .
interestingly , the cmulti distribution ( orange ) was found to be 
 significantly higher than the cmono_in_multi one ( blue ) , indicating that , within multipartner proteins , multipartner 
 residues are on average more correlated than monopartner ones .
no 
 significant differences were instead found between cmulti ( orange ) and cmono ( cyan ) 
 distributions . in this section
, we will investigate 
 if the higher intrinsic flexibility found for multipartner residues 
 correlates with the variability observed in the experimental structures . 

indeed , many proteins in the s data set have a relatively 
 large number of occurrences in the pdb ( > 10 for 50% of 
 the 
 proteins ) .
these pdb ensembles contain information on the protein 
 conformational variability that can be 
 extracted and compared with that derived from the simulated ensembles .
the overall variability within the pdb ensemble of each protein 
 was decomposed into three different contributions according to the 
 binding state of the structures that are compared : unbound 
unbound 
 ( u  u ) , bound  bound ( b  b ) , and unbound  bound 
 ( u 
b ) . the highest structural variability ( as measured by the 
 maximum c rmsd calculated over all the structure 
 pairs of a given protein ) was found on average for the u 
b 
 pairs ( 2.66  0.36  ) , followed by u  u ( 1.81  
 0.28  ) and b  b ( 1.34  0.15  ) . the u  u 
 rmsd values can be considered as related to the intrinsic plasticity 
 of the isolated protein , which can have different accessible states 
 (  pre - existing equilibrium ) , while the u 
b values 
 include also the structural changes caused by the binding of the partners 
 (  induced fit ) .
the smaller variability observed within 
 bound structures ( b  b ) reflects the higher number of constraints 
 in the complexes .
for each of these different contributions , 
 multipartner proteins 
 showed on average a higher conformational variability ( u  umulti , b 
bmulti , and u  bmulti , dark colors in figure 6a ) than monopartner 
 ones ( u  umono , b  bmono , and u  bmono , light colors ) .
this suggests that the higher intrinsic 
 flexibility observed for multipartner residues in isolated proteins 
 could be used to enhance pre - existing equilibrium or to assist induced - fit 
 changes .
indeed , the single contributions of the residues to the overall 
 c fluctuation within the pdb ensemble ( figure 6b  c ) showed a similar picture as that observed 
 for the simulated rmsf ( see also si figure s4 for a direct comparison ) . a higher conformational variability ( figure 6b , light colors ) and a weaker dependence from conservation 
 ( figure 6c )
were found for multipartner residues , 
 resulting in an increased rmsf difference between highly conserved 
 mono and multipartner residues ( figure 6b , 
 dark colors ) .
the cmono and cmono_in_multi groups were instead more similar to each 
 other than in the simulated case , indicating a comparable degree of 
 variation within the pdb structures .
consistently with the simulated 
 ensembles , reduced differences were found between the different binding 
 classes when the side chain variability was considered ( si figure s8 and s4 ) .
analysis of the conformational 
 variability within pdb ensembles 
 in s. ( a ) averages of maximum c rmsd 
 values calculated over all the possible pairs of pdb structures within 
 the unbound ( pink ) and bound ( green ) ensembles , and between bound 
 and unbound structures ( magenta ) of each mono- ( light color ) and multi- 
 ( dark color ) partner protein .
the significance levels from pairwise wilcoxon 
 comparison tests between ( mono , multi ) pairs of distributions are 
 reported with a star code ( see figure 2 legend ) . 

( b ) average c rmsf z - scores ( pdb ensembles ) calculated 
 over cmono ( cyan ) , cmono_in_multi ( blue ) , and cmulti ( orange ) residues in the s data set .
the significance 
 levels from pairwise wilcoxon comparison tests are reported with a 
 star code ( see figure 2 legend ) .
averages and 
 significance levels calculated considering only the residues with 
 the highest consurf conservation grade ( gcons = 9 ) are also reported in dark colors .
( c ) dependence of pdb average 
 c rmsf z - scores ( dots ) from evolutionary conservation 
 for cmono ( cyan ) , cmono_in_multi ( blue ) , cmulti ( orange ) 
 residues .
residues are partitioned into 9 groups according to their 
 consurf conservation grade ( gcons ) .
a 
 best - fit linear regression is also reported for each binding class . in this section ,
we provide further 
 insight on the functional relevance 
 of promiscuous residues by analyzing the distribution of single nucleotide 
 polymorphisms ( snps ) across the different binding classes considered 
 in the previous sections .
recent large - scale studies of the 
 human genome such as the international hapmap project and the 1000 genomes project have produced a large number of data that can be used to 
 accurately assess the relationship between genotype and phenotype . 

in particular , snps are single nucleotide variations observed at a 
 specific location of the genome in at least 1% of the population .
the mapping of snps , and in particular nonsynonymous 
 snps , to protein regions is currently exploited in a wide range of 
 applications , from disease association studies to pharmacogenomics . in this study , only missense nonsynonymous 
 snps are considered . in order to study the relationship between 
 binding promiscuity 
 and human snps , human homologues of s and s proteins were identified .
we recorded the occurrence of nonsynonymous 
 snps from the dbsnp database and of 
 snps with a known association with disease ( dissnps ) from the omim 
 database ( methods ) .
the snp and dissnp positions were then mapped back to the original 
 s and s proteins .
the comparison 
 of the snp propensities of the different binding 
 classes with respect to the surface ( figure 7a , light colors ) showed that promiscuous positions in s ( orange ) are less rich in snps than both classes of monopartner 
 residues ( cyan and blue ) .
even if with reduced statistical significance , 
 this trend was confirmed by s proteins ( green and red ) . 
 as an example , the survivin protein , already in the human form in 
 the s data set , is shown in figure 7b . while promiscuous ( red shades surface ) and nonpromiscuous 
 ( light blue surface ) residues are present in survivin in almost equal 
 proportion , 5 out of 6 interface snps ( spheres )
were found in nonpromiscuous 
 regions , mainly involved in the interaction with the survivin partner 
 borealin ( green cartoon ) .
( a ) propensities of snps ( light colors ) 
 and dissnps ( dark 
 colors ) relative to the surface for cmono ( cyan ) , cmono_in_multi ( blue ) , cmulti ( orange ) , cmono_in_multi(soc ) ( green ) and cmulti(soc ) ( red ) interface 
 residues .
the reported values 
 are averages over s monopartner proteins for cmono , s multipartner proteins 
 for cmono_in_multi and cmulti , and s proteins for cmono_in_multi(soc ) and cmulti(soc ) .
the significance 
 levels from pairwise student s t tests are 
 reported with a star code ( see figure 2 legend ) . 

snps found in the interface 
 region of survivin are labeled and represented as van der waals spheres . 

the protein surface ( transparent ) is colored according to the binding 
 multiplicity b from blue ( b = 0 , 
 non - interface ) to red ( b = 3 ) .
a survivin binding 
 partner ( borealin ) is also represented as green cartoon ( pdb i d : 2raw ) .
a possible explanation of these findings is that 
 the human equivalent 
 of the promiscuous positions considered here tend to be less tolerant 
 to variation , probably due to the higher number of constraints that 
 they are experiencing to preserve effective binding .
mutations at 
 these sites might be more prone to yield a lethal phenotype and are 
 thus less likely to be viable .
the analysis of dissnps ( figure 7a , dark 
 colors ) was strongly affected by the small number of observations 
 ( si table s12 ) .
the large uncertainty associated 
 with the average propensities , especially for the s proteins , 
 prevented the observation of statistically significant differences 
 between the binding classes .
it is thus evident that these results , 
 while indicating an interesting and unexplored relationship between 
 binding diversity and variability in the human genome , will need to 
 be confirmed on a larger data set of human proteins . in this section
, 
 we will exemplify the relationship observed between binding and flexibility 
 analyzing two related multipartner proteins : neddylin and the small 
 ubiquitin - related modifier 2 ( sumo-2 ) .
neddylin is a ubiquitin 
 homologue ( 56% sequence identity , figure 8a ) 
 with the characteristic ubiquitin - like fold ( figure 8b / d ) .
the analysis of its complexes highlighted that , as for 
 ubiquitin , the neddylin main interface 
 is centered at an hydrophobic patch at the c - terminus of the 5 
 strand ( figure 8b , green circle ) , involving 
 different hydrophobic residues from the 1-2 loop ( l8 ) 
 and the 35 strands ( i44 , v70 , l71 ) .
this region is 
 highly promiscuous and rich of hot spots ( sticks in figure 8b ) .
a secondary interaction site , previously observed 
 also in ubiquitin , was found at the 
 1 c - terminus ( figure 8b , yellow circle ) . 

this does not show any hot spot and is mainly composed of monopartner 
 residues .
( a ) sequence alignment 
 of ubiquitin ( uniprot / kb i d : p0cg48 ) , neddylin ( q15843 ) , and sumo-2 
 ( p61956 ) .
the dssp secondary 
 structure of ubiquitin ( pdb i d : 1ubi ) is shown ( top ) . conserved residues are 
 highlighted in red .
stars below neddylin and sumo-2 sequences indicate cmono_in_multi ( blue ) and cmulti ( orange ) residues identified as hotspots by at least 
 two prediction methods .
( b / c ) cartoon representation of neddylin ( b , 
 pdb i d : 1ndd ) and sumo-2 ( c , pdb i d : 1wm3 ) structures .
residues are colored according to the 
 binding multiplicity b from blue ( b = 0 , non - interface ) to red ( b = 4 ) .
residues identified 
 as hotspots by at least two prediction methods are represented as 
 sticks .
the approximate boundaries of two representative interfaces 
 of neddylin are indicated in b with a green ( main interface ) and yellow 
 ( secondary interaction site ) dashed line .
the insets show the structures 
 of the representative complexes ( pdb i d 1xt9 , chains b / a , and pdb i d 1r4n , chains j / c in the 
 green and yellow insets , respectively ) .
( d / e ) mapping of tconcoord 
 c rmsf values onto neddylin ( d ) and sumo-2 ( e ) structures . 

the structures are colored according to rmsf values from blue ( 0  ) 
 to red ( 2  ) .
selected residues ( highlighted with yellow arrows 
 in panel a ) are represented as yellow sticks , while hydrogen bonds 
 in d are indicated as dotted blue lines .
( f / g ) profiles of binding 
 multiplicity ( b , blue ) and c rmsf 
 from tconcoord ( orange ) for neddylin ( f ) and sumo-2 ( g ) .
dssp annotation 
 of the secondary structure is reported as magenta ( helices ) and yellow 
 ( strands ) blocks .
a positive correlation 
 was observed between the profiles of binding 
 multiplicity ( orange ) and c rmsf from the tconcoord 
 ensemble ( blue ) of neddylin ( figure 8f ) , with 
 a correlation coefficient calculated over the exposed residues of 
 0.36 .
correspondingly , a higher average c rmsf z - score 
 ( si table s13 ) was found for multipartner 
 residues ( 0.79  1.48 ) than monopartner ones ( 0.03  0.72 ) . 

indeed , many multipartner residues ( figure 8b ) are located at or in close proximity of high - flexibility regions 
 ( figure 8d ) , namely the 1-2 , 
 1-3 , and 3-4 loops , and the c - terminus . 

the intrinsic backbone mobility of promiscuous locations in the simulated 
 ensemble is paralleled by a high conformational variability at these 
 same positions in the different pdb complexes where neddylin interacts 
 with different partners ( si table s13 ) . 

additionally , the fma analysis ( figure 9a ) 
 showed that the backbone flexibility of each of the two representative 
 interfaces mainly correlates with one specific collective motion . 
 in particular , pc4 ( figure 9b ) accounts for 
 54% of the rg variance of the main interface ( green bars in figure 9a ) . a similar collective motion , involving a  pincer - like 
 movement of the 1-2 and 1-3 loops , has 
 been observed in nmr ensembles of ubiquitin representing solution 
 dynamics up to the s time scale . fma 
 analysis of neddylin . ( a ) decomposition of the rg variance 
 of the two representative interfaces of neddylin into percentage contributions 
 ( pvar ) from the first 17 pcs of the tconcoord ensemble .
green and 
 yellow bars represent the pvar values for the interfaces shown in 
 the green and yellow insets in figure 8b .
( b ) 
 porcupine representation of the 4th pc of the neddylin tconcoord ensemble . 

direction and relative amplitude of the motion of each c atom along the pc is represented by orange spikes .
the residues 
 involved in the main neddylin interface ( figure 8b , green inset ) are color - mapped according to their binding multiplicity b onto the tube representation of the average tconcoord 
 structure .
multipartner residues of neddylin 
 seemed to rely less on side chain 
 flexibility than backbone mobility to adapt to different environments , 
 since the average side chain rmsf was comparable to that of monopartner 
 residues for both the tconcoord and pdb ensembles ( si table s13 ) . indeed , while a few multipartner residues adopted 
 different rotamers in different complexes ( e.g. , r42 in si figure s9b ) , others relied either on the backbone 
 flexibility ( l8 in si figure s9a ) or on 
 their capacity to support different interaction geometries without 
 changing their conformation ( e.g. , h68 in si figure 
 s9c , where it interacts with an aromatic ring in either a t - shaped 
 or parallel stacking interaction ) . the sumo-2 protein , which 
 is a homologue of ubiquitin in spite 
 of its low sequence identity ( 14% , figure 8a ) , showed important differences with respect to neddylin in binding 
 modes and dynamics .
this is consistent with this protein belonging 
 to a separate sequence subgroup of ubiquitin 
 homologues , characterized by the replacement of the key ubiquitin 
 residues q41 and y59 ( yellow arrows and sticks in figure 8a and d , respectively ) with two hydrophobic residues 
 that are no longer able to form hydrogen - bonds with the nearby loops 
 ( sticks in figure 8e ) .
this is reflected in 
 both the protein structure and dynamics . in particular , compared to 
 neddylin the tconcoord sumo-2 ensemble presented a higher c flexibility at the 1 c - terminus and the 1/3 
 loop ( figure 8e ) .
the increment in flexibility 
 was paralleled by an increment in 
 binding promiscuity in this region .
a multipartner interaction site , 
 typical of sumo proteins , was found at the  groove 
 defined by the 1 c - terminus and the 2 strand , which 
 is correspondingly enriched in multipartner hot spots with respect 
 to neddylin ( figure 8c and g ) .
the mapping 
 of tconcoord c and side chain rmsf onto the sumo-2 
 surface ( figure 10c / d ) , shows that the / 
 groove is flanked by residues with high flexibility either at c ( 1-2 loop ) or side chain atoms ( 2 ) 
 or both ( 1 c - terminus ) .
moreover , the water density map from 
 the md simulation showed that the region surrounding the groove is 
 richer in high - density hydration sites than the nonpromiscuous regions 
 ( si figure s10 ) .
( a / b ) cartoon representation of two different complexes 
 of sumo-2 ( light green ) .
selected multipartner hot spots of sumo-2 
 ( blue ) and interacting residues from partners ( yellow ) are shown as 
 sticks .
( b ) complex with thymine - dna glycosylase ( orange , pdb 
 i d : 2d07 ) . ( c / d ) 

c ( c ) and side chain ( d ) rmsf values mapped onto 
 sumo-2 surface ( pdb i d : 1wm3 ) .
residues are colored according to their rmsf values 
 from blue ( 0  ) to red ( 2  ) .
the analysis of the sumo-2 pdb structures ( si table s14 ) highlighted that , differently from neddylin , 
 side chain conformational changes are more important for multipartner 
 residues ( average side chain rmsf z - score = 0.56  1.13 ) than 
 for monopartner ones ( 0.06  0.98 ) .
particularly relevant 
 is the contribution of the r34 and q35 hot spots at the 1 c - terminus , 
 with a side chain rmsf z - score of 2.63 and 1.32 , respectively .
correspondingly , 
 these residues adopt significantly different rotamers in the different 
 pdb complexes ( figure 10a / b ) .
this example 
 suggests that a fine - tuning of the intrinsic dynamical 
 properties of the interface can be central in modulating the binding 
 specificity in evolutionary related proteins .
in this work , we studied the intrinsic conformational flexibility 
 of multipartner proteins to assess its role in promoting diversity 
 of binding . through the generation of simulated conformational ensembles 
 from a starting experimental structure , we measured the tendency of 
 the isolated proteins to sample different conformations independently 
 from the interactions with their partners .
the conformational flexibility 
 that we considered here is a different concept from the intrinsic 
 disorder analyzed in other works on hub proteins , characterized by the absence of a definite structure 
 in all or part of the protein .
it is also different from the conformational 
 plasticity as measured by the conformational changes observed 
 when multiple experimental structures are available .
indeed , if these 
 structures correspond to different states , such as bound and unbound 
 conformations , the conformational variability is likely to include 
 also induced fit and allosteric effects in addition to intrinsic 
 flexibility .
the choice to consider the isolated protein is justified by many 
 studies indicating that the intrinsic dynamics of a protein is correlated 
 with its function .
indeed , even when isolated , 
 proteins have been shown to sample functionally relevant states , which 
 can then be stabilized or selected by interactions with the environment 
 such as post - translational modifications or interactions with ligands .
many of these works used elastic network models 
 ( enm ) to characterize the equilibrium fluctuations within the native 
 structure basin . in spite of their simplicity
, 
 these methods have been shown to provide results in good agreement 
 with both experiments and more sophisticated approaches such as md . 

for this study , we chose a method of intermediate complexity and computational 
 cost , tconcoord .
while allowing for anharmonicities and providing 
 a full atom description of the protein where both backbone and side 
 chains are included , it is still faster than md simulations .
the interaction data on each protein of 
 our s data 
 set were derived from a structural ppi database , pisite , where partners from all the pdb complexes involving 
 homologues within a high ( > 90% ) sequence identity family are mapped 
 onto a family representative .
structural ppi databases generally contain 
 higher - confidence interaction data than interactome networks derived 
 from a range of different experimental methods , resulting in degree 
 distributions with shorter tails . in 
 our s data set , this is reduced also by the fact that 
 only monodomain proteins were considered , so that a maximum number 
 of 12 nonredundant partners per protein was found .
the investigation 
 was limited to monodomain proteins to exclude cases where multiple - partner 
 binding is simply achieved by using different domains or a different 
 relative arrangement of the domains .
the relationship between conformational flexibility and ability 
 to bind multiple partners was analyzed primarily at the residue level . 

indeed , the composition of either proteins or interfaces in terms 
 of binding multiplicity turned out to be highly heterogeneous , with 
 40% of multipartner proteins in s presenting both mono- 
 and multipartner interfaces and 91% of multipartner interfaces including 
 both mono- and multipartner residues .
thus , we aimed mostly at identifying 
 the dynamical features specific of multipartner residues and not at 
 characterizing multipartner proteins as a whole . even if we 
 did not attempt a classification of the s multipartner 
 proteins , more than 60% of them have at least 10% of 
 their interface residues involved in interactions with different partners , 
 while only 8% multipartner proteins have no overlapping interfaces . 
 on the other hand , as found in different data sets , very large overlaps between interfaces
seem to be avoided 
 ( only 5% of the multipartner proteins have more than 60% of their 
 interface interacting with multiple partners ) .
the alternate data 
 set s used to validate our results , is composed of proteins 
 previously classified as  sociable or 
transient 
 hubs , that can be related , even if with 
 a different definition , to date hubs .
we found a significant 
 difference in the intrinsic conformational 
 flexibility of mono- ( cmono and cmono_in_multi ) and multi- ( cmulti ) partner residues .
in particular , the comparison 
 between cmono_in_multi and cmulti residues suggests a nonuniform distribution of the 
 flexibility across the interface of multipartner proteins , with multipartner 
 residues being on average more flexible than monopartner ones .
this 
 holds when considering either backbone or side chain motions , suggesting 
 that an enhanced ability to sample different conformations , either 
 globally or locally , can indeed be exploited to support diversity 
 of interactions .
remarkably , a more pronounced tendency of promiscuous 
 residues to sample different conformations is observed in all the 
 different types of conformational ensembles analyzed in this study 
 and in two different protein data sets .
the flexibility 
 of multipartner proteins was found to modulate the shape of single 
 interfaces in a highly specific way , by using different collective 
 motions for different interfaces within the same protein .
moreover , 
 within multipartner proteins a stronger correlation was found between 
 the local motions of multipartner residues , suggesting that promiscuous 
 regions are connected by preferential communication pathways .
similarly to disorder , a higher flexibility 
 at the interface of isolated molecules is associated with a higher 
 entropic penalty upon binding .
this entropy increment has been suggested 
 to be exploited by disordered multipartner proteins to decouple binding 
 affinity from binding specificity .
on 
 the other side , residues considered to be important in the interaction 
 with the partner , either because of their high evolutionary conservation 
 or for being predicted as hot spots , have been shown to be more rigid 
 than the average . here , we found 
 that ( a ) there is a clear anticorrelation between evolutionary conservation 
 and intrinsic flexibility and ( b ) this trend is reduced in the case 
 of multipartner residues , suggesting that the higher flexibility of 
 conserved multipartner residues is the result of a balance between 
 the counteracting effects of preserving binding strength and allowing 
 for binding diversity . in this context , the higher propensity for 
 conservation found for multipartner residues in ordered proteins could 
 reflect the higher number of evolutionary constraints deriving from 
 the necessity to optimize this balance .
the findings on intrinsic 
 flexibility are confirmed by the analysis 
 on the conformational variability observed in the experimental structures . 

previous works using pdb ensembles of different data sets have shown 
 that overlapping regions in date hub interfaces tend to visit more 
 different side chain rotamers than nonoverlapping ones and that sociable / hub proteins sample more different 
 overall backbone conformations than nonsociable / nonhub ones . here
, we unified these results by performing 
 a systematic analysis on both side chain and backbone changes at the 
 residue level .
moreover , we decomposed the overall variability observed 
 in the pdb into contributions from pairs of structures with the same 
 and with different binding states .
interestingly , we found a higher 
 plasticity for multipartner proteins not only when considering changes 
 between unbound and bound structures but also when analyzing ensembles 
 of unbound and bound structures separately .
this is highly consistent 
 with the larger intrinsic flexibility found for the multipartner residues in isolated proteins , since it implies a more pronounced tendency to explore different 
 conformations independently from their binding state .
the higher 
 propensity to conservation and surface burial upon complexation 
 of multipartner residues suggests that they have an important role 
 in defining the binding energetics .
the comparison of the intrinsic 
 flexibility of predicted hot spot residues from different binding 
 classes confirmed the results obtained on the whole set of interface 
 residues .
indeed , multipartner hot spots turned out to be on average 
 more flexible than monopartner ones independently from the specific 
 method used for hot spot prediction .
this rules out the possibility 
 that the observed higher flexibility of multipartner residues in the 
 whole interface originates from residues that are only marginally 
 involved in the interaction with partners .
moreover , the observation 
 that promiscuous regions are depleted in snps , even if to be confirmed 
 on larger data sets of human proteins , provides further evidence to 
 the essentiality of these regions .
the findings presented in 
 this work have potential applications 
 to methods for the prediction of ppis , whose accuracy has been recently brought to levels comparable to 
 high - throughput experiments . in particular
, 
 the introduction of per - residue flexibility indices could be used 
 for the identification of promiscuous regions in protein interfaces . 

this could also be relevant for the detection of druggable regions , 
 which could be targeted by small molecules or other proteins to inhibit 
 interactions with specific partners .
a further possible application is in protein design , where the specificity 
 or promiscuity of proteins could be enhanced by modifying the distribution 
 of flexibility across the interfaces .
 the 
 ability to interact with different partners is one of the most 
 important features in proteins . 
 proteins that bind a large number 
 of partners ( hubs ) have been often associated with intrinsic disorder 
 . 
 however , many examples exist of hubs with an ordered structure , and 
 evidence of a general mechanism promoting promiscuity in ordered proteins 
 is still elusive . 
 an intriguing hypothesis is that promiscuous binding 
 sites have specific dynamical properties , distinct from the rest of 
 the interface and pre - existing in the protein isolated state . here 
 , 
 we present the first comprehensive study of the intrinsic dynamics 
 of promiscuous residues in a large protein data set . 
 different computational 
 methods , from coarse - grained elastic models to geometry - based sampling 
 methods and to full - atom molecular dynamics simulations , were used 
 to generate conformational ensembles for the isolated proteins . the 
 flexibility and dynamic correlations of interface residues with a 
 different degree of binding promiscuity were calculated and compared 
 considering side chain and backbone motions , the latter both on a 
 local and on a global scale . 
 the study revealed that ( a ) promiscuous 
 residues tend to be more flexible than nonpromiscuous ones , ( b ) this 
 additional flexibility has a higher degree of organization , and ( c ) 
 evolutionary conservation and binding promiscuity have opposite effects 
 on intrinsic dynamics . 
 findings on simulated ensembles were also validated 
 on ensembles of experimental structures extracted from the protein 
 data bank ( pdb ) . additionally , the low occurrence of single nucleotide 
 polymorphisms observed for promiscuous residues indicated a tendency 
 to preserve binding diversity at these positions . 
 a case study on 
 two ubiquitin - like proteins exemplifies how binding promiscuity in 
 evolutionary related proteins can be modulated by the fine - tuning 
 of the interface dynamics . 
 the interplay between promiscuity and flexibility 
 highlighted here can inspire new directions in protein  
 protein 
 interaction prediction and design methods . 
Introduction
Methods
Results
Discussion