Accurate and Efficient Calculation of Protein ... - ACS Publications

Nov 15, 2018 - NYU−ECNU Center for Computational Chemistry at NYU Shanghai, ... Department of Chemistry, New York University, New York, New York ...
1 downloads 0 Views 596KB Size
Subscriber access provided by Bibliothèque de l'Université Paris-Sud

Computational Biochemistry

Accurate and Efficient Calculation of Protein-Protein Binding Free EnergyInteraction Entropy with Residue Type Specific Dielectric Constants Xiao Liu, Long Peng, and John Z.H. Zhang J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.8b00248 • Publication Date (Web): 15 Nov 2018 Downloaded from http://pubs.acs.org on November 18, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Accurate and Efficient Calculation of Protein-Protein Binding Free Energy-Interaction Entropy with Residue Type-Specific Dielectric Constants Xiao Liu1#, Long Peng1#, John Z.H. Zhang1,2,3* 1State

Key Laboratory for Precision Spectroscopy, Shanghai Engineering Research Center of Molecular Therapeutics & New Drug Development, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai 200062, China 2NYU-ECNU

Center for Computational Chemistry at NYU Shanghai, Shanghai 200062, China

3Department

of Chemistry, New York University, NY, NY 10003, USA

#These

authors contributed equally to this work. whom correspondence should be addressed: [email protected] Key words: PPI, binding free energy, interaction entropy, dielectric constant, GBSA *To

Abstract Accurate and efficient computation of protein-protein binding free energy remains a grand challenge. In this study, we develop a new strategy to achieve efficient calculation for total protein-protein binding free energies with improved accuracy. The new method combines the recently developed interaction entropy method for efficient computation of entropic change together with the use of residue type-specific dielectric constants in the framework of MM/GBSA to achieve optimal result for protein-protein binding free energies. The new strategy is shown to be computationally efficient and accurate than that using standard MM/GBSA method in which the entropic computation is performed by the normal model approach and the protein interior is represented by the standard dielectric constant (typically set to 1), both in terms of accuracy and computational efficiency. Our study using the new strategy on a set of randomly selected 20 protein-protein binding systems produced an optimal dielectric constant of 2.7 for charged residues and 1.1 for non-charged residues. Using this new strategy, the mean absolute error in computed binding free energies for these 20 selected protein-protein systems is significantly reduced by more than threefold while the computational cost is reduced by more than two orders of magnitude, compared to the result using standard MM/GBSA method with the normal mode approach. Similar improvement in accuracy is confirmed for a test set consisting of 10 1

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

protein-protein systems.

2

ACS Paragon Plus Environment

Page 2 of 28

Page 3 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

INTRODUCTION Protein-protein interaction (PPI) plays a critical role in biology and is of vital importance to many biological processes including control of gene expression, immune response, signal transduction, enzyme inhibition, etc.1-11 Examples include binding of antibody-antigen, protease-inhibitor, hormone-receptor, and other PPIs whose disorders are closely related to human disease. 12 Protein-protein interactions have become increasingly important targets for drug discovery and attracted intensive experimental and computational investigations in biomedical sciences. The increasing numbers of high-resolution structures of protein-protein complexes have provided excellent opportunities for understanding PPI at atomic level. Thus, accurate computation of protein-protein binding free energies is highly desirable as it can lead to reliable prediction of protein-protein associations.13-27 Protein-protein binding is very complicated; it involves complex balance of various interactions including electrostatic interaction, van der Waals interaction, hydrogen bonds, hydrophobic and hydrophilic effects. One can experimentally measure the equilibrium dissociation constant (Kd) for protein-protein association and obtain the Gibbs free energy of binding. However, in theory, accurate and efficient calculation of protein-protein binding free energy is still a grand challenge due to complex nature of various interactions at play. There exists various methods for free energy calculations including free energy perturbation (FEP) 28-31, thermodynamics integration (TI) 32, 33 and MM/PB(GB)SA3440 etc. However, due to the complexity of protein-protein interaction and the size of the interacting proteins, implicit solvent based MMPB/GBSA approach is currently the only viable method for practical computation of protein-protein binding free energy due to its lower computation cost (provided entropy calculation by normal mode method is not activated). The MMPB/GB method for solvation has been successfully applied to predicting the binding free energies 41-44 and the hot spots at protein-protein interfaces4548. MMPB/GBSA uses implicit solvent model to evaluate the solvation energy of the solute but it does not include entropic contribution to the binding free energy. To calculate entropic contribution to the binding free energy, the current approach is to use the standard normal mode method49. However, since the normal mode method is approximate, it often does not give accurate result. Furthermore, computation of entropy using normal mode method is extremely expensive computationally, especially for large systems such as protein-protein interaction systems. As a result, most MMPB/GBSA calculation of protein-protein binding free energy simply skips the calculation of entropy. However, without entropic contribution, the calculated binding energies are generally much large in magnitude than measured experimental values, especially for protein-protein systems. In order to overcome the difficulty in calculating entropic contribution to binding 3

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

