the specificity of interactions is mostly determined by structural 
 and physico - chemical properties of interfaces of interacting proteins .
proteins highly similar in sequence can participate 
 in different interactions , and even a small number of amino acid substitutions 
 or insertions / deletions at binding interfaces can lead to different 
 interaction partners or different binding modes . introducing mutations at protein  protein interfaces 

can quite effectively modulate binding affinity or specificity , and 
 site directed mutagenesis has been successfully used in rational protein 
 design and directed evolution to create novel protein complexes . in some cases , the redesign of protein  protein interactions 
 may yield new binding partners with differences in specificities of 
 at least 300-fold between the cognate and the noncognate complexes .
protein interactions and lead to mendelian diseases 
 or cancer . one possible way to assess the effect of mutation on protein 
 binding 
 affinity is to experimentally measure it .
however , although site - directed 
 mutagenesis methods are inexpensive and fast , surface plasmon resonance , 
 isothermal titration calorimetry , fret and other methods used to measure 
 binding affinity can be time - consuming and costly .
therefore , the 
 development of reliable computational approaches to predict changes 
 in binding affinity upon mutation is urgently required .
it is especially 
 true in light of the rapid development of next generation sequencing 
 technologies , which are producing overwhelming amounts of data on 
 polymorphisms and disease mutations .
there have been many studies 
 aiming to develop potentials for protein ligand binding ( see references within ) . according to some studies ,
the burial of 
 polar groups should not be penalized in the case of hydrophilic protein 
, the contribution of electrostatic interactions 
 might be more pronounced when switching from protein folding to protein  protein 
 interaction potentials .
many models have been proposed that try to estimate binding 
 affinities of protein  protein complexes , yet they can not reproduce 
 experimental binding data with satisfactory accuracy .
this happens 
 mostly due to the difficulties associated with modeling of conformational 
 changes upon binding , estimating the entropy changes and long - range 
 noninterfacial interactions .
moreover , 
 most of the approaches are largely tuned to predict the effect of 
 alanine - scanning mutations defined as substitutions of residues into 
 alanine .
energy functions used to describe and predict interactions 
 governing 
 protein  protein binding linearly combine several energy terms 
 and the coefficients in front of each term can be optimized using 
 datasets of experimentally measured binding affinities or changes 
 in binding affinity due to missense mutations .
several approaches 
 have recently been proposed that predict the effects of point mutations 
 on binding energy .
they include coarse - grained predictors based on 
 statistical or empirical potentials , molecular mechanics 
 force fields with different solvation models and others . 

for example , the molecular mechanics poisson  boltzmann surface 
 area ( mm - pbsa ) method has been shown to yield good agreement with 
 experimental studies in determining protein stability , binding affinity 
 and ranking docking templates .
certain methods specifically 
 investigate the impact of mutations on binding in relation to function 
 or different environmental conditions .
for example , recently an in 
 silico mutagenesis method was proposed to calculate the effect of 
 mutation at various ph , whereas another 
 approach predicted the mutations ability to shift the allosteric 
 conformational equilibrium with consequences for functional regulation .
although many of these methods report reasonable 
 root - mean - square error ( rmse ) values , their applicability is biased 
 toward very limited sets of proteins and mutations ( on an order of 
 a dozen protein complexes and up to several hundreds of mutations ) 
 used for training and parametrization . despite their limitations , 

methods that estimate changes in affinity produced by mutations show 
 somewhat better performance compared to computational approaches to 
 predict absolute affinities . it can be partially attributed to counterbalancing 
 of errors , produced by energy calculations or conformational sampling , 
 between the wild - type and mutant structures , under an assumption that 
 structures do not undergo large changes upon mutations . here
we introduce a new approach to assess the impact of single 
 and multiple missense mutations on protein binding affinity .
it uses 
 a modified mm - pbsa energy function , statistical scoring function and 
 efficient energy minimization protocol .
our protocol is parametrized 
 on a large set of several thousands of mutations and their corresponding 
 experimental binding free energy changes assembled from the scientific 
 literature .
we find that the choice of 
 water model is very important for the quality of prediction and a 
 simulation protocol with explicit water without restraints on the 
 backbone atoms gives the best results .
we show that the conformational 
 sampling by molecular dynamics simulations does not help to achieve 
 better quality predictions for nm set ( see methods ) , and further analyses are needed to prove or disprove this observation 
 for larger sets of mutations .
in addition , we describe how prediction 
 accuracy depends on the type of amino acid substitution and its location 
 in the complex .
compared to several available computational protocols , 
 our approach achieves very good accuracy and speed .
it is benchmarked 
 against a large and independent experimental dataset and shows high 
 correlation coefficients between predicted and experimental values 
 of 0.69 and 0.63 with root - mean - square errors of 1.20 and 1.90 kcal 
 mol for single and multiple mutations , respectively . 
 at the same time