free energy, an “Interaction Entropy (IE)” approach has recently been proposed for computing entropic contribution to the binding free energy50. The IE approach is highly efficient in computation as compared to the standard normal mode approach. In addition, it is found that the use of different dielectric constants in computational alanine scanning yields more accurate prediction of hot spot residues than that using the standard dielectric constant of unit51-56. In particular, a recent work has shown that by combining the use of IE approach with residue-specific dielectric constants, more accurate prediction of residue-specific free energies in protein-protein interaction is achieved57. In these studies, larger dielectric constants are needed for charged residues than for neutral residues. The work in ref.57 uses alanine scanning to calculate ΔΔG, i.e., the difference in binding free energy ΔG before and after mutation of a specific residue to alanine. Since only the side chain of interest is affected, the corresponding dielectric constant is optimized to obtain the best agreement in ΔΔG with experiment in ref.57. In the current study, however, we are interested in calculating the total binding free energy ΔG which depends on more properties than in the case of relative binding free energy ΔΔG. Therefore more balanced values of dielectric constants (smaller values) are found to be optimal than in ref.57. In this work, we explore the use of different dielectric constants using the IE approach in direct calculation of total protein-protein binding free energies in order to final optimal values of dielectric constants for use with different residues within the framework of MM/GBSA method for solvation energy. For that purpose, we randomly selected twenty protein-protein interaction systems what have known experimental binding affinities for comparison. Besides comparing computed binding free energies with experiment, we also compared the current results with those computed by using the standard MM/GBSA approach with the normal mode method for entropy computation. RESULTS 1) Optimal dielectric constants for protein residues Proper use of dielectric constant to represent the interior electrostatic environment of protein is important to obtaining more realistic interaction energy in computational study of biomolecule interaction. However, there is no consensus on what the optimal values are for the problems of interest. Here we first investigate the effect of using a single dielectric constant to represent protein interior on the accuracy of the calculated total protein-protein binding free energies. For that purpose, we randomly selected 20 protein-protein interaction systems from the SKEMPI58 database with known experimental binding affinities. Within this dataset, the trypsin-BPTI protein complex (PDB entry: 2FTL) has the highest binding affinity (Kd = 5.88 × 10 ―14 mol/L), and 4

ACS Paragon Plus Environment

Page 4 of 28

Page 5 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

the UEV domain-Ubiquitin complex (PDB entry: 1S1Q) has the lowest binding affinity (Kd =6.35 × 10 ―4 mol/L ).The experimental binding affinities cover an extremely broad range of 14 orders of magnitude. The corresponding coordinates of the studied complexes were downloaded from the SKEMPI database (PDB entries:2J0T, 1FY8, 2FTL, 3TGK, 1UUZ, 1KAC, 1SGP, 2SGP, 2CSO, 1CT4,1Z96, 1F5R,1GL0,1GL1,1KTZ,1LFD, 1S0W, 1S1Q, 1SBB, 1XXM). Table 1 shows that the optimal choice of a common dielectric constant is around ε=1.7, which gives a mean absolute error (MAE) of about 8.2 kcal/mol. in computed binding free energies with respect to experimental values. Obviously, this is not very satisfactory and indicates the need for improvement. Realistically, the optimal value of dielectric constant to use for protein interior should be problem- or method- dependent. In addition, the use of a single dielectric constant to represent everywhere inside protein may not be good enough and the use of protein site-dependent dielectric constant should do better in theoretical calculation. It is already known from earlier works on computational alanine scanning that by using residue-specific dielectric constants, one can significantly improve the accuracy of the calculated residue-specific free energies in PPI systems53-57. Yan, et al.57 recent found that the use of dielectric constants of 1, 3 and 5 (or higher) for nonpolar, polar and charged residues, respectively, are preferred for alanine scanning calculation for PPI systems using the interaction entropy method. Since charged residues generate more polarization nearby and therefore require relatively larger dielectric constants than neutral ones. However, these values from alanine scanning studies are not directly transferrable to the present study in which total protein-protein binding free energies are directly calculated. In order to find optimal values, we tested various combinations of dielectric constants for charged and neutral residues. To reduce the amount of computation, we used a single dielectric constant εn for all neutral (polar and nonpolar) residues and another dielectric constant εc for all charged (positive and negative) residues. A series tests of computation were carried out and the optimal values of dielectric constants are found to be εn=1.1 for neutral residues and εc =2.7 for charged residues. The evidence for this choice of dielectric constants is give in Tables 2 and 3. Table 3 showed that increasing the dielectric constant for polar residue does not improved the result. Thus, for a fixed dielectric constant of 2.7 for charged residues, optimal dielectric residues for neutral residues should not deviate much from unity. This is understood since charged residues are much more sensitive to the use of dielectric constant compared to neutral ones. As shown in these tables, this optimal choice of dielectric constants produces a MAE (mean absolute error) of 2.8 kcal/mol. This is significantly improved from the value of 8.2 kcal/mol. obtained by using a single dielectric constant of 1.7 for all residues as shown in Table 1. Although there is a slight improvement of 5

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

correlation coefficient RP (Pearson Correlation Coefficient) by increasing the values of dielectric constants further, but the MAE also increases as well as shown in Tables 2 & 3. 2) Comparison between IE and normal mode calculations Since the IE is determined by the fluctuation of the interaction energy, which is directly extracted from the MD trajectory, it does not require additional computational cost as required in normal mode calculation. However, the thermal averaging of the int ) is sensitive to the fluctuation of the energy. The proteinexponential term exp( E AB protein interaction energy includes both van der Waals and electrostatic energies and the latter usually accounts for bulk of the energy fluctuation. Thus, the use of dielectric constant will have a large impact on the fluctuation of electrostatic energy. Figure 1 shows the effect of using the optimal value of dielectric constants on interaction energy for the 2J0T system. As shown in Fig. 1, both the interaction energy and its fluctuation are significantly reduced by the dielectric constant compared to the standard case with unit dielectric constant. Of course, this decrease in electrostatic interaction energy is largely compensated by similar decrease in electrostatic solvation energy. In contrast, the interaction energy is only slightly affected by the dielectric constants for the 2S0W system, which indicates that the PPI in this system is dominated by van der Waals energy as shown in Fig. 2. Figure 3 plots time integration of the interaction entropy (-TΔS) as a function of MD time for the 1CSO system. For comparison with the normal mode approach for entropy calculation, we ran five independent MD trajectories for each system and used optimal dielectric constant to calculate binding free energy for every trajectory. We took the average of results from five trajectories as the result and calculated the standard deviation from these five trajectories. The result is listed Table 4. We also show results of individual trajectories for 20 systems in table S1 in the supporting information. For comparison, we also performed standard normal mode calculation (unit dielectric constant) for all 20 PPI systems using 25 snapshots for each system, the detail in table5. More snapshots are tested and the result remains similar to that of 50 &100 snapshots. The detailed result is given in tables S2 & S3. The calculated IE over different MD time is shown in Table S7 in the Supporting Information. For a more detailed examination in error comparison between these two approaches, we plotted in Fig. 4 the signed errors for each of these 20 PPI systems. As shown in Fig. 4, the errors are significantly reduced for all the systems except two. It should be mentioned that the normal mode result is based on the use of unit dielectric constant. A more direct comparison would have been between the two methods with individually optimized dielectric constants. But this could be extremely costly as huge amount of normal mode computations would be required to optimize residue or charge specific 6

ACS Paragon Plus Environment

Page 6 of 28

Page 7 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

dielectric constants. Finally, we compare the computational cost between the present IE approach and the normal mode method. The normal mode calculation involves the average over 25 snapshots evenly extracted from 10 ns MD simulation. The IE calculation of entropy is based on time averaging over 10,000 snapshots evenly extracted from the MD trajectory. Figure 5 plots the computational cost for entropy calculation between the IE and NM methods for each of the 20 PPI systems. The average cost of NM calculation using 25 snapshots for each system is about 184 CPU hours, compared to an average of 1.2 hours using the IE approach. 3) Validation on more systems In order to validate the reliability of our method, we randomly selected additional 10 protein-protein systems from SKEMPI database as a test set (PDB entries: 2C0L, 2NU2, 2NU0, 2OOB, 1HE8, 1FCC, 1FR2,2VLQ, 3BTQ, 2VLO). This test set has no overlaps with the training set. We should mention that some of the systems in our study belong to the same protein types but with different mutations. But that are not different PPI systems with different binding free energies. The calculation result of binding free energies of those 10 PPI systems are listed in Table 6 & 7. From table 4 & 6 we found that the MAE (3.0 kcal/mol.) of the test is close to that of the training set (3.2 kcal/mol.) using the IE method with dielectric constants obtained from the training set. Comparison with the normal mode method shows similar trend as shown in tables 5 & 7 in terms of MAE, but the correlation in Table 7 is better than in Table 5, likely due to the smaller number of systems used in Table 7. CONCLUSION AND DISCUSSION In this study, we presented a new approach using IE with residue-specific dielectric constants within the framework of GBSA for accurate and efficient calculation of protein-protein binding free energy. Based on the computational study on 20 randomly selected PPI systems with known experimental binding affinities, optimal values of dielectric constants are found to be 2.7 for charged residues and 1.1 for neutral residues. Compared to the normal mode method used in the standard MM/GBSA protocol with unit dielectric constant, the accuracy of the computed PPI binding free energies for these 20 PPI systems using the present approach is significantly improved. Test calculation for 10 additional PPI systems validated this conclusion. Furthermore, the present approach is computationally more efficient than the normal mode method, resulting in computational savings by two orders of magnitude. There are numerous sources for errors in the calculation and it is almost impossible 7

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 28

to predict very accurate binding free energies for PPI systems without finding and fixing all the errors. Our current study aims to reduce the errors in predicted binding free energies within the practical MM/GBSA approach. Of course, there are limitations to what we can do with the current approach. For example, the standard error is still relatively large because the total interaction energy is very large due to large electrostatic energies.

METHODS The free energy for protein-protein binding

Gbind

can be expressed as the sum of

“gas-phase” and solvation components,

Gbind  Ggas  Gsolv

(1)

where the gas phase free energy ΔGgas give given by

Ggas  Ei nt  T S ,

(2)

where ΔEint is the protein-protein interaction energy, TΔS is the entropic term. ΔGsolv is the solvation free energy whose calculation will be discussed later34-40.

The

1) Calculation of the entropic term of the gas phase free energy The calculation of the energy term ΔEint in Eq. (2) is relatively straightforward but the computation of the entropic term TΔS is somewhat challenging. The standard approach to computing this term is by the Normal Mode method49. In the normal mode method, one first optimize the protein-protein system to a local minima; secondly, one diagonalizes the Hessian matrix (second order derivatives) to obtain eigenvalues; finally, using these eigenvalues to obtain the partition function analytically and obtain the entropy of system. However, the normal mode calculation for systems with large number of degrees of freedom such as protein-protein system is computational expensive since it involves both diagonalization of large matrices and optimization of the system to local minimum. Recently, an interaction entropy (IE) approach50, 57, 59, 60 was proposed to compute entropic term TΔS. In this approach, the entropic contribution to binding free energy is given by

- T S  KT l n e Ei nt ,

(3)

Where β is 1/kT with k being the Boltzmann constant and Ei nt  Ei nt  Ei nt . The int

thermal average of e E AB

can be carried out by time-averaging along the MD

trajectory, 8

ACS Paragon Plus Environment

Page 9 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

e Ei nt



1

N

e E  Ni 

i nt ( t i

)

,

1

(4) As will be shown later, the above IE method for computing entropic term of the proteinprotein interaction free energy is highly efficient as compared to normal mode method. 2) Calculation of solvation energies At present, a widely popular method for calculating solvation free energy is the implicit solvent MM/GBSA approach 61, 62. We employ the MM/GBSA method in our current study mainly because it is computationally efficient and convenient to modify the dielectric constants using the analytical form. In the MM/GBSA approach, the solvation energy is split into the polar and nonpolar components, GB SASA Gsol  Gsol  Gsol ,

(5)

where the GB (polar term) represents the electrostatic solvation free energy which can be obtained by using the analytic GB model 61, 62 GB Gsol  

qi qj e  f GB  1 1    out  i j f GB( r i j , Ri , Rj ) 2  i n

(6)

where



f GB (rij2 , Ri , R j )  rij2  Ri R j exp(rij2 / 4 Ri R j )



1/ 2

(7)

Here εin is the internal dielectric constant of the solute, εout is the dielectric constant of the solvent, rij is the distance between an atomic pair i and j, Ri and Rj represent, respectively, their effective Born radii and κ is Debye-Huckel electrostatic screening parameter. The second term in Eq. (5) is a surface term which accounts for the cavitation and nonpolar part of the solvation energy and is proportional to the solventaccessible surface area (SASA) SASA Gsol   SASA  

(8)