, it can efficiently handle the large amount of data 
 generated by modern genomics techniques .
we use two datasets of experimentally determined 
 values of changes in standard binding free energy upon mutations ( denoted 
 as g in the paper ) , which are calculated 
 as the difference between binding affinities of mutant ( gmut ) and wild - type ( gwt ) complexes and are derived 
 from the scientific literature for protein  protein complexes 
 with experimentally determined wild - type structures .
one set is taken 
 from the paper of bockmann et al . and 
 includes 242 single mutations and 123 multiple mutations from nine 
 wild - type protein  protein complexes (
we retained only structures of wild - type proteins , eliminated redundant 
 entries and made sure that the skempi set did not overlap with the 
 nm set .
we also excluded ten single and eight multiple mutants having 
 initial vmd models ( see next section ) 
 with unrealistically large steric clashes between the substituted 
 and adjacent residues that could not be fixed by the minimization 
 procedure ( table s1 , supporting information ) . as a result
, our second experimental set included 1844 single 
 and 574 multiple mutations from 81 wild - type protein  protein 
 complexes ( it will be referred to as the  skempi set  

results for our and other methods and main scripts for 
 running the programs are accessible through ftp://ftp.ncbi.nih.gov/pub/panch/mutation_binding
. an energy minimization procedure 
 is used to compute the equilibrium configuration and find a minimum 
 of the potential energy of the system . here
we use a conjugate gradient 
 algorithm implemented in namd , which 
 finds a nearest minimum to the starting point , repairs distorted geometries 
 and removes steric clashes .
the initial crystal structures of wild - type 
 protein  protein complexes were obtained from the protein data 
 bank ( pdb ) ; only protein chains were 
 retained in the calculations and missing heavy side chain atoms and 
 hydrogen atoms were added using the vmd ( version 1.9.1 ) program . the procedure for adding atoms in vmd
we tested several 
 simulation protocols ( figure s1 , supporting information ) ( 1 ) with and without restrains on the backbone atoms upon minimization 
 using explicit tip3p water model ( 2 ) with an unconstrained molecular 
 dynamics ( md ) simulations with tip3p explicit water model and ( 3 ) with the minimization using the generalized 
 born implicit solvent model for both 
 wild - type and mutant structures .
a simulation protocol that used 
 minimization without restraints on the backbone atoms in explicit 
 water showed the best agreement with the experiments ( see results ) .
all models were immersed into rectangular boxes of water 
 molecules extending up to 10  from the protein in each direction . 
 to ensure an ionic concentration of 150 mm and zero net charge ,
the distribution 
 of the system size including water molecules and ions is shown in 
 figure s2 ( supporting information ) .
all 
 wild - type protein complexes were minimized using a 40  000-step 
 energy minimization procedure to make all wild - type protein complexes 
 fully optimized for further binding energy calculation and the preparation 
 of the initial mutant structures .
this 40  000-step energy minimization 
 procedure included an initial 5000-step minimization that was carried 
 out using harmonic restraints ( with the force constant of 5 kcal mol  ) applied on the backbone 
 atoms of all residues , followed by a 35  000-step energy minimization 
 on the whole system .
we used the final minimized models of wild - type 
 protein complexes 
 to produce the mutant structures ( figure s1 , supporting 
 information ) ; namely , we introduced an amino acid substitution 
 using the  mutator  plugin of the vmd and kept the rest 
 of the system invariant including the water box and ions .
then we 
 performed an additional 1000-step minimization of the mutant ( altogether 
 40  000 steps for wild - type and 1000 steps for mutant minimization ) 
 and showed that the agreement with the experimental data increases 
 as the number of steps increases but reaches a plateau after 300500 
 steps for single mutations .
figure s3 ( supporting 
 information ) shows that increasing the number of minimization 
 steps for a mutant above 300500 does not improve the agreement 
 with the experiments .
moreover , we showed that even very large complexes 
 did not require more than 300500 minimization steps for mutants 
 ( figure s4 , supporting information ) .
the 
 minimized models of mutant protein complexes after 300 minimization 
 steps were used in our further analysis .
minimization was done only 
 for the protein complexes , and protein structures of monomers were 
 extracted from the minimized complexes for the following energy calculation . 

consistent with the previous studies , this approach was more accurate 
 and faster than when minimization was done for complexes and monomers 
 separately .
the energy minimization and md simulation were carried 
 out with 
 the namd program version 2.9 using the 
 charmm27 force field .
periodic boundary 
 conditions and a 12  cutoff distance for nonbonded interactions 
 were applied to the systems .
the particle mesh ewald ( pme ) method was used to calculate the long - range electrostatic 
 interactions .
lengths of hydrogen - containing bonds were constrained 
 by the shake algorithm . unlike minimization , 

md simulation generates atomic trajectories and may be used to generate 
 conformational ensemble and determine macroscopic properties of the 
 system by ensemble averages .
in addition , md simulations can in principle 
 explore conformations which are not accessible via energy minimization 
 protocol because of the energy barriers between different conformational 
 states .
binding energies were calculated 
 based on the mm - pbsa method that combines the molecular mechanical 
 terms with the poisson  boltzmann continuum representation of 
 the solvent calculated using the charmm 
 force field .
the total free energy in mm - pbsa is expressed as the 
 following:1here egasmm corresponds to the molecular 
 mechanical energy and is calculated as a sum of the internal energy 
 of the molecule , the van der waals and coulomb electrostatic interactions 
 in the gas - phase .
gsolvp is the polar contribution to the solvation 
 free energy of the molecule calculated using the poisson 
boltzmann 
 ( pb ) equation , which is the difference between the electrostatic energy 
 of solute in the solvent environment and in the reference environment 
 ( gas phase in our study).gsolvnp is the nonpolar 
 solvation energy , which is proportional to the solvent - accessible 
 surface area ( sa ) of a molecule and ts 
 corresponds to the entropy of solute
. as entropy calculations are 
 very computationally expensive and might not be very accurate , we did not perform entropy calculations in our 
 analysis . moreover , as was discussed in previous studies , entropy 
 terms may cancel each other later in eq 3 , especially 
 for cases when structures do not undergo large changes upon mutations 
 in the minimization simulation protocol .
the binding energy in mm - pbsa approach is usually 
 calculated as a difference between the average energies of the complex 
 and each monomer:2then the change of the binding energy due 
 to a mutation can be calculated as the following:3 we performed a 
 multiple regression 
 fitting procedure to model the linear relationship between the response 
 ( dependent ) variable , which is , in our case , experimental values of 
 changes in binding affinity , and different dependent variables .
we 
 tried different types of energy functions composed of different energy 
 components ( see results section ) and estimated 
 the optimal model parameters ( the weight factors for each energy term ) 
 using the nm and skempi experimental sets .
we assessed the accuracy 
 of predictions ( tables 2 and 3 ) and found that energy function pred1 had the best agreement 
 with experiments for both datasets and had only four parameters and 
 three energy terms , all of which were derived from charmm:4here evdw is the change 
 of van der waals interaction energy and gsolv is the change of polar solvation energy 
 of solute in water .
samut represents a term proportional 
 to the interface area of the mutant complex ( if we choose the wild - type 
 complex interface area , the results are not significantly different 
 from those where mutant interface area is used ) .
nonpolar solvation 
 energy and other terms did not contribute significantly to the model 
 quality ( p - value > 0.01 ) and were not used in 
 our 
 model . for the pb calculation
dielectric 
 constants 
  = 1 , 2 , 4 and different ion concentrations were tested on 
 single mutations of nm set using the optimized minimization protocol 
 and energy function pred1 ( the testing results are shown in table 
 s2 , supporting information ) . as a result
, 
  = 2 for the protein interior ,  = 80 for the exterior 
 aqueous environment and ion concentration of zero were used for further 
 energy calculations .
all pb calculations were performed with the pbeq 
 module ; the van der 
 waals interaction energy was calculated with the energy module , and 
 molecular surface area was obtained with the sasa module of the charmm 
 program .
the atomic born radii were previously 
 calibrated and optimized to reproduce the electrostatic free energy 
 of the 20 amino acids in md simulations with explicit water molecules .
the foldx method calculates the effect 
 of mutations on protein stability using an empirical force field . 

it optimizes the side chain configurations without taking into account 
 the backbone conformational movements .
we tested two different protocols 
 using foldx : the repairpdb module applied to protein crystal structures 
 and the repairpdb module applied to the minimized structures .
we choose 
 the first foldx protocol for comparison as it gave us better correlation 
 with experimental values .
to calculate binding energy by foldx , we 
 used eq 2 and calculated the unfolding free 
 energy of the whole complex and each monomer separately .
the 
 second method , beatmusic is specifically 
 trained to calculate the change in binding affinity produced by single 
 missense mutations .
it uses residue based statistical potentials that 
 were previously optimized to discriminate native from decoy complexes .
in addition , the beatmusic energy function has 
 terms accounting for packing defects and for the effect of mutations 
 on the unfolding free energy of the whole complex in cases where unbound 
 monomers could be considered disordered .
the 
 third method , cc / pbsa , applies the concoord 
 approach to sample the protein configurational ensemble and the pb 
 approach to calculate the solvent contribution to the binding free 
 energy .
we used two main measures to estimate 
 the quality of the agreement between experimental and predicted values . 

first , we calculated the pearson correlation coefficient and tested 
 a null hypothesis about the equality of the correlation coefficient 
 to zero .
all correlation coefficients reported in the paper were significantly 
 different from zero with p - values of less than 0.01 . 

another measure was the root - mean - square error ( rmse ) , which is the 
 square root of the sum of the differences between values predicted 
 by a model and experimental values divided by the number of cases . 

five - fold cross validation was done on training sets where the data 
 was randomly partitioned into five subsets of approximately equal 
 size , and one subset was used for testing the model and the remaining 
 four subsets were used as training data . to diminish similarity between 
 training and testing sets and ensure the accuracy of cross - validated 
 correlation coefficients , protein complexes from the same similarity 
 cluster , provided by the skempi database , were used only for training 
 and not for testing .
our goal is 
 to construct a computational protocol that would yield good prediction 
 accuracy for a diverse and large set of missense mutations .
we obtained 
 a new energy function with weighting factors in front of each terms 
 based on eq 3 ( see methods ) by applying a multiple linear regression procedure .
this model 
 was fit to the experimental data of the differences in binding affinities 
 produced by mutations .
because several mutations can interact with 
 each other and exhibit cooperativity and positive epistasis , we performed the fitting procedure separately 
 for single and multiple mutations .
our results showed that the separate 
 parameter optimization for single and multiple mutations yielded a 
 better agreement with experiments than if we combined them .
table 1 ( table s2 , supporting information , has more detailed information ) shows results for different simulation 
 protocols for 242 single mutants from the nm set with the 5-fold cross - validation 
 using the pred1 energy function .
as can be seen from table 1 , minimization in explicit water without restrains 
 on the backbone atoms shows the highest correlation which indicates 
 how important the water model is for the quality of prediction .
moreover , 
 minimizing the structures of mutants results in a better quality of 
 g estimates than a protocol where the 
 vmd in silico model is used without minimization ( see zero minimization 
 step in figure s3a  d , supporting information ) .
however , increasing the number of minimization steps after about 
 300 does not significantly improve the average prediction performance . 

the exception is the minimization for multiple mutations from the 
 skempi set where the shorter minimization procedure in general produces 
 better results ( figure s3d , supporting information ) .
r , the pearson correlation 
 coefficient between experimental and predicted g values , r , the five - fold 
 cross - validated correlation and rmse , the root - mean squared error 
 are each shown for the case of training / testing on single mutations 
 of nm set .
this table only shows the results with dielectric constant 
 of two and ion concentration of zero .
we executed 1 ns md simulations 
 for 242 single mutants from nine 
 protein complexes of the nm set .
no significant conformational changes 
 were observed in most cases in terms of rmsd from the starting mutant 
 complex that was minimized for 1000 steps prior to md simulation . 

relatively large conformational changes between the wild - type and 
 mutant proteins were observed in a few cases ( rmsd of more than 2.5 
  ) , although the average backbone rmsd was less than 1.5  
 ( figure s5 , supporting information ) .
we 
 observed even smaller backbone conformational changes during the minimization 
 procedure ( data not shown ) .
as can be seen from table 1 and figure s3(e , f ) ( supporting information ) , md simulation over 1 ns yields a maximum correlation of 0.40 between 
 experiments and predictions for cases when md simulations were done 
 for both mutant and wild - type complexes ( figure s3f , supporting information ) . nevertheless , it is considerably 
 lower compared to the protocol when only minimization procedure is 
 used .
tables 2 and 3 compare 
 our approach 
 with other methods for single and multiple mutations , respectively . 

first , we tested two energy functions : pred3 and pred4 with 11 and 
 13 parameters for single and multiple mutations , respectively ( energy 
 functions are defined in the supporting information ) .
pred4 and pred3 were applied to single and multiple mutations 
 of the nm set ( tables 2 and 3 ) and produced very good agreement with the experiment ( cross - validated 
 correlation coefficient of 0.68 and 0.85 for pred4 and pred3 , respectively , 
 applied to single and multiple mutations ) , but they showed rather 
 poor results on the large representative set of skempi mutations with 
 the correlation coefficients of 0.37 and 0.46 for single and multiple 
 mutations , respectively .
it should be mentioned that , unlike the correlation 
 coefficient , rmse values are scale dependent and can only be compared 
 for different methods for the same test set , not between different 
 sets .
it is known that if parameters are tuned to overminimize mean 
 squared errors and overfit the model , this can lead to decreased generalization 
 performance .
r is the pearson 
 correlation coefficient between experimental and predicted g values and r is the five - fold 
 cross - validated correlation coefficient .
the last column shows the 
 slope of the regression line between predicted and experimental values . 

r is the pearson 
 correlation coefficient between experimental and predicted g values and r is the five - fold 
 cross - validated correlation coefficient .
the last column shows the 
 slope of the regression line between predicted and experimental values . 

taking all this into account , we have constructed 
 a model and ensured 
 that it was trained on one large representative set ( skempi ) and yielded 
 a good agreement with experiments when tested on another independent 
 set ( nm ) .
our energy function pred1 ( see methods , eq 4 ) had four parameters and three energy 
 terms .
energy terms with statistically significant contribution to 
 the quality of multiple regression model are listed in table s3 ( supporting information ) together with their corresponding 
 coefficients / weights . as can be seen from tables 2 and 3
, the pred1 function gives a 
 correlation coefficient of 0.62 and 0.63 and rmse values of 1.27 and 
 1.90 kcal mol , if trained on skempi and tested 
 on the nm set for single and multiple mutations respectively .
interestingly , 
 the energy function pred1 trained on a set of skempi multiple mutations 
 showed better results for multiple mutations compared to using the 
 energy function pred1 trained on skempi single mutations set ( correlation 
 coefficient of 0.55 versus 0.49 ) .
we compared our approach 
 with three other independent methods , foldx , beatmusic and cc / pbsa . 

beatmusic can not be applied to multiple mutations whereas cc / pbsa 
 is not very efficient for large sets like skempi because of its slow 
 computation time . as can be seen from table 2 , beatmusic overall yields better results than foldx whereas the 
 pred1 model gives a better agreement with the experiments compared 
 to both of these methods ( correlation of 0.62 , 0.52 and 0.47 for pred1 , 
 beatmusic and foldx , respectively , for testing on the nm set ) .
limited 
 foldx capacity can be attributed to the fact that the foldx empirical 
 energy function was parametrized on the experimental changes of unfolding 
 free energy , not the binding affinity .
the comparison of their computational 
 speed will be reported in another section . given the relatively high 
 algorithm efficiency of the foldx and beatmusic methods , we decided 
 to combine their scoring schemes with our pred1 energy function.5the performance of the combined model ( pred2 ) 
 exceeded the performances of the individual pred1 ( p - value < 0.01 , see table s3 , supporting information ) , foldx and beatmusic methods and returned the correlation coefficient 
 of 0.69 with an rmse of 1.20 kcal mol .
importantly , 
 the slope of the regression line for pred2 is very close to 1 ( 0.97 ) , 
 which indicates that predicted values are on the same scale as experimental 
 ones .
in addition , we used the cc / pbsa server and managed to obtain 
 the results for 66 single mutations randomly selected from the skempi 
 set . after excluding five mutations with very high predicted energies 
 of gcc / pbsa > 100 kcal 
 mol ,
the correlation between experimental and 
 calculated 
 values for remaining 61 mutations was found to be r = 0.23 ( p - value = 0.07 ) using the cc / pbsa method 
 and r = 0.47 ( p - value < 0.01 ) 
 using the pred2 energy function .
we further compared the accuracy 
 separately for mutations that stabilized ( gpred < 0 ) or destabilized ( gpred > 0 ) the protein  protein binding .
mutations 
 were also subdivided into binding hot spots ( |gpred|  2 kcal mol ) or others ( |gpred| < 
 2 kcal mol ) .
the prediction accuracy was defined 
 as a percentage of correctly identified mutations out of the total 
 number of mutations ( tp + tn)/total , where tp denotes true positives 
 and tn corresponds to true negatives .
sensitivity was defined as tp/(tp 
 + fn ) and specificity was calculated as tn/(tn + fp ) ( fn , false negative ; 
 fp , false positive ) .
table 4 shows that the 
 pred2 energy function has very high prediction accuracy for destabilizing 
 mutants but low accuracy for stabilizing mutants . stabilizing mutants 

are not well predicted by any of the methods used in this study and 
 the foldx method provides the highest prediction accuracy for stabilizing 
 mutants .
the specificity of predictions of binding hot spots is almost 
 100% for all methods although the sensitivity is compromised for all 
 of them .
we did not have enough data to perform a similar analysis 
 for multiple mutations .
figure 1 shows experimental 
 and predicted g for single and multiple 
 mutants for the skempi training set and the nm test set .
we can see 
 that only a few stabilizing mutants are predicted as stabilizing whereas 
 the majority of them are predicted as destabilizing .
training is 
 performed on skempi single ( a , b ) and multiple 
 ( c , d ) mutation sets . a : testing on skempi single mutation set .
beatmusic has the shortest processing 
 time and calculation for one mutation takes less than a second on 
 its webserver ( http://babylone.ulb.ac.be/beatmusic/ ) .
foldx 
 is also very fast and takes about 5 min for a protein of about 300 
 residues when calculations are performed on an intel core duo2 2.8 
 ghz processor .
if we consider the parallel cpu calculation ability , 
 our minimization protocol with namd uses 15 min of cpu time for a 
 protein of 160 residues and 10  water box , altogether of about 
 23  900 atoms on 16 processors ( with a 2.8 ghz intel emt64 ) . 

there is no need to redo the long wild - type minimization procedure 
 if several mutations of one protein are analyzed . in our study
, we 
 utilized the high computational capabilities of the biowulf linux 
 cluster at the national institutes of health .
for the binding free 
 energy calculation using charmm , it takes about 10 min on one processor 
 for about 160 residues protein ( 2.6 ghz amd opteron ) .
the cc / pbsa 
 approach uses the conformational sampling method and has long processing 
 time ( 249 min for 149 residue protein on a 3.2 ghz intel xeon processor , 
 as reported in ref ( 22 ) ) .
however , the cc / pbsa webserver ( http://ccpbsa.biologie.uni-erlangen.de/ccpbsa/ ) takes much longer to process one mutation and is not practically 
 applicable for large datasets . on the basis of our previously mentioned results
, we 
 decided to use skempi as a representative training set for the multiple 
 regression fitting procedure .
the resulting coefficients / weights in 
 front of each energy terms and standardized coefficients in front 
 of standardized variables ( variable with the variance equal to one ) 
 are listed in table s3 ( supporting information ) .
as can be seen from this table , the polar solvation 
 energy term calculated using poisson  boltzmann equation ( gsolv ) and the van der waals ( evdw ) terms have largest contributions to the 
 total change in binding energy gpred1 for single mutants ( standardized coefficients of 0.40 
 and 0.34 , respectively ) .
interestingly , for multiple mutations , gsolv becomes twice as important as for single 
 mutations whereas the role of interface area diminishes .
if we use 
 the pred2 model , then the largest contribution also comes from the 
 polar solvation energy term ( standardized coefficient of 0.31 ) .
the 
 beatmusic score ( gbm ) has 
 a relatively high impact , whereas foldx and interface area make smaller 
 contributions to the quality of the model ( standardized coefficients 
 of about 0.15 ) .
however , all above terms contribute significantly 
 to the agreement between experiments and predictions with p - values of less than 0.01 . in this section , 

we consider several factors that might influence the quality of the 
 prediction for single mutations .
importantly , all of these factors 
 represent features of studied mutation sites that can be derived either 
 from sequences or protein structures a priori . first , 
 we analyzed the effect of the type of residue substitutions on prediction 
 accuracy .
figure 2 shows the boxplots of prediction 
 errors for wild - type and mutant residue types , respectively . in each 
 box ,
the central line is the median , the edges of the box are the 
 25th and 75th percentiles , the whiskers extend to the most extreme 
 data points and outliers are plotted with circles .
the median signed 
 error for most wild - type amino acids is close to zero with two notable 
 exceptions of overprediction , for met and pro . on the other side of 
 the spectrum
are glu , tyr and thr , which produce underprediction errors 
 but not as prominently as met and pro . as to the mutant residue type , 
 pro
the cyclic structure of proline s side chain introduces constraints 
 on the main - chain dihedral angles and can be structurally important 
 for stability or binding .
as can be seen from figure s6 ( supporting information ) , substitutions from / into 
 pro might cause relatively large local conformational changes compared 
 to other substitution types .
moreover , substitutions of other amino 
 acids into proline can result in significant destabilization of complexes 
 ( figure s7 , supporting information ) .
boxplots of 
 prediction error ( residual ) for different wild - type 
 ( a ) and mutant ( b ) residue types for single mutants from skempi training 
 set .
the residual is calculated as a difference between experimental 
 gexp and predicted gpred2 values . as was previously shown
, mutations resulting in amino acid substitutions 
 with similar physico - chemical properties may not drastically alter 
 the stability of a protein or a complex .
we examined 
 the performance of our model pred2 in relation to the change in charge 
 and volume of side chains produced by mutations ( tables s4 and s5 , supporting information ) .
one can see from table 
 s4 ( supporting information ) that the majority 
 of substitutions are neutral - into - neutral , whereas the minority of 
 them change the charge of amino acids .
positive - into - negative and 
 neutral - into - negative amino acid substitutions are characterized by 
 notable correlations ( not enough statistics to estimate negative - into - positive 
 charge changing substitutions ) .
if we consider changes in the amino 
 acid volume , then the highest correlation is observed for substitutions 
 of small - into - large amino acids ( table s5 , supporting 
 information ) .
as can be seen from figure s6 ( supporting information ) , small - into - large substitutions also 
 cause largest conformational changes unlike large - into - small substitutions , 
 which usually do not produce steric clashes .
the 
 location of mutations in protein structures is another important feature , 
 and even neighboring mutations might produce very different effects 
 on binding affinity .
all mutations can 
 be classified into five types depending on their locations with respect 
 to interface and surface ( see ref ( 61 ) and the supporting information for definition ) .
figure 3 shows how mutation 
 location may influence the experimental and predicted ( not shown , 
 but the effect is very similar ) g values . 

mutations located in the core of the interface ( cor ) produce the largest 
 changes in binding energy followed by the interface mutations partially 
 exposed to the solvent ( rim and sup ) .
on the basis of these observations , 
 one might conclude that majority of single mutations from the skempi 
 set do not lead to long - range allosteric effects and only those located 
 directly at the interface regions make the largest contribution to 
 the change in binding energy . nevertheless , the effect of noninterface 
 mutations is not negligible .
cor , rim and sup are the 
 core , rim and support regions of the interface and int and sur are 
 the solvent accessible regions .
the 
 effect of protein flexibility on binding can not be ignored and in 
 some cases can improve or deteriorate the predictions . here
we describe 
 two examples where we account for conformational changes in complex 
 structures when the mutation is introduced .
figure 4 compares two structures obtained by minimization 
 and md simulation for mutant alpha - chymotrypsin a / turkey ovomucoid 
 third domain protein complex ( pdb code : 1cho ) ; the leu15glu mutation is introduced 
 into turkey ovomucoid third domain ( chain i ) .
if we compare the minimized 
 structure and the complex conformation after md simulation , we can 
 see the changed conformation after md simulation where side chains 
 of ser217 and ser190 from alpha - chymotrypsin a protein flip over to 
 form four extra hydrogen bonds with mutated residue glu15 ; moreover , 
 ser195 moves closer to glu15 to form another hydrogen bond ( figure 4 and table s6 , supporting information ) .
this might enhance the interaction between the two partners and 
 results in a smaller destabilizing gpred1 values ( 1.13 kcal mol ) for the 
 md simulation protocol compared to minimization protocol ( 4.32 kcal 
 mol ) and experimental value ( 6.6 kcal mol ) .
this example illustrates that although md simulation 
 does allow for larger conformational changes , these transitions might 
 not take place in a native protein .
difference between conformations of leu15glu 
 mutant for -chymotrypsin 
 a / turkey ovomucoid third domain protein complex ( pdb code : 1cho ) obtained by minimization 
 ( mm ) and md simulation ( md ) .
complex between turkey ovomucoid third 
 domain ( chain i ) and -chymotrypsin a ( chains f and g ) for mm 
 structure are shown in yellow and green , respectively .
complex between 
 chains i and f / g for md are shown in red and blue , respectively .
black 
 dotted lines correspond to the hydrogen bonds formed between glu15 
 of chain i and ser217 , ser190 and ser195 of chain g ( ball - and - stick 
 model ) in the md simulated structure . in the case of interleukin-4/receptor  chain protein 
 ( pdb 
 code : 1iar )
md simulation protocol made a relatively better prediction compared 
 to the minimization protocol ( gexp = 0.43 kcal mol , gpred1_mm = 4.57 and gpred1_md = 1.89 kcal mol ) .
the mutant 
 minimized structure loses seven hydrogen bonds between arg85 and asp67 
 and asp125 compared to the wild - type minimized structure ( table s6 , supporting information ) whereas the mutant md 
 simulated structure loses only three hydrogen bonds compared to the 
 md simulated wild - type structure ( table s6 , supporting 
 information ) .
in this study , we attempt 
 to design an efficient computational 
 protocol in order to estimate the impact of single and multiple missense 
 mutations on protein binding affinity . our analysis showed that the 
 choices of simulation procedure and energy function are both important 
 for achieving accurate predictions .
we demonstrated that using an 
 unconstrained minimization procedure in explicit solvent yielded a 
 better agreement with experiments compared to implicit solvent , which 
 is consistent with a previous study on a smaller dataset we used a modified mm - pbsa energy function and 
 optimized its parameters on the experimental data of several thousands 
 of mutations .
we showed that even with the 5-fold cross - validation , 
 it is possible to attain a very high correlation coefficient ( r = 0.85 ) with the experimental data using a model with 
 eleven parameters .
however , these energy functions do not perform 
 well when applied to an independent set of mutations .
it points to 
 the fact that such a model can be overtrained and biased toward certain 
 groups of proteins or mutations .
therefore , one should be very cautious 
 in interpreting the results using only one limited dataset for fitting 
 and testing even if cross - validation is done .
our approach with 
 an all - atom force field model , explicit water 
 minimization protocol and energy function with only four parameters 
 yielded a reasonable correlation coefficient ( r = 
 0.62 ) and rmse values ( rmse = 1.27 kcal mol ) between 
 experimental and predicted values after training was done on a representative 
 skempi set and testing was performed on an independent and nonoverlapping 
 set of single mutations .
consistent with another study , the mm - pbsa 
 energy function displayed superior performance compared to statistical 
 and empirical potentials .
moreover , including 
 statistical scores from foldx and beatmusic into the model produced 
 significantly better agreement with experiments ( p - value < 0.01 ) with a correlation coefficient of 0.69 and rmse 
 of 1.20 kcal mol .
we showed that the largest significant 
 contribution that explained the largest proportion of experimental 
 data variation came from the polar solvation energy term .
this result 
 once more points to the extreme importance of the solvent model and 
 solvation effects .
we also would like to spotlight the drawback of 
 our model that underestimates the contribution of stabilizing mutations . 

this , in turn , could be the result of an insufficient number of stabilizing 
 amino acid substitutions in the experimental training set .
to 
 assess the effect of multiple mutations , we investigated two 
 models trained on single and multiple mutations , respectively .
interestingly , 
 we found that the model trained on a set of multiple mutations showed 
 better results for multiple mutations compared to the model trained 
 on single mutations .
indeed , nonadditivity of the effect of multiple 
 mutations on binding was observed previously for alanine scanning 
 and disease mutations .
it was partially attributed to 
 the modular structure of the binding interface and cooperativity of 
 amino acids within each binding site cluster .
it is difficult to assess the conformational changes of monomers 
 occurring upon binding and the effect of flexibility on prediction 
 accuracy .
although we did not estimate conformational changes upon 
 binding ( in many cases , no unbound states were available ) , we tried 
 to account for conformational changes in complex structures produced 
 by mutations .
as we showed , conformational sampling of mutant structures 
 by molecular dynamics simulations on a 1 ns time scale did not help 
 to achieve better agreement with experiments compared to the protocol 
 when only one minimized mutant complex structure was used .
this is 
 consistent with several previous studies demonstrating that averaging 
 over an md - generated ensemble of conformations had a negligible effect 
 on the quality of binding affinity predictions , on estimating the 
 effect of mutations on binding and on homology model 
 refinement .
it was previously shown that 
 backbone sampling can do more harm than good when estimating the effects 
 of mutations on protein stability in cases where structural changes 
 are negligible .
it implies that either 
 our single minimized conformation dominated a native conformational 
 ensemble or we sampled inaccessible conformations during the unconstrained 
 md simulations . in any case , even without conformational sampling , 
 we were able to achieve similar or better prediction accuracy on a 
 much faster time scale compared to one of the best currently available 
 method ( cc / pbsa ) .
this is especially 
 important for high - throughput virtual screening studies because the 
 conformational sampling is very time - consuming .
finally , we 
 found that the highest correlation was observed for 
 substitutions that changed either charge or side chain volume .
we 
 also observed that missense mutations located directly in the core 
 region of interfaces had the largest effect on binding affinity and 
 long - range effects from noninterface mutations were relatively minor . 

all these distinctive realistic patterns and the reasonable agreement 
 with experiments validates the use of atomic force fields , statistical 
 potentials or combinatorial scoring schemes for estimating the effects 
 of single and multiple missense mutations on protein 
 the 
 crucial prerequisite for proper biological function is the 
 protein s ability to establish highly selective interactions 
 with macromolecular partners . a missense mutation that alters the 
 protein binding affinity may cause significant perturbations or complete 
 abolishment of the function , potentially leading to diseases . 
 the 
 availability of computational methods to evaluate the impact of mutations 
 on protein  protein binding is critical for a wide range of 
 biomedical applications . here , we report an efficient computational 
 approach for predicting the effect of single and multiple missense 
 mutations on protein  protein binding affinity . 
 it is based 
 on a well - tested simulation protocol for structure minimization , modified 
 mm - pbsa and statistical scoring energy functions with parameters optimized 
 on experimental sets of several thousands of mutations . 
 our simulation 
 protocol yields very good agreement between predicted and experimental 
 values with pearson correlation coefficients of 0.69 and 0.63 and 
 root - mean - square errors of 1.20 and 1.90 kcal mol1 for single and multiple mutations , respectively . compared with other 
 available methods , 
 our approach achieves high speed and prediction 
 accuracy and can be applied to large datasets generated by modern 
 genomics initiatives . in addition , we report a crucial role of water 
 model and the polar solvation energy in estimating the changes in 
 binding affinity . 
 our analysis also reveals that prediction accuracy 
 and effect of mutations on binding strongly depends on the type of 
 mutation and its location in a protein complex . 
Introduction
Methods
Results
Discussion