The coefficients γ and β take the standard values of 0.00542 kcal/(mol•Å2) and 0.92 kcal/mol63. In the present study, a modification of the standard GB method is introduced to account for the use of residue-specific dielectric constants for the protein-protein system. 3) Use of residue-specific dielectric constants The internal dielectric constant of protein is usually set to 1 when we use implicit model such as MM/PBSA and MM/GBSA20, 64, 65, but some studies found that the use of internal dielectric constants greater than 1 could lead to better result 54-56. Here, we explore the use of different internal dielectric constants for different residues in the 9

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 28

hope to find optimal ones for protein-protein binding free energy calculations. To employ residue-specific dielectric constants in GB model, the GB formula in Eq. (6) is modified to take the general form GB sol

G

 1 qi qj 1 e  f GB      2 i j f GB( r i j , Ri , Rj )  i  j out 

 ,  

(9)

where εi and εj are the dielectric constants of the interacting atomic pair i and j whose values depend on the residues the atoms belong to. Similarly, the electrostatic component

Eielnte of

Eielnte 

1 2

the gas phase interaction energy

 a Ab B 



Ei nt

also needs to be modified as

qaqb ,  a b r ab

(13)

where a and b represents, respectively, atom a from residue a in protein A and atom b from residue b in protein B. 4) Molecular dynamics simulation In this study, we used 20 protein-protein binding systems as training set and another 10 systems as test set (thus two data set has no overleap.) from the SKEMPI database for computational study58. All MD simulations were performed using pmemd.MPI and pmemd.cuda in AMBER14 programs package66 with ff14SB67 force field for the protein complex. All the complex are solvated with TIP3P68 water, with the initial closest distance from protein surface to the edge of the box being no less than 12 angstroms, and Na+ or Cl- are added to neutralize the system. The particle mash Ewald (PME)69 is used to treat the long-range electrostatic interactions with periodic boundary conditions (PBC). The non-bonded interactions are truncated with 8 angstroms cutoff. Langevin dynamics was used to control the temperature with a collision frequency of 2.0. The time step is set to 2fs and SHAKE70 was used to constrain bonds involving hydrogen atoms. First, minimizing the system by 5000 steps steepest descent and 5000 steps conjugate gradient with constrained the backbone of protein, and then minimizing the whole system by same flow but free protein’s backbone. Second, the system is slowly heated from 0K to 300K over a period of 50ps with constrained protein. Third, a 20 ns MD is run, the first 10ns was used to equilibrate the system and the remaining 10ns was used for statistical sampling at an interval of every 100fs with a total of 100k snapshots. All these samples were used for IE computation and equally spaced 200 snapshots within this time range is also used for calculating MM/GBSA using a modified version of MMPBSA.py63 in AMBER1466 to include residue-specific dielectric constants. In our calculation, all histidine residues are protonated as HIE In order to compare the computational cost between the IE and normal mode method for entropy calculation, 25 snapshots from 10ns MD simulation are used for entropic calculation using the normal mode49 method. 10

ACS Paragon Plus Environment

Page 11 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Acknowledgement This work was supported by National Key R&D Program of China (Grant no. 2016YFA0501700), National Natural Science Foundation of China (Grant nos. 21433004, 91753103), Shanghai Putuo District (Grant 2014-A-02), Innovation Program of Shanghai Municipal Education Commission (201701070005E00020), and NYU Global Seed Grant. We thank the Supercomputer Center of East China Normal University for providing us computer time.

Supporting Information Detailed individual components of the protein-protein binding energies and other technical details of the calculation.

11

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References 1. Keskin, O.; Gursoy, A.; Ma, B.; Nussinov, R., Principles of protein-protein interactions: What are the preferred ways for proteins to interact? Chem. Rev. 2008, 108, 1225-1244. 2. Pazos, F.; HelmerCitterich, M.; Ausiello, G.; Valencia, A., Correlated mutations contain information about protein-protein interaction. J. Mol. Biol. 1997, 271, 511-523. 3. Sowmya, G.; Breen, E. J.; Ranganathan, S., Linking structural features of protein complexes and biological function. Protein Sci. 2015, 24, 1486-1494. 4. Yook, S. H.; Oltvai, Z. N.; Barabasi, A. L., Functional and topological characterization of protein interaction networks. Proteomics 2004, 4, 928-942. 5. Jones, S.; Thornton, J. M., Principles of protein-protein interactions. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 13-20. 6. Lo Conte, L.; Chothia, C.; Janin, J., The atomic structure of protein-protein recognition sites. J. Mol. Biol. 1999, 285, 2177-2198. 7. Muegge, I.; Schweins, T.; Warshel, A., Electrostatic contributions to protein-protein binding affinities: Application to Rap/Raf interaction. Proteins. 1998, 30, 407-423. 8. Pawson, T.; Nash, P., Protein-protein interactions define specificity in signal transduction. Genes Dev. 2000, 14, 1027-1047. 9. Schreiber, G.; Fersht, A. R., Energetics of protein-protein interactions: analysis of the barnasebarstar interface by single mutations and double mutant cycles. J. Mol. Biol. 1995, 248, 478-86. 10. Schreiber, G.; Fersht, A. R., Rapid, electrostatically assisted association of proteins. Nat. Struct. Biol. 1996, 3, 427-431. 11. Schreiber, G.; Haran, G.; Zhou, H. X., Fundamental Aspects of Protein-Protein Association Kinetics. Chem. Rev. 2009, 109, 839-860. 12. Ivanov, A. A.; Khuri, F. R.; Fu, H. A., Targeting protein-protein interactions as an anticancer strategy. Trends Pharmacol. Sci. 2013, 34, 393-400. 13. Dong, F.; Vijayakumar, M.; Zhou, H. X., Comparison of calculation and experiment implicates significant electrostatic contributions to the binding stability of barnase and barstar. Biophys. J. 2003, 85, 49-60. 14. Elcock, A. H.; Sept, D.; McCammon, J. A., Computer simulation of protein-protein interactions. J. Phys. Chem. B 2001, 105, 1504-1518. 15. Kundrotas, P. J.; Alexov, E., Electrostatic properties of protein-protein complexes. Biophys. J. 2006, 91, 1724-1736. 16. Massova, I.; Kollman, P. A., Computational alanine scanning to probe protein-protein interactions: A novel approach to evaluate binding free energies. J. Am. Chem. Soc. 1999, 121, 8133-8143. 17. Sheinerman, F. B.; Honig, B., On the role of electrostatic interactions in the design of proteinprotein interfaces. J. Mol. Biol. 2002, 318, 161-177. 18. Uetz, P.; Giot, L.; Cagney, G.; Mansfield, T. A.; Judson, R. S.; Knight, J. R.; Lockshon, D.; Narayan, V.; Srinivasan, M.; Pochart, P.; Qureshi-Emili, A.; Li, Y.; Godwin, B.; Conover, D.; Kalbfleisch, T.; Vijayadamodar, G.; Yang, M. J.; Johnston, M.; Fields, S.; Rothberg, J. M., A comprehensive analysis of protein-protein interactions in Saccharomyces cerevisiae. Nature 2000, 403, 623-627. 19. Wang, X. Y.; Ji, C. G.; Zhang, J. Z. H., Glycosylation Modulates Human CD2-CD58 Adhesion 12

ACS Paragon Plus Environment

Page 12 of 28

Page 13 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

via Conformational Adjustment. J. Phys. Chem. B 2015, 119, 6493-6501. 20. Sheinerman, F. B.; Norel, R.; Honig, B., Electrostatic aspects of protein-protein interactions. Curr. Opin. Struct. Biol. 2000, 10, 153-159. 21. von Mering, C.; Krause, R.; Snel, B.; Cornell, M.; Oliver, S. G.; Fields, S.; Bork, P., Comparative assessment of large-scale data sets of protein-protein interactions. Nature 2002, 417, 399-403. 22. Yao, X. X.; Ji, C. G.; Xie, D. Q.; Zhang, J. Z. H., Interaction specific binding hotspots in Endonuclease colicin-immunity protein complex from MD simulations. Sci. China: Chem. 2013, 56, 1143-1151. 23. Brock, K.; Talley, K.; Coley, K.; Kundrotas, P.; Alexovy, E., Optimization of electrostatic interactions in protein-protein complexes. Biophys. J. 2007, 93, 3340-3352. 24. Dong, F.; Zhou, H. X., Electrostatic contribution to the binding stability of protein-protein complexes. Proteins. 2006, 65, 87-102. 25. Hendsch, Z. S.; Sindelar, C. V.; Tidor, B., Parameter dependence in continuum electrostatic calculations: A study using protein salt bridges. J. Phys. Chem. B 1998, 102, 4404-4410. 26. Lee, L. P.; Tidor, B., Barstar is electrostatically optimized for tight binding to barnase. Nat. Struct. Biol. 2001, 8, 73-76. 27. Schreiber, G., Kinetic studies of protein-protein interactions. Curr. Opin. Struct. Biol. 2002, 12, 41-47. 28. Jorgensen, W. L.; Thomas, L. L., Perspective on free-energy perturbation calculations for chemical equilibria. J. Chem. Theory Comput. 2008, 4, 869-876. 29. Kita, Y.; Arakawa, T.; Lin, T. Y.; Timasheff, S. N., Contribution of the surface free energy perturbation to protein-solvent interactions. Biochemistry 1994, 33, 15178-89. 30. Rao, S. N.; Singh, U. C.; Bash, P. A.; Kollman, P. A., Free energy perturbation calculations on binding and catalysis after mutating Asn 155 in subtilisin. Nature 1987, 328, 551-4. 31. Reddy, M. R.; Singh, U. C.; Erion, M. D., Ab initio quantum mechanics-based free energy perturbation method for calculating relative solvation free energies. J. Comput. Chem. 2007, 28, 491-494. 32. Beveridge, D. L.; Dicapua, F. M., Free energy via molecular simulation: applications to chemical and biomolecular systems. Annu. Rev. Biophys. Biophys. Chem. 1989, 18, 431-492. 33. Zacharias, M.; Straatsma, T. P.; Mccammon, J. A., Separation‐shifted scaling, a new scaling method for Lennard‐Jones interactions in thermodynamic integration. J. Chem. Phys. 1994, 100, 9025-9031. 34. Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S. H.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; Donini, O.; Cieplak,P.; Srinivasan, J.; Case, D. A.; Cheatham, T. E. Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc.Chem.Res. 2000, 32, 889-897. 35. Chen, Z.; Baker, N. A.; Wei, G. W., Differential geometry based solvation model I: Eulerian formulation. J. Comput. Phys. 2010, 229, 8231-8258. 36. Kuhn, B.; Gerber, P.; Tanja Schulzgasch, A.; Stahl, M., Validation and Use of the MM-PBSA Approach for Drug Discovery. J. Med. Chem. 2005, 48, 4040-8. 37. Luo, R.; David, L.; Gilson, M. K., Accelerated Poisson-Boltzmann calculations for static and dynamic systems. J. Comput. Chem. 2002, 23, 1244-1253. 38. Massova, I.; Kollman, P. A., Combined molecular mechanical and continuum solvent approach 13

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(MM-PBSA/GBSA) to predict ligand binding. Perspect.Drug Discovery Des. 2000, 18, 113-135. 39. Nicholls, A.; Honig, B., A rapid finite difference algorithm, utilizing successive overrelaxation to solve the Poisson–Boltzmann equation. J. Comput. Chem. 1991, 12, 435–445. 40. W. Rocchia; E. Alexov, a.; B. Honig, Extending the Applicability of the Nonlinear Poisson−Boltzmann Equation:  Multiple Dielectric Constants and Multivalent Ions†. J. Phys. Chem. B 2001, 105, 6507-6514. 41. Wang, W.; Kollman, P. A., Free energy calculations on dimer stability of the HIV protease using molecular dynamics and a continuum solvent model 1. J. Mol. Biol. 2000, 303, 567-582. 42. Hou, T.; Chen, K.; Mclaughlin, W. A.; Lu, B.; Wang, W., Computational Analysis and Prediction of the Binding Motif and Protein Interacting Partners of the Abl SH3 Domain. PLoS Comput. Biol. 2006, 2, e1. 43. Gohlke, H.; Case, D. A., Converging free energy estimates: MM-PB(GB)SA studies on the protein-protein complex Ras-Raf. J. Comput. Chem. 2004, 25, 238-50. 44. Chen, F.; Liu, H.; Sun, H.; Pan, P.; Li, Y.; Li, D.; Hou, T., Assessing the performance of the MM/PBSA and MM/GBSA methods. 6. Capability to predict protein–protein binding free energies and re-rank binding poses generated by protein–protein docking. Phys. Chem. Chem. Phys. 2016, 18, 22129-22139. 45. Ylilauri, M.; Pentikäinen, O. T., MMGBSA as a tool to understand the binding affinities of filamin-peptide interactions. J. Chem. Inf. Model. 2013, 53, 2626-33. 46. Bradshaw, R. T.; Patel, B. H.; Tate, E. W.; Leatherbarrow, R. J.; Gould, I. R., Comparing experimental and computational alanine scanning techniques for probing a prototypical proteinprotein interaction. Protein Eng., Des. Sel. 2011, 24, 197-207. 47. Sørensen, J.; Palmer, D. S.; Schiøtt, B., Hot-spot mapping of the interactions between chymosin and bovine κ-casein. J. Agric. Food Chem. 48. Fulle, S.; Withers-Martinez, C.; Blackman, M. J.; Morris, G. M.; Finn, P. W., Molecular determinants of binding to the Plasmodium subtilisin-like protease 1. J. Chem. Inf. Model. 2013, 53, 573-583. 49. Nguyen, D. T.; Case, D. A., On finding stationary states on large-molecule potential energy surfaces. J. Phys. Chem. 1985, 89, 4020-4026. 50. Duan, L. L.; Liu, X.; Zhang, J. Z. H., Interaction Entropy: A New Paradigm for Highly Efficient and Reliable Computation of Protein-Ligand Binding Free Energy. J. Am. Chem. Soc. 2016, 138, 5722-5728. 51. Hou, T.; Wang, J.; Li, Y.; Wang, W., Assessing the performance of the MM/PBSA and MM/GBSA methods: I. The accuracy of binding free energy calculations based on molecular dynamics simulations. J. Chem. Inf. Model. 2011, 51, 69-82. 52. Sun, H.; Li, Y.; Shen, M.; Tian, S.; Xu, L.; Pan, P.; Guan, Y.; Hou, T., Assessing the performance of MM/PBSA and MM/GBSA methods. 5. Improved docking performance using high solute dielectric constant MM/GBSA and MM/PBSA rescoring. Phys. Chem. Chem. Phys. 2014, 16, 22035-45. 53. Xu, L.; Sun, H.; Li, Y.; Wang, J.; Hou, T., Assessing the Performance of MM/PBSA and MM/GBSA Methods. 3. The Impact of Force Fields and Ligand Charge Models. J. Phys. Chem. B 2013, 117, 8408-21. 54. Martins, S. A.; Perez, M. A.; Moreira, I. S.; Sousa, S. F.; Ramos, M. J.; Fernandes, P. A., Computational Alanine Scanning Mutagenesis: MM-PBSA vs TI. J. Chem. Theory Comput. 2013, 14

ACS Paragon Plus Environment

Page 14 of 28

Page 15 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

9, 1311. 55. Simões, I. C.; Costa, I. P.; Coimbra, J. T.; Ramos, M. J.; Fernandes, P. A., New Parameters for Higher Accuracy in the Computation of Binding Free Energy Differences upon Alanine Scanning Mutagenesis on Protein-Protein Interfaces. J. Chem. Inf. Model. 2016. 56. And, I. M.; Kollman, P. A., Computational Alanine Scanning To Probe Protein−Protein Interactions:  A Novel Approach To Evaluate Binding Free Energies. J.Am.Chem.Soc. 1999, 121, 8133-8143. 57. Yan, Y.; Yang, M. Y.; Ji, C. G.; Zhang, J. Z. H., Interaction Entropy for Computational Alanine Scanning. J. Chem. Inf. Model. 2017, 57, 1112-1122. 58. Moal, I. H.; Fernandez-Recio, J., SKEMPI: a Structural Kinetic and Energetic database of Mutant Protein Interactions and its use in empirical models. Bioinformatics 2012, 28, 2600-2607. 59. Sun, Z. X.; Yan, Y. N.; Yang, M. Y.; Zhang, J. Z. H., Interaction entropy for protein-protein binding. J. Chem. Phys. 2017, 146. 60. Qiu, L.; Yan, Y.; Sun, Z.; Song, J.; Zhang, J. Z. H., Interaction entropy for computational alanine scanning in protein–protein binding. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2018,8,e1342. 61. Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T., Semianalytical treatment of solvation for molecular mechanics and dynamics. J. Am. Chem. Soc. 1990, 112, 6127-6129. 62. Srinivasan, J.; Trevathan, M. W.; Beroza, P.; Case, D. A., Application of a pairwise generalized Born model to proteins and nucleic acids: inclusion of salt effects. Theor. Chem. Acc. 1999, 101, 426-434. 63. Miller, B. R.; McGee, T. D.; Swails, J. M.; Homeyer, N.; Gohlke, H.; Roitberg, A. E., MMPBSA.py: An Efficient Program for End-State Free Energy Calculations. J. Chem. Theory Comput. 2012, 8, 3314-3321. 64. Schutz, C. N.; Warshel, A., What are the dielectric "constants" of proteins and how to validate electrostatic models? Proteins. 2001, 44, 400-417. 65. Wagoner, J.; Baker, N. A., Solvation forces on biomolecular structures: a comparison of explicit solvent and Poisson-Boltzmann models. J. Comput. Chem. 2004, 25, 1623-9. 66. Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J., The Amber biomolecular simulation programs. J. Comput. Chem. 2005, 26, 1668-1688. 67. Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C., ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696-3713. 68. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L., Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926-935. 69. Darden, T.; York, D.; Pedersen, L., Particle mesh Ewald: An N⋅log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089-10092. 70. Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C., Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n -alkanes. J. Comput. Phys. 1977, 23, 327-341.

15

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 28

Table 1 MAE (mean absolute error) of computed binding free energies of 20 proteinprotein systems computed by the IE method with MM/GBSA for solvation using different values of DC (dielectric constant) for all residues. The MAE is with respect to experimental values in unit of kcal/mol.

DC

1.5

1.6

1.7

1.8

1.9

2.0

2.1

MAE

12.9

9.6

8.2

8.8

10.7

13.1

15.8

16

ACS Paragon Plus Environment

Page 17 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Table 2 Three different error estimates for binding free energies of 20 protein-protein systems between experimental and computed values using IE method with MM/GBSA for solvation using different values of DCc (dielectric constant for charged residues) while neutral residues were given the dielectric constant of 1.1. MAE is mean absolute error; RMSE, root means square error; and RP, the Pearson correlation coefficient. The MAE and RMSE are all in unit of kcal/mol. DCC

MAE

RMSE

RP

2.2

6.3

7.1

0.68

2.3

5.0

5.8

0.71

2.4

4.0

4.7

0.73

2.5

3.3

4.0

0.75

2.6

2.9

3.8

0.77

2.7

2.8

4.0

0.78

2.8

3.3

4.4

0.79

2.9

4.0

5.0

0.80

3.0

4.6

5.6

0.80

3.1

5.2

6.2

0.80

3.2

5.8

6.8

0.80

17

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 28

Table 3 The same as in Table 2 but for variations of DCp (dielectric constant of polar residues) or DCn (dielectric constant of nonpolar residues) while keeping the dielectric constant of 2.7 for charged residues. DCP (DCN=1.1)

MAE

RMSE

RP

1.0

3.7

4.3

0.71

1.1

2.8

4.0

0.78

1.2

4.8

5.9

0.81

1.3

6.9

7.9

0.81

1.4

8.9

9.8

0.81

1.5

10.4

11.5

0.80

1.0

3.2

3.8

0.75

1.1

2.8

4.0

0.78

1.2

3.6

4.7

0.80

1.3

4.6

5.5

0.81

1.4

5.4

6.4

0.81

1.5

6.2

7.1

0.81

DCN (DCP=1.1)

18

ACS Paragon Plus Environment

Page 19 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Table 4. Decomposition of binding free energies of 20 protein-protein systems computed by IE method and comparison to the experimental binding free energies. The optimal values of 2.7 and 1.1 of dielectric constants are used for all charged and neutral residues, respectively. vdw E AB

ele E AB

is the electrostatic interaction energy between proteins;

the van der Waals energy; ΔGGB the GB solvation energy; ΔGcal the calculated

binding free energy; ΔGexp the experiment binding free energy from ref 58. The experiment binding affinities and temperatures measured at are listed in TableS4 in SI. MAE is mean absolute error; RMSE root means square error; and RP is the Pearson correlation coefficient; STD is the sample standard deviation; SEM is the standard error of the mean and MEAN STD is the mean sample standard deviation values of the 20 systems. All values are in unit of kcal/mol. SEM

PDB

ele E AB

vdw E AB

ΔGGB

ΔH

-TΔS

ΔGcal

ΔGexp

STD

2J0T

-380.3

-90.6

395.3

-75.5

65.2

-10.3

-12.1

5.7

2.5

1FY8

-225.1

-72.9

233.1

-64.9

56.4

-8.5

-6.9

4.6

2.1

2FTL

57.6

-88.0

-38.0

-68.4

45.9

-22.5

-18.2

3.1

1.4

3TGK

-193.0

-76.7

204.0

-65.7

55.0

-10.6

-6.8

8.3

3.7

1UUZ

-281.3

-105.3

306.2

-80.3

67.6

-12.7

-12.4

3.6

1.6

1KAC

-124.4

-75.2

151.7

-47.9

43.0

-4.9

-9.0

3.0

1.3

1SGP

-88.7

-80.4

108.4

-60.6

46.1

-14.6

-11.7

7.7

3.4

2SGP

-70.8

-76.8

92.9

-54.8

43.8

-11.0

-6.3

8.1

3.6

1CSO

-95.4

-86.9

113.0

-69.3

51.6

-17.8

-10.2

3.7

1.7

1CT4

-96.4

-84.9

113.6

-67.8

51.6

-16.2

-11.7

8.6

3.8

1E96

-125.5

-64.6

133.9

-56.1

45.6

-10.5

-9.5

7.5

3.4

1F5R

-228.4

-69.4

235.5

-62.4

55.0

-7.4

-5.8

4.7

2.1

1GL0

-0.9

-95.4

31.0

-65.3

50.3

-15.0

-11.9

14.6

6.5

1GL1

-10.7

-90.1

35.5

-65.3

47.4

-17.9

-13.3

9.1

4.1

1KTZ

-164.8

-46.9

163.6

-48.1

36.8

-11.3

-8.9

3.6

1.6

1LFD

-140.0

-61.2

143.4

-57.7

50.8

-6.9

-7.8

5.3

2.4

1S0W

-65.9

-110.0

107.3

-68.6

56.9

-11.7

-9.3

1S1Q

-74.5

-64.7

94.6

-44.6

44.2

-0.4

-4.4

5.7

2.5

1SBB

-46.5

-65.1

70.8

-40.8

41.7

1.0

-5.3

3.6

1.6

1XXM

-80.5

-98.8

116.4

-62.8

51.2

-11.6

-9.0

4.1

1.8

MAE

3.2

RMSE

3.7

RP

0.8

19

ACS Paragon Plus Environment

MEAN STD

7.3

6.1

3.3

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 28

Table 5 Same as table 4 but using normal mode method (25 snapshots) and unit dielectric constant to calculating binding free energies. PDB 2J0T 1FY8 2FTL 3TGK 1UUZ 1KAC 1SGP 2SGP 1CSO 1CT4 1E96 1F5R 1GL0 1GL1 1KTZ 1LFD 1S0W 1S1Q 1SBB 1XXM

ele E AB

vdw E AB

ΔGGB

ΔH

-TΔS

ΔGcal

ΔGexp

STD

-731.5 -442.6 233.7 -396.4 -625.7 -263.2 -118.3 -112.5 -155 -187.9 -306.9 -564.6 24.6 35.0 -444.7 -293.4 -52.3 -154 -61.8 -130.4

-83.6 -66.7 -92.2 -80.0 -113.5 -75.6 -81.4 -76.2 -83.7 -93.2 -65.8 -64.8 -101.8 -89.5 -46.6 -60.1 -108.6 -63.3 -69.3 -92.7

736.3 449.1 -217.1 417.6 656.7 294.5 144.4 134.1 167.0 204.4 305.0 564.3 2.0 -2.3 437.4 281.0 112.3 172.4 96.3 153.8

-67.9 -52.2 -64.8 -48.8 -67.8 -33.1 -45.8 -46.0 -62.5 -66.1 -57.8 -56.6 -63.2 -46.4 -47.0 -62.6 -32.8 -36.0 -25.5 -55.1

43.5 42.8 42.5 39.3 52.6 39.7 33.7 32.7 37.0 37.9 36.8 42.1 44.1 44.2 33.4 50.5 57.3 38.4 36.9 60.4

-24.4 -9.4 -22.3 -9.5 -15.3 6.5 -12.1 -13.3 -25.4 -28.2 -21.1 -14.5 -19.1 -2.3 -13.6 -12.1 24.5 2.4 11.4 5.3

8.1 8.3 8.3 9.8 9.9 8.3 9.6 7.2 7.8 7.0 7.9 8.1 10.7 8.2 6.1 8.0 10.0 7.0 5.8 9.2

MAE

9.9

-12.1 -6.9 -18.2 -6.8 -12.4 -9.0 -11.7 -6.3 -10.2 -11.7 -9.5 -5.8 -11.9 -13.3 -8.9 -7.8 -9.3 -4.4 -5.3 -9.0 MEAN STD

RMSE RP

12.4 0.4

20

ACS Paragon Plus Environment

8.3

Page 21 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Table 6 Same as Table 4 but using an independent test set which has 10 protein-protein systems random selected from SKEMPI database to validate our optimal dielectric is suitable or not. PDB

ele E AB

vdw E AB

ΔGGB

ΔH

-TΔS

ΔGcal

ΔGexp

2C0L 2NU2 2NU0 2OOB 1HE8 1FCC 1FR2 2VLQ 3BTQ 2VLO

-258.0 -113.0 -75.1 -16.1 -52.3 -151.6 -438.5 -569.9 19.2 -543.1

-103.2 -90.1 -87.5 -45.5 -72.2 -64.6 -83.7 -82.9 -88.7 -79.9

276.6 132.2 94.4 29.4 67.6 153.9 442.7 568.7 3.3 547.6

-84.7 -70.9 -68.2 -32.2 -56.8 -62.3 -79.4 -84.2 -66.2 -75.4

72.6 58.0 58.0 28.6 52.3 57.0 61.0 65.8 53.4 63.4 MAE RMSE RP

-12.07 -12.87 -10.2 -3.55 -4.54 -5.3 -18.45 -18.39 -12.81 -11.94 3.0 3.2 0.8

-9.56 -11.4 -12.9 -5.9 -8.1 -9.1 -16.8 -14.6 -8.7 -16.5

21

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 28

Table 7 Same as table 6 but using normal mode method and unit dielectric constant to calculating binding free energies. PDB

ele E AB

vdw E AB

ΔGGB

ΔH

-TΔS

ΔGcal

ΔGexp

2C0L 2NU2 2NU0 2OOB 1HE8 1FCC 1FR2 2VLQ 3BTQ 2VLO

-571.1 -204.6 -109.1 -24.2 -159.4 -309.9 -1070.3 -1414.6 174.2 -1245.6

-103.2 -90.2 -87.5 -45.5 -72.2 -64.6 -83.7 -82.9 -88.7 -79.9

595.5 237.2 139.1 47.2 180.8 321.0 1082.7 1416.0 -142.1 1255.6

-78.9 -57.6 -57.5 -22.4 -50.8 -53.5 -71.3 -81.5 -56.6 -69.9

60.3 39.4 36.4 25.5 58.5 43.1 43.4 47.2 41.0 44.1 MAE RMSE RP

-18.6 -18.2 -21.1 3.1 7.6 -10.4 -27.9 -34.2 -15.6 -25.8 9.7 10.8 0.8

-9.6 -11.4 -12.9 -5.9 -8.1 -9.1 -16.8 -14.6 -8.7 -16.5

22

ACS Paragon Plus Environment

Page 23 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 1 Inter-protein interaction energy in gas phase (left) and its distribution (right) for the protein-protein system 2J0T from the final 10 ns MD simulation. The upper figures are the result using standard dielectric constant of unit for all the residues and the lower figures are the result using optimal dielectric constant of 2.7 for charged residues and 1.1 for neutral residues.

23

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2 Same as figure 1 but for the PPI system of 2S0W.

24

ACS Paragon Plus Environment

Page 24 of 28

Page 25 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 3. The integration curve of IE (kcal/mol) as a function of MD time for the PPI systems 1CSO, 1E96, 1F5R, 1KTZ and 1LFD.

25

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4. Deviations of calculated binding free energies from their experimental values using the IE method with optimal dielectric constants and the standard normal mode method with unit dielectric constant for 20 PPI systems. The deviation is defined as ΔΔG= ΔGexp - ΔGcal.

26

ACS Paragon Plus Environment

Page 26 of 28

Page 27 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Figure 5 Computation cost (core hours) of normal mode and IE methods to calculate entropy. In this study, the CPU cost of normal mode calculation is for 25 snapshots informally extracted from 10ns MD simulation. The CPU cost is for the time integration of 100k configurations uniformly extracted from 10ns MD simulation. The X20 in the legend represent the cpu hours of IE are multiplied by a factor of 20 before they were plotted. The computation is performed on Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.40GHz platform.

27

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11

ACS Paragon Plus Environment

Page 28 of 28