Improved Computation of Protein–Protein Relative Binding Energies

Aug 8, 2016 - A MMGBSA variant (here referred to as Nwat-MMGBSA), based on the inclusion of a certain number of explicit water molecules (Nwat) during...
0 downloads 12 Views 3MB Size
Subscriber access provided by Northern Illinois University

Article

Improved Computation of Protein-Protein Relative Binding Energies with the Nwat-MMGBSA Method Irene Maffucci, and Alessandro Contini J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.6b00196 • Publication Date (Web): 08 Aug 2016 Downloaded from http://pubs.acs.org on August 16, 2016

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 free 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 accessible to all readers and 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.

Journal of Chemical Information and Modeling 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 46

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

Improved Computation of Protein-Protein Relative Binding Energies with the Nwat-MMGBSA Method Irene Maffucci and Alessandro Contini* Dipartimento di Scienze Farmaceutiche – Sezione di Chimica Generale e Organica “Alessandro Marchesini”, Università degli Studi di Milano, Via Venezian, 21 20133 Milano, Italy KEYWORDS. Protein-protein interaction, molecular dynamics, MMGBSA, water, force field, interface, implicit solvent model.

ABSTRACT. A MMGBSA variant (here referred to as Nwat-MMGBSA), based on the inclusion of a certain number of explicit water molecules (Nwat) during the calculations, has been tested on a set of 20 protein-protein complexes, using the correlation between predicted and experimental binding energy as the evaluation metric. Beside the Nwat parameter, the effect of the force field, the molecular dynamics simulation length, and the implicit solvent model used in the MMGBSA analysis have been also evaluated. We found that considering 30 interfacial water molecules improved the correlation between predicted and experimental binding energies by up to 30%, compared to the standard approach. Moreover, the correlation resulted rather sensitive to the force field and, to a minor extent, to the implicit solvent model, and to the length of the MD simulation.

Introduction

ACS Paragon Plus Environment

1

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 2 of 46

Protein-protein interactions (PPIs) are central to the regulation of important biological processes,1,2 such as cell proliferation, growth, differentiation, signal transduction, and apoptosis.3–6 As a consequence, they can be involved in many disease states, such as cancer, neurodegeneration, and viral and bacterial infections.7–9 Therefore, the modulation of PPIs has a great therapeutic potential and, in the last decades, targeting PPIs known to be involved in pathologic states became an appealing goal for drug design/discovery.1,3,7,10–14 However, the design of effective PPI modulators still represents a challenge, since PPIs are mediated by surfaces which are much larger and flatter than “classic” binding sites found in enzymatic receptors.1,2,15–19 Among the computational techniques usually applied for drug design/discovery, those based on molecular dynamics (MD) and allowing the estimation of binding energies are becoming more and more exploited.20,21 In particular MMGBSA22,23 has aroused high interest, because it provides an acceptable balance between speed and reliability.24–31 Many studies report a good correlation between MMGBSA predicted binding energies and experimental activities in classic ligand-receptor complexes (i.e. systems where the binding site is well defined and the natural substrates are small molecules/peptides).32–41 Sometimes, a correct reproduction of experimental results has only been obtained by including one or more explicit water molecule during the MMGBSA analysis, although the selection of these waters can be non-trivial.37,38,40–44 The consideration of the crystallographic water molecules mediating H-bonds between the two partner proteins represents the most intuitive approach. However, it has been observed that the crystallographic bridging water can be replaced by a neighboring one during the simulation time, since more than one water can contribute to total water-mediated H-bond occupancy.40 This observation does not conflict with data obtained

ACS Paragon Plus Environment

2

Page 3 of 46

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

from X-ray analyses, since crystallographic residues are identified as mean electron densities and these densities can be determined by many water molecules concurring for a well-defined position.45 Therefore, to fully consider the effect of water on binding energy, all those water molecules that act as H-bond bridges during a MD should be included in the MMGBSA calculation, making the process (selection of waters to be included and preparation of corresponding topological files) long and laborious. To overcome this problem, we recently proposed a MMGBSA variant, hereafter called NwatMMGBSA, based on the inclusion, as part of the receptor, of a fixed number of water molecules (Nwat) which are the closest to a particular residue or to a selection of residues (the ligand or, in case of PPIs, the residues at the interface) in each frame of the MD simulation.40 This method allows to automatically consider in the calculation important water molecules, for example those bridging the ligand and the receptor through H-bonds, even if a water is replaced by a neighboring one during the MD simulation. The application of Nwat-MMGBSA to four different case studies, based on classic ligand-receptor complexes, increased the correlation coefficient between predicted binding energies and experimental activities by about 15% to 78%, depending on the model system.40 Some examples on the ability of MMGBSA to predict relative binding energies of PPI systems can be found in literature,46–48 but the effect of the inclusion of explicit solvent molecules during MMGBSA calculations has been poorly investigated,46 although it is known that water molecules can stabilize PPIs, for example by bridging the protein partners.49–52 Indeed, although protein-protein interfaces are rich with apolar residues, the systematic analysis of known hot spots (i.e. residues accounting for most of the binding energy) revealed that their composition is

ACS Paragon Plus Environment

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 4 of 46

not random53 and that tryptophan, tyrosine, and arginine, which are able to form H-bonds with their side chains, are the most frequently occurring and conserved amino acids.54 In the light of this, also considering our interest in the design of PPI modulators,55–57 we investigated the ability of Nwat-MMGBSA to estimate relative binding energies of a set of 20 heterodimeric PPI complexes having known ∆‫ܩ‬௕௜௡ௗ (Table 1) and one or more water molecules at the interface.58 Table 1. PDB ID, X-ray resolution, and experimental ∆‫ܩ‬௕௜௡ௗ of the selected PPI complexes. PDB ID

Resolution (Å)

Exp. ∆‫ܩ‬௕௜௡ௗ (kcal/mol)

PDB ID

Resolution (Å)

Exp. ∆‫ܩ‬௕௜௡ௗ (kcal/mol)

1ACB59

2.0

-13.05

1ZHI60

2.7

-9.08

1AVX61

1.9

-12.50

2HLE62

2.05

-10.09

1AY763

1.7

-13.23

2HRK64

2.05

-10.98

1BVN65

2.5

-15.06

2OOB66

1.9

-5.66

1EMV67

1.7

-18.58

2OUL68

2.2

-11.96

1FLE69

1.9

-12.28

2SIC70

1.8

-13.84

1GLA71

2.6

-6.76

2SNI72

2.1

-15.96

1KAC73

2.6

-10.68

2UUY74

1.15

-11.26

1R0R75

1.1

-14.17

3BZD76

2.3

-9.57

1YVB77

2.7

-11.17

3SGB78

1.8

-14.51

In detail, we evaluated how the inclusion of 0 to 50 explicit water molecules in the MMGBSA calculations affected the linear correlation between predicted and experimental data, using the coefficient of determination (r2) as a metric. In addition, we assessed the robustness of the method and looked for optimal MD and MMGBSA protocols by testing two recent AMBER force fields, ff14SB79 and ff99SBildn,80 and two implicit solvent models, GB-OBC(II)81 and GB-

ACS Paragon Plus Environment

4

Page 5 of 46

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

Neck2.82 We also compared the correlations obtained by analyzing the fourth and the twelfth ns of the MD trajectory, in order to check the dependency of the method on simulation length. Finally, a “fixed distance” approach to select explicit water molecules was also evaluated and compared to Nwat-MMGBSA. Methods. Structure preparation. A strength of the Nwat-MMGBSA method is that no prior knowledge on structural water molecules is actually needed.40 For this reason, the method can be applied to any kind of 3D complex coming from X-ray, NMR, and homology modelling or docking. Consequently, crystallographic water molecules were removed from the PDB files of the PPI complexes. Successively, the structure preparation and protonate 3D tools of MOE have been used prepare the starting geometries for the MD simulations.83 In particular, the N- and C-termini of those protein chains having more than three missing residues were capped with an acetyl and a methyl-amino group, respectively. All complexes were then protonated at physiological conditions (pH = 7, T = 300 K, salinity = 0.1 M). MD simulations. MD simulations were performed with the pmemd module of Amber14 package,84 using either the ff99SBildn80 or the ff14SB79 force fields. In each complex, the total charge was neutralized by adding an adequate number of Na+/Cl- ions, and the systems were solvated with an octahedral box of TIP3P85 water added up to a distance of 10 Å from the solute. Then, accordingly to the protocol used in previous studies, the systems were relaxed by optimizing the geometry of hydrogens (1000 cycles of steepest descent and 5000 cycles of conjugated gradient), ions and waters (2000 cycles of steepest descent and 5000 cycles of conjugated gradient). The solvent box was equilibrated at 300 K by 100 ps of NVT and 100 ps of NPT simulation using a Langevin thermostat with a collision frequency of 2.0 ps-1. Successively,

ACS Paragon Plus Environment

5

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 6 of 46

a minimization of the geometry of side chains, water and ions with backbone restraints of 25 kcal/mol and a total minimization with backbone restraints of 10 kcal/mol (2500 cycles of steepest descent and 5000 cycles of conjugated gradient) were performed. The systems were then heated up to 300 K in 6 steps of 5 ps each (∆T = 50 K), where backbone restraints were reduced from 10.0 kcal/mol to 5 kcal/mol. Full equilibration was performed in the NVT ensemble (100 ps, backbone restraints = 5.0 kcal/mol) and in the NPT ensemble (1 step of 200 ps, backbone restraints = 5 kcal/mol; 3 steps of 100 ps each, reducing the backbone restraints from 5.0 kcal/mol to 1.0 kcal/mol, and 1 step of 1 ns with 1.0 kcal/mol of backbone restraints). Finally, unrestrained production runs were conducted at 300K for 4 to 12 ns. An electrostatic cutoff of 8.0 Å, a Berendsen barostat, PME for long-range electrostatic interactions, and the SHAKE algorithm were applied to all the calculations. For each run, the root mean squared displacement (RMSD) from the starting coordinates, computed on the backbone atoms, has been used as a measure of convergence (Figure S10 and S11). H-bonds analyses (donor – acceptor distance = 4.0 Å, angle = 150°) were performed with cpptraj, together with grid analyses for the generation of water density maps (50 Å × 50 Å × 50 Å, mesh = 0.5 Å), which were visualized with Chimera.86 Nwat-MMGBSA. MMGBSA analyses were performed with the MMPBSA.py script included in the AmberTools15 package.84 The analyses were conducted on either the 4th or the 12th ns of the production runs by selecting 100 evenly spaced out snapshots, accordingly to our previous experience.27,28,30,40 Either the GB-Neck282 or the GB-OBC(II)81 implicit solvent models were used for the GB calculations (using mbondi3 and mbondi2 radii set, respectively), and the salt molar concentration was set at 0.15 M.

ACS Paragon Plus Environment

6

Page 7 of 46

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

Explicit water molecules that will be included in both the complex and the receptor topologies and trajectories (here referred as Nwat) are then selected by using the cpptraj command “closest”,84 which strips an MD trajectory by all solvent molecules except those that are the closest to the atoms or residues defined in a specific mask, in this case formed by the interface residues belonging to both proteins. For example, supposing that residue 38 and residue 302, belonging to the first and second protein chain, respectively, are identified as interface residues, the cpptraj command “closest 10 :38,:302” will delete from the trajectory every water molecule except the ten that are the closest to the center of mass of residues number 38 and 302. We evaluated Nwat = 0, 10, 20, 30, 40 or 50 and the largest chain of the PDB file, including the coordinates of selected water molecules (Nwat > 0), was considered as the receptor in MMGBSA calculations. Indeed, we previously observed that considering explicit water as part of the ligand (the smallest protein, ins this case) provides comparable results in terms of binding energies, but with a notable increase in standard deviations.40 Cpptraj was also used to select the desired number of snapshots and to generate the updated topology files needed for MMGBSA calculations, while interface residues were identified in a robust and reproducible way using the “InterfaceResidues” Pymol script (see Discussion).87,88 All the steps needed to perform Nwat-MMGBSA calculations were automatized through shell scripts that are provided as Supporting Information (SI), together with the MD protocol input files starting complex coordinate files in PDB format and a tutorial file to reproduce the calculations here reported. Results and discussion. The Nwat-MMGBSA protocol has been principally developed as a method to efficiently rank compounds (small molecules, peptides or protein) in terms of binding affinity, taking into

ACS Paragon Plus Environment

7

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 46

account important water molecules in a fully automatic and unbiased manner. In this method, the entropic term is neglected since the computational cost of normal modes analyses of large systems can be prohibitive when tens or more complexes need to be evaluated, and the benefit given by its inclusion is controversial.33,34,41,89 Similarly, although the Nwat method can be applied to MMPBSA calculations as well,40,90 PB calculations are extremely time-consuming and strongly dependent on MD simulation protocol.21 Moreover, when relative binding free energies are computed, the benefit of PB, compared to GB, is not always significant.26,33 Thus, in this article, only the MMGBSA approach will be discussed. Since the entropic term is neglected and the calculated energies include the contribution given by the selected explicit water molecules (when Nwat > 0), only relative ∆‫ܩ‬௕௜௡ௗ are meaningful. For this reason, the correlation between predicted and experimental energy (and not other metrics such as the mean absolute error) has been considered to assess the quality of prediction. In particular, in this work we report our results on the ability of Nwat-MMGBSA to rank different protein complexes in terms of binding affinity, and we provide optimized variables for its application on PPIs. Selection of the interface residues. When the Nwat-MMGBSA method is applied to classical ligand-receptor complexes, the selection of explicit water molecules is done on the basis of their closeness to the ligand atoms.40 Conversely, when Nwat-MMGBSA is applied to PPIs, the water molecules are selected by evaluating their closeness to all the interfacial residues, or to a subset of them. The selection of the interface residue might thus be critical for the robustness and reproducibility of Nwat-MMGBSA calculations, since it influences the position of the water molecules that will be explicitly considered. The “InterfaceResidues” pymol script provides an unbiased method of selection,88 since it defines as interface residues those amino acids whose

ACS Paragon Plus Environment

8

Page 9 of 46

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

difference in solvent accessible surface area (dASA) from the complex to a single chain is greater than a given cutoff (when the cutoff is set to “0”, all residues are selected).88 Since an optimal cutoff needed to be determined, we evaluated values of 1.0, 0.75, 0.50, 0.25, 0.10 and 0.05 Å2. Moreover, we ought to find out if the optimal residue mask to be passed to the cpptraj “closest” command should include all the amino acids at the interface or just the polar ones (ALL or POLAR setup, respectively), whose surroundings are more likely to be populated by water. We observed that the correlation between experimental and predicted binding energies was not significantly affected either by the cutoff value or by the nature of interfacial residues (Figure 1). Comparable results were also obtained when the Nwat-MMGBSA calculations were performed on MD simulations conducted with the ff99SBildn force field (Figure S1).

ACS Paragon Plus Environment

9

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 46

Figure 1. Variation of r2 on dependency of the number of explicit water molecules included in the calculation (Nwat). Nwat-MMGBSA analyses were performed on the 4th ns of the ff14SB MD simulation, using GB-Neck2 as the implicit solvent model. The poor dependency of r2 on both the chosen cutoff and the nature of residues defined in the mask can be explained by observing that most of the water molecules included in the NwatMMGBSA analysis are common to all methods (water molecules with red oxygen, Figure 2). Indeed, only few water molecules (water molecules with blue or yellow oxygen, Figure 2) were not matching when different selection conditions were used. Cutoff values of 0.25 Å2, 0.10 Å2 and 0.05 Å2 led to the same selection obtained with a cutoff of 0.50 Å2. Analogously, cutoff values of 0.75 Å2 and 1.0 Å2 gave the same result. In the light of this, the following discussion will be limited to the 0.50, POLAR setup, while results obtained with the 0.75, ALL setup are given in the SI.

Figure 2. Snapshot of 2OUL complex (ff14SB, GB-Neck2, Nwat = 20) obtained by defining in the cpptraj input mask: A) All or polar residues with a dASA cutoff = 0.50 Å2, b) Polar residues with dASA cutoff = 0.50 or 0.75 Å2. Water molecules with red oxygen are those identified by

ACS Paragon Plus Environment

10

Page 11 of 46

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

both methods compared in A or B; water molecules with yellow or blue oxygen are those peculiar of one selection method, as explained in the legend. Role of Nwat. The r2 values obtained from the Nwat-MMGBSA analyses showed that the inclusion of 20 – 30 interfacial water molecules during the calculation positively affected the correlation between predicted and experimental binding energies (Figure 1). Indeed, r2 increased by about 30% if compared to that of the standard approach (Nwat = 0). However, the positive effect was partially lost when Nwat > 30, probably because of “background noise” given by energy fluctuations due to the presence of more solvent molecules than needed, as previously observed for some ligand-receptor complexes.40 Interestingly, although protein-protein interfaces are much larger than the classic ligandreceptor interaction surfaces,91 the optimal number of explicit waters (i.e. the value of Nwat providing the highest correlation between Nwat-MMGBSA and experimental binding energies) found here is comparable to that previously observed for ligand-receptor complexes.40 This might be explained by considering the high frequency of apolar residues at the protein-protein interfaces, although the residues accounting for most of the binding energy (i.e. the hot spots) are polar residues that are more likely to be involved in water-bridged H-bonds, such as tryptophan, tyrosine and arginine.54 Therefore, although the protein-protein interfaces are usually in the range of 1500-3000 Å2,15–18 only few water molecules actually affect the PPI stability.

ACS Paragon Plus Environment

11

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 12 of 46

Figure 3. Correlation between experimental ∆‫ܩ‬௕௜௡ௗ and predicted binding energies obtained from the Nwat-MMGBSA analysis of the 4th ns of the ff14SB MD trajectory with GB-Neck2 and Nwat = 0 (A) and 30 (B). Although all the complexes benefited from the inclusion of 20 – 30 interfacial water molecules, the Nwat-MMGBSA approach was particularly advantageous for some of them, such as 3BZD, 1KAC, 1YVB and 2OUL. Indeed, when these complexes are removed from the analysis, the correlation between predicted and experimental binding energies raises to 0.70 even with Nwat = 0. Conversely, Nwat-MMGBSA only weakly affected the results obtained for other complexes, such as 2UUY, 2SNI and 1BVN (Figure 3). In detail, the binding energies of 3BZD and, to minor extent, 1KAC were underestimated when Nwat = 0, suggesting a direct participation of water in the interaction between the two partner proteins. Indeed, a high number of frames (> 200) with water-mediated H-bonds between the two partners were detected by Hbond analysis for these complexes (Table 2). Conversely, a complete absence, or a low number, of water-mediated H-bonds between the two interacting proteins were found for 1BVN and 2UUY or 2SNI, respectively. Moreover, the occupancies of the water-mediated H-bonds of these two latter complexes were lower compared to those observed for 3BZD and 1KAC (Table 2). In addition, the different importance of water in facilitating the PPI of the considered complexes

ACS Paragon Plus Environment

12

Page 13 of 46

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

can also be highlighted by the water density plots obtained from the grid analysis of the MD trajectories of 3BZD and 2UUY, chosen as examples (Figure 4).

ACS Paragon Plus Environment

13

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 14 of 46

Table 2. Water Mediated H-bonds Bridging the Two Partner Proteins of 3BZD, 1KAC, 2UUY and 2SNI.a Occupancies are Reported as the Total Number of Frames (NF) where the Bridging is Observed. 3BZD

1KAC

Bridging residues

NF

Bridging residues

NF

55:GLU 56:LYS 129:THR

754

13:ASP 221:GLU

685

55:GLU 313:ASP

686

85:GLN 219:ASP

675

70:GLN 71:GLU 169:ASN

683

13:ASP 295:ASN

637

55:GLU 315:PHE

675

27:LYS 233:ASP

516

51:ALA 199:TYR 200:VAL 317:GLN 642

50:GLY 248:TYR

390

68:PRO 135:TYR

594

48:VAL 235:VAL

264

71:GLU 169:ASN

522

23:GLU 245:TYR

261

71:GLU 168:LYS

453

16:PRO 240:SER

225

71:GLU 168:LYS 169:ASN

394

23:GLU 49:LYS 248:TYR

216

55:GLU 129:THR

382

63:TYR 284:ASN

324

55:GLU 129:THR 313:ASP

282

46:HID 129:THR

238

2UUY

2SNI

Bridging residues

NF

Bridging residues

NF

131:TYR 256:GLU

366

218:ASN 317:TYR

700

40:HIE 238:CYX

306

99:ASP 323:ARG

548

194:GLY 237:LEU

273

62:ASN 314:THR

362

78:SER 237:LEU

210

99:ASP 306:VAL

282

155:ASN 316:GLU

236

a

th

H-bond analysis were conducted on the 4 ns of the ff14SB MD simulation. 3BZD: chain A, residues 1-110; chain B, residues 111-345. 1KAC: chain A, residues 1-185; chain B, residues 186-309. 2UUY: chain A: residues 1-223; chain B, residues 224-275. 2SNI chain A, residues 1275; chain B, residues 276-339.

ACS Paragon Plus Environment

14

Page 15 of 46

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 4. Water density plots obtained by grid analysis of the 3BZD (left) and 2UUY (right) complexes visualized with Chimera (step 1, level 15). Conversely, 1YVB and 2OUL binding affinities were overestimated with Nwat = 0, whereas with Nwat = 20 – 30 their predicted binding energies correlated well with experimental ∆‫ܩ‬௕௜௡ௗ . When the binding affinity is overestimated during standard MMGBSA calculations, the positive effect on the correlation with experiments given by the inclusion of explicit water molecules cannot be directly ascribed to the contribution of water-mediated interactions between the two partner proteins. Indeed, if the complex is stabilized by explicit waters, the result is a strengthening of the binding between the two partners. Although, it is possible that the overestimation of binding affinities is caused by an underestimation of the receptor stability. Therefore, the inclusion of explicit waters can lower the estimated energy of the receptor, thus leading to an increase in the calculated binding energy. This is well represented by water density plots obtained for 1YVB, where many high water density areas are found within the receptor (Figure S2), and for 3BZD, where the majority of high water density areas are found at the interface between the two partners (Figure 4).

ACS Paragon Plus Environment

15

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 46

Simulation length. It has been reported that the MD simulation length does not necessarily affect the quality of MMPBSA/GBSA calculations in classic ligand-receptor complexes.26,33,90 To find out if this still holds for PPI complexes, we performed Nwat-MMGBSA calculations by analyzing either the 4th or the 12th ns of the MD trajectories. When considering the ff14SB simulation, we observed that there is not a significant difference (< 0.10) in correlations between experimental binding energies and those predicted in the two above conditions (Figure 5). As expected, the same considerations are true whatever criteria are used for the selection of interfacial residues (Figure S3).

Figure 5. Dependency of r2 from Nwat in Nwat-MMGBSA analyses performed on the 4th (blue bars) and 12th (grey bars) ns of the ff14SB MD simulation using GB-Neck2 as implicit solvent model. The criteria for the selection of the interfacial residues are shown in the legend.

ACS Paragon Plus Environment

16

Page 17 of 46

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 6. Correlation between experimental ∆‫ܩ‬௕௜௡ௗ and predicted binding energies obtained from the Nwat-MMGBSA analysis of the 12th ns of the ff14SB MD simulations with Nwat = 0 (A) and 30 (B). Although these findings are in-line with some results reported in the literature,26 part of the lack of increase in correlation for the longer MD simulation can be ascribed to the behavior of 1AVX, whose binding affinity is overestimated when the analysis is conducted on the 12th ns (Figures 3 and 6 and Tables S6 and S7). Indeed, the removal of this complex rises r2 from 0.76 to 0.84 and from 0.77 to 0.78 for the analysis conducted on the 12th and 4th ns of the MD trajectories, respectively. As expected, the overestimation of 1AVX binding affinity is neither improved nor worsened by the inclusion of explicit water, since the water density plot obtained for this system showed the presence of a relatively low number of high water density regions at the 1AVX interface (Figure 7). However, H-bond analyses conducted on the 4th and on the 12th ns showed in the latter trajectory the presence of four additional H-bonds, not mediated by water (Table 3 and Figure 7). Since these H-bonds were not found in the crystal structure,61 they might be responsible of the observed overestimation.

ACS Paragon Plus Environment

17

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 46

Table 3. H-bondsa between the two Partner Proteins of the 1AVX Complex at the 4th and 12th ns of MD Simulation.b 4th ns

12th ns

Acceptor

Donor

Occ% Acceptor

Donor

Occ%

H23@O(BB)

R288@NH2(SC)

91.9

S192@O(BB)

R286@NH(BB)

90.0

S192@O(BB)

R286@NH(BB)

91.3

G78@O(BB)

Y285@OH(SC)

88.4

E235@O(SC)

S129@OH(SC)

91.2

H23@O(BB)

R288@NH2(SC)

88.2

F24@O(BB)

R288@NH(BB)

79.8

R286@O(BB)

G175@NH(BB)

71.3

S172@O(SC)

R286@NH2(SC)

76.4

D171@O(SC)

R286@NH2(SC)

69.2

R286@O(BB)

S177@NH(BB)

68.4

F24@O(BB)

R288@NH(BB)

67.0

G196@O(BB)

R286@NH2(BB)

64.0

E235@O(SC)

Y131@OH(SC)

64.7

P284@O(BB)

G194@NH(BB)

50.2

P284@O(BB)

G194@NH(BB)

59.5

D171@O(SC)

R286@NH2(SC)

48.4

G196@O(BB)

R286@NH2(SC)

49.3

R286@O(BB)

G175@NH(BB)

41.6

S172@O(SC)

R286@NH2(SC)

49.2

D365@O(BB)

N79@NH(SC)

31.6

E235@O(SC)

S129@OH(SC)

44.6

R286@O(BB)

S177@NH(BB)

29.2

E235@O(SC)

Y131@OH(SC)

22.7

E235@O(SC)

S129@OH(SC)

20.9

C359@O(BB)

N79@NH(SC)

20.9

a

Only the H-bonds with percentage occupancy (occ%) > 20% are reported. b BB = backbone, SC = side chain.

ACS Paragon Plus Environment

18

Page 19 of 46

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 7. Water density plots obtained by grid analysis of the 1AVX complex visualized with Chimera (step 1, level 15). The residues involved in the additional H-bonds found at the 12th ns of MD simulations are labelled. Analogous results were obtained by analyzing the ff99SBildn MD simulations using the same MMGBSA conditions, although slightly better r2 values were reached when the 12th ns was analyzed, compared to the 4th ns (Figures S4 and S5). Force field. Concerning the effect of the force field used for the generation of solvated MD trajectories, the best correlations between predicted and experimental binding energies were obtained with the ff14SB force field, whatever value of Nwat was used in Nwat-MMGBSA calculations (Figure 8 and Figure S6).

ACS Paragon Plus Environment

19

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 46

Figure 8. Dependency of r2 from Nwat in Nwat-MMGBSA analyses performed on the 4th ns of the ff14SB (blue bars) and ff99SBildn (grey bars) MD simulations using GB-Neck2 as implicit solvent model. The criteria for the selection of the interfacial residues are showed in the legend.

Figure 9. Correlation between experimental ∆‫ܩ‬௕௜௡ௗ and predicted binding energies obtained from the Nwat-MMGBSA analysis of the 4th ns of the ff99SBildn MD simulations with GBNeck2 and Nwat = 0 (A) or 30 (B). As previously observed, the analysis of the 12th ns of the ff99SBildn MD trajectories gave slightly better results as those performed on the 4th ns (Figure S5), levelling the results to those

ACS Paragon Plus Environment

20

Page 21 of 46

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

obtained with ff14SB. This suggests that the simulations conducted with ff99SBildn need to be longer than those using ff14SB for a proper evaluation of the binding energy. However, the same trend in the correlation index can be observed when applying NwatMMGBSA on trajectories obtained with either ff14SB or ff99SBildn force fields, with an increase in the r2 value of about 30% when passing from Nwat = 0 to Nwat = 30 (Figure 8). While confirming the robustness of the Nwat-MMGBSA approach against different setups, these results also suggest that Nwat = 30 can be considered as a default when calculating relative binding energies of PPI complexes, similarly to what has previously been observed when calculating relative binding energies of drug-like molecules.40 Implicit solvent model for MMGBSA. Two different solvent models, namely the GB-OBC(II)81 and the most recent GB-Neck2,82 were also evaluated within Nwat-MMGBSA calculations on both ff14SB and ff99SBildn trajectories (Figure 10 and S7). For both simulations, the overall best results were obtained with the most recently developed GB-Neck2 model, where the benefit of including explicit water molecules in MMGBSA calculations was also evident. Conversely, this benefit is somehow dampened when the GBOBC(II) model is used. Moreover, even if the highest correlation between predicted and experimental energy was obtained for Nwat = 10, a relatively high r2 was found when the GBOBC(II) model and Nwat = 0 were used to process the ff14SB, but not the ff99SBildn, MD trajectories.

ACS Paragon Plus Environment

21

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 46

Figure 10. Dependency of r2 from Nwat in Nwat-MMGBSA analyses performed on the 4th ns of MD simulations (top: ff14SB; bottom: ff99SBildn) using GB-Neck2 (blue bars) or GB-OBC(II) (grey bars) as the implicit solvent model.

Figure 11. Correlation between experimental ∆‫ܩ‬௕௜௡ௗ and predicted binding energies obtained from the Nwat-MMGBSA analysis of the 4th ns of the ff14SB MD simulations with Nwat = 0 (A) and 30 (B) and by using the GB-OBC(II) as implicit solvent model.

ACS Paragon Plus Environment

22

Page 23 of 46

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

Taking into account that several recent articles found GB-Neck2 better than GB-OBC(II) in reproducing experimental results,82,92,93 we hypothesize that the relatively high r2 observed for Nwat = 0 when GB-OBC(II) was used to process ff14SB trajectories is due to a fortuitous cancellation of errors. Moreover, the higher detrimental effect on the correlation index obtained from the MMGBSA analyses with GB-OBC(II) and Nwat ≥ 20 performed on both the ff14SB (Figures 10, 11 and S7) and ff99SBildn (Figures 10 and S9) MD simulations might be due to the stabilizing effect of explicit waters combined with the overestimation of salt bridges by GBOBC(II). This assumption can be supported by observing the behavior of the 1GLA complex, which is the main actor responsible for the fall of correlation observed with GB-OBC(II) and Nwat ≥ 20 (Figure 11). This complex, together with a particularly high number of salt bridges found at the interface (Table S2), showed a particularly extended area of high water density (Figure S8). The combination of these two features can then lead to an overestimation of its binding affinity when the GB-OBC(II) implicit solvent model was used (Figures 11 and S9). Conversely, when GB-Neck2 is used, the predicted binding energy of this complex is correctly computed or only slightly overestimated, depending on the adopted conditions (Figures 3 and 6). Therefore, we think that the combination of the most recent ff14SB force field and GB-Neck2 model should be the proper choice when using either the classic MMGBSA or the NwatMMGBSA variant. Selecting a fixed number of waters or waters at a fixed distance. Nwat-MMGBSA is based on the selection of a number of explicit water molecules which is kept constant for all snapshots selected for the MMGBSA analysis and for all the evaluated complexes. Some authors reported the use of a distance criterion to select explicit water molecules to be included in MM-PB/GBSA calculations.42 In this case, the distance from the ligand atoms (or from interface residues, for

ACS Paragon Plus Environment

23

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 24 of 46

PPIs) is fixed, while the number of selected waters is variable among the snapshots of each complex. This method might be better than ours in selecting only relevant water molecules, since an arbitrary number of water molecules is not defined a priori. However, the variable number of water molecule selected for each snapshot might be questionable if systematic errors that depend on the explicit water molecules are introduced in the calculations, as hypothesized by other authors for similar questions.94 Moreover, here and elsewhere,40 we observed that considering a number of water molecules greater than the optimal only marginally influence the correlation of computed and experimental binding energies. On the basis of the above assumptions, we decided to compare the Nwat-MMGBSA method with MMGBSA calculations conducted by including, as part of the receptor, explicit water molecules selected among those that are at a fixed distance from any heavy atom belonging to the protein chain considered as the receptor and that considered as the ligand (hereafter referred as Dwat-MMGBSA). Accordingly to reference 42, we initially considered a threshold distance of 3.5 Å, but other values were also evaluated (2.0, 2.5, 3.0 and 4.0 Å). The detailed protocol followed to perform Dwat-MMGBSA calculations is reported in the SI, together with the minimum and maximum number of water molecules selected for each complex (Table S19, SI) energy results (Table S20, SI), and correlations to experiments (Figure S12, SI). To evaluate the reproducibility of both Nwat-MMGBSA and Dwat-MMGBSA, the MD simulations, using the optimal ff14SB force field, were repeated twice (see Table S21 and Figure S13 for Nwat-MMGBSA and Table S22 and Figure S14 for DwatMMGBSA). From the results reported in Tables S20-S22 and Figures S12-S14 appears that, in this specific case, the Nwat-MMGBSA provided a better correlation with experiments and a better reproducibility (top r2 = 0.77 and 0.72 for the two independent run, both obtained for Nwat = 30) than the Dwat-MMGBSA method (top r2 = 0.30 and 0.50 for the two independent run; the

ACS Paragon Plus Environment

24

Page 25 of 46

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

former and the latter are obtained by using a distance of 2.0 Å and 3.0 Å, respectively). The low performance here obtained by Dwat-MMGBSA calculations might be ascribed to the very large difference between the minimum and maximum number of water molecules selected for each complex (see Table S19), but this discussion goes beyond the purpose of this article. Conclusions. Nowadays, the modulation of PPIs represents one of the main goals of drug design/discovery, thanks to the biological relevance of these systems and their link to several disease states.1–4,14 Therefore, a fast, but reliable, method able to rank different PPI complexes (i.e. the natural complex against complexes with potential inhibitors, which might be small organic molecules, peptides or proteins) can be useful in computational drug design. Although other approaches have been developed to specifically predict PPIs binding energies,95 a method which could be widely applied on either PPIs or classic receptor-ligand complexes might represent an interesting alternative. In this framework, we tested a MMGBSA variant, called Nwat-MMGBSA, that was previously shown to improve the correlation between predicted binding energies and experimental affinities for a selection of ligand-receptor complexes,40 on a set of PPI complexes having known binding affinities. The Nwat-MMGBSA approach consists in including, as part of the receptor, a constant number of water molecules which are the closest to the ligand or, in the case of PPI, to interfacial residues, in each of the MD trajectory frames used in the calculation. In this way, even if a given water molecule (for example a water bridging the two partners of a complex) is replaced by a neighboring one during the MD simulation, its presence will always be correctly considered in the MMGBSA analysis. The selection of water molecules can thus be done without any particular knowledge of the protein physicochemical properties. Thus, NwatMMGBSA calculations can be easily done in a fully automatic way, even on multiple structures,

ACS Paragon Plus Environment

25

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 26 of 46

by using standard shell scripts (available from the authors upon request) and tools included in the AmberTools distribution, that can be obtained at no charge for non-profit applications.84 Therefore, in the present study we tested the ability of the Nwat-MMGBSA method to improve the correlation between predicted and experimental binding energies of 20 heterodimeric PPI complexes with known ∆‫ܩ‬௕௜௡ௗ . To optimize the Nwat-MMGBSA protocol and to test its robustness, we investigated how the force field (ff14SB79 and ff99SBildn80), the MD simulation length (4 and 12 ns) and the implicit solvent models used in MMGBSA calculations (GB-Neck282 and GB-OBC(II)81) affected the correlation between predicted and experimental data. We found that the inclusion of 30 interfacial water molecules in Nwat-MMGBSA calculations can increase the r2 correlation index by up to 30% as compared to the standard approach. We also found that the best results are obtained by using the most recent ff14SB force field and the newest implicit solvent model GBNeck2. Moreover, we observed that going beyond the 4 ns of MD simulation did not provide any advantage, although it should be considered that all starting geometries were derived from X-ray structures, and longer simulations might be needed when starting from homology models and/or from complexes obtained by docking. Finally, a remark concerning the computational cost can also be given. Nwat-MMGBSA calculations can be conducted with a negligible increase in computation time respect to classic MMGBSA calculations. In the example reported here, the difference between MMGBSA (Nwat = 0) and Nwat-MMGBSA calculations in terms of total computation time per complex was between 3.4 and 28.0 seconds (see Table S23, SI). On the basis of results here reported, as well as on those derived by on-going studies, we believe that Nwat-MMGBSA might represent a valid and robust alternative to predict relative binding energies both in classic ligand-receptor and in PPI complexes. Nevertheless, further protocol

ACS Paragon Plus Environment

26

Page 27 of 46

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

optimization is currently being attempted to make Nwat-MMGBSA faster and even more appealing for drug design/discovery purposes. ASSOCIATED CONTENT Supporting Information. Details of the PPI complexes, details of MMGBSA calculation results, selected water density maps, correlation plots, selected Hbond analyses. Scripts, starting coordinates in PDB format and input files necessary to reproduce the calculations here described. This material is available free of charge via the Internet at http://pubs.acs.org. AUTHOR INFORMATION Corresponding Author * e-mail [email protected] Tel. +39.02.503.14480; Fax. +39.02.503.1447 Author Contributions The manuscript was written through contributions of both authors. Both authors have given approval to the final version of the manuscript. Notes The authors declare no competing financial interest. ACKNOWLEDGMENT We acknowledge Università degli Studi di Milano for maintaining the convention with the CINECA high performance computer service and the Italian Ministry of Education, University and Research (MIUR) for the PhD scholarship that supported Dr. Irene Maffucci. We also thanks Luisa Collini and Giacomo Cislaghi for their excellent work during their MS degree thesis.

ACS Paragon Plus Environment

27

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 28 of 46

ABBREVIATIONS MD, molecular dynamics; MMGBSA, molecular mechanics Generalized Born surface area; PPI, protein-protein interaction; H-bond, hydrogen bond; dASA, difference in solvent accessible area. REFERENCES (1)

Fischer, P. M. Protein-Protein Interactions in Drug Discovery. Drug Des. Rev. - Online 2005, 2, 179–207.

(2)

Metz, A.; Ciglia, E.; Gohlke, H. Modulating Protein-Protein Interactions: From Structural Determinants of Binding to Druggability Prediction to Application. Curr. Pharm. Des. 2012, 18, 4630–4647.

(3)

Zinzalla, G.; Thurston, D. E. Targeting Protein–protein Interactions for Therapeutic Intervention: A Challenge for the Future. Future Med. Chem. 2009, 1, 65–93.

(4)

Cory, S.; Adams, J. M. The Bcl2 Family: Regulators of the Cellular Life-or-Death Switch. Nat. Rev. Cancer 2002, 2, 647–656.

(5)

Villunger, A.; Scott, C.; Bouillet, P.; Strasser, A. Essential Role for the BH3-Only Protein Bim but Redundant Roles for Bax, Bcl-2, and Bcl-W in the Control of Granulocyte Survival. Blood 2003, 101, 2393–2400.

(6)

Xenarios, I.; Eisenberg, D. Protein Interaction Databases. Curr. Opin. Biotechnol. 2001, 12, 334–339.

(7)

Wells, J. A.; McClendon, C. L. Reaching for High-Hanging Fruit in Drug Discovery at Protein-Protein Interfaces. Nature 2007, 450, 1001–1009.

(8)

Blazer, L. L.; Neubig, R. R. Small Molecule Protein-Protein Interaction Inhibitors as CNS Therapeutic Agents: Current Progress and Future Hurdles. Neuropsychopharmacology 2009, 34, 126–141.

(9)

Ryan, D. P.; Matthews, J. M. Protein-Protein Interactions in Human Disease. Curr. Opin. Struct. Biol. 2005, 15, 441–446.

(10)

Arkin, M. R.; Whitty, A. The Road Less Traveled: Modulating Signal Transduction Enzymes by Inhibiting Their Protein-Protein Interactions. Curr. Opin. Chem. Biol. 2009, 13, 284–290.

(11)

Chène, P. Drugs Targeting Protein-Protein Interactions. ChemMedChem 2006, 1, 400– 411.

(12)

Gerrard, J. A.; Hutton, C. A.; Perugini, M. A. Inhibiting Protein-Protein Interactions as an Emerging Paradigm for Drug Discovery. Mini-Reviews Med. Chem. 2007, 7, 151–157.

(13)

Mullard, A. Protein–protein Interaction Inhibitors Get into the Groove. Nat. Rev. Drug Discov. 2012, 11, 173–175.

(14)

Fry, D. C. Protein-Protein Interactions as Targets for Small Molecule Drug Discovery. Biopolymers 2006, 84, 535–552.

(15)

Horton, N.; Lewis, M. Calculation of the Free Energy of Association for Protein

ACS Paragon Plus Environment

28

Page 29 of 46

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

Complexes. Protein Sci. 2008, 1, 169–181. (16)

Janin, J.; Chothia, C. The Structure of Protein-Protein Recognition Sites. J. Biol. Chem. 1990, 265, 16027–16030.

(17)

Jones, S.; Thornton, J. M. Principles of Protein-Protein Interactions. Proc. Natl. Acad. Sci. U. S. A. 1996, 93, 13–20.

(18)

Lo Conte, L.; Chothia, C.; Janin, J. The Atomic Structure of Protein-Protein Recognition Sites. J. Mol. Biol. 1999, 285, 2177–2198.

(19)

Hopkins, A. L.; Groom, C. R. The Druggable Genome. Nat. Rev. Drug Discov. 2002, 1, 727–370.

(20)

Pearlman, D. A.; Charifson, P. S. Are Free Energy Calculations Useful in Practice? A Comparison with Rapid Scoring Functions for the p38 MAP Kinase Protein System †. J. Med. Chem. 2001, 44, 3417–3423.

(21)

Srivastava, H. K.; Sastry, G. N. Molecular Dynamics Investigation on a Series of HIV Protease Inhibitors: Assessing the Performance of MM-PBSA and MM-GBSA Approaches. J. Chem. Inf. Model. 2012, 52, 3088–3098.

(22)

Massova, I.; Kollman, P. A. Combined Molecular Mechanical and Continuum Solvent Approach (MM-PBSA/GBSA) to Predict Ligand Binding. Perspect. Drug Discov. Des. 2000, 18, 113–135.

(23)

Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; 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, 33, 889–897.

(24)

Ferro, S.; De Luca, L.; Agharbaoui, F. E.; Christ, F.; Debyser, Z.; Gitto, R. Optimization of Rhodanine Scaffold for the Development of Protein-Protein Interaction Inhibitors. Bioorg. Med. Chem. 2015, 23, 3208–3214.

(25)

Xu, D.; Lin, S. L.; Nussinov, R. Protein Binding versus Protein Folding: The Role of Hydrophilic Bridges in Protein Associations. J. Mol. Biol. 1997, 265, 68–84.

(26)

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–8421.

(27)

Ferri, N.; Corsini, A.; Bottino, P.; Clerici, F.; Contini, A. Virtual Screening Approach for the Identification of New Rac1 Inhibitors. J. Med. Chem. 2009, 52, 4087–4090.

(28)

Ferri, N.; Radice, T.; Antonino, M.; Beccalli, E. M.; Tinelli, S.; Zunino, F.; Corsini, A.; Pratesi, G.; Ragg, E. M.; Gelmi, M. L.; Contini, A. Synthesis, Structural, and Biological Evaluation of Bis-Heteroarylmaleimides and Bis-Heterofused Imides. Bioorg. Med. Chem. 2011, 19, 5291–5299.

(29)

Pellegrino, S.; Contini, A.; Clerici, F.; Gori, A.; Nava, D.; Gelmi, M. L. 1H-Azepine-4Amino-4-Carboxylic Acid: A New Α,α-Disubstituted Ornithine Analogue Capable of Inducing Helix Conformations in Short Ala-Aib Pentapeptides. Chem. – A Eur. J. 2012, 18, 8705–8715.

(30)

Contini, A.; Cappelletti, G.; Cartelli, D.; Fontana, G.; Gelmi, M. L. Molecular Dynamics

ACS Paragon Plus Environment

29

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 30 of 46

and Tubulin Polymerization Kinetics Study on 1,14-Heterofused Taxanes: Evidence of Stabilization of the Tubulin Head-to-Tail Dimer-Dimer Interaction. Mol. Biosyst. 2012, 8, 3254–3261. (31)

Pellegrino, S.; Ronda, L.; Annoni, C.; Contini, A.; Erba, E.; Gelmi, M. L.; Piano, R.; Paredi, G.; Mozzarelli, A.; Bettati, S. Molecular Insights into Dimerization Inhibition of cMaf Transcription Factor. Biochim. Biophys. Acta 2014, 1844, 2108–2115.

(32)

Hou, T.; Wang, J.; Li, Y.; Wang, W. Assessing the Performance of the Molecular mechanics/Poisson Boltzmann Surface Area and Molecular Mechanics/generalized Born Surface Area Methods. II. The Accuracy of Ranking Poses Generated from Docking. J. Comput. Chem. 2011, 32, 866–877.

(33)

Hou, T.; Wang, J.; Li, Y.; Wang, W. Assessing the Performance of the MM/PBSA and MM/GBSA Methods. 1. The Accuracy of Binding Free Energy Calculations Based on Molecular Dynamics Simulations. J. Chem. Inf. Model. 2011, 51, 69–82.

(34)

Yang, T.; Wu, J. C.; Yan, C.; Wang, Y.; Luo, R.; Gonzales, M. B.; Dalby, K. N.; Ren, P. Virtual Screening Using Molecular Simulations. Proteins 2011, 79, 1940–1951.

(35)

Genheden, S.; Ryde, U. Comparison of End-Point Continuum-Solvation Methods for the Calculation of Protein-Ligand Binding Free Energies. Proteins 2012, 80, 1326–1342.

(36)

Hayes, J. M.; Skamnaki, V. T.; Archontis, G.; Lamprakis, C.; Sarrou, J.; Bischler, N.; Skaltsounis, A.-L.; Zographos, S. E.; Oikonomakos, N. G. Kinetics, in Silico Docking, Molecular Dynamics, and MM-GBSA Binding Studies on Prototype Indirubins, KT5720, and Staurosporine as Phosphorylase Kinase ATP-Binding Site Inhibitors: The Role of Water Molecules Examined. Proteins 2011, 79, 703–719.

(37)

Freedman, H.; Huynh, L. P.; Le, L.; Cheatham, T. E.; Tuszynski, J. A.; Truong, T. N. Explicitly Solvated Ligand Contribution to Continuum Solvation Models for Binding Free Energies: Selectivity of Theophylline Binding to an RNA Aptamer. J. Phys. Chem. B 2010, 114, 2227–2237.

(38)

Greenidge, P. A.; Kramer, C.; Mozziconacci, J.-C.; Wolf, R. M. MM/GBSA Binding Energy Prediction on the PDBbind Data Set: Successes, Failures, and Directions for Further Improvement. J. Chem. Inf. Model. 2013, 53, 201–209.

(39)

Checa, A.; Ortiz, A. R.; de Pascual-Teresa, B.; Gago, F. Assessment of Solvation Effects on Calculated Binding Affinity Differences: Trypsin Inhibition by Flavonoids as a Model System for Congeneric Series. J. Med. Chem. 1997, 40, 4136–4145.

(40)

Maffucci, I.; Contini, A. Explicit Ligand Hydration Shells Improve the Correlation between MM-PB/GBSA Binding Energies and Experimental Activities. J. Chem. Theory Comput. 2013, 9, 2706–2717.

(41)

Wallnoefer, H. G.; Liedl, K. R.; Fox, T. A Challenging System: Free Energy Prediction for Factor Xa. J. Comput. Chem. 2011, 32, 1743–1752.

(42)

Zhu, Y.-L.; Beroza, P.; Artis, D. R. Including Explicit Water Molecules as Part of the Protein Structure in MM/PBSA Calculations. J. Chem. Inf. Model. 2014, 54, 462–469.

(43)

Mikulskis, P.; Genheden, S.; Ryde, U. Effect of Explicit Water Molecules on LigandBinding Affinities Calculated with the MM/GBSA Approach. J. Mol. Model. 2014, 20, 2273–2283.

ACS Paragon Plus Environment

30

Page 31 of 46

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

(44)

Genheden, S.; Ryde, U. The MM/PBSA and MM/GBSA Methods to Estimate LigandBinding Affinities. Expert Opin. Drug Discov. 2015, 10, 449–461.

(45)

Schiffer, C.; Hermans, J. Promise of Advances in Simulation Methods for Protein Crystallography: Implicit Solvent Models, Time- Averaging Refinement, and Quantum Mechanical Modeling. In Methods in Enzymology; Carter, C. W., J., Sweet, R. M., E., Eds.; New York, 2003; pp 412–461.

(46)

Wong, S.; Amaro, R. E.; McCammon, J. A. MM-PBSA Captures Key Role of Intercalating Water Molecules at a Protein-Protein Interface. J. Chem. Theory Comput. 2009, 5, 422–429.

(47)

Gohlke, H.; Kiel, C.; Case, D. A. Insights into Protein–Protein Binding by Binding Free Energy Calculation and Free Energy Decomposition for the Ras–Raf and Ras–RalGDS Complexes. J. Mol. Biol. 2003, 330, 891–913.

(48)

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–2633.

(49)

Ahmad, M.; Gu, W.; Geyer, T.; Helms, V. Adhesive Water Networks Facilitate Binding of Protein Interfaces. Nat. Commun. 2011, 2, 261.

(50)

Ahmed, M. H.; Habtemariam, M.; Safo, M. K.; Scarsdale, J. N.; Spyrakis, F.; Cozzini, P.; Mozzarelli, A.; Kellogg, G. E. Unintended Consequences? Water Molecules at Biological and Crystallographic Protein–protein Interfaces. Comput. Biol. Chem. 2013, 47, 126–141.

(51)

Ahmed, M. H.; Spyrakis, F.; Cozzini, P.; Tripathi, P. K.; Mozzarelli, A.; Scarsdale, J. N.; Safo, M. A.; Kellogg, G. E. Bound Water at Protein-Protein Interfaces: Partners, Roles and Hydrophobic Bubbles as a Conserved Motif. PLoS One 2011, 6, e24712.

(52)

Janin, J. Wet and Dry Interfaces: The Role of Solvent in Protein–protein and protein– DNA Recognition. Structure 1999, 7, R277–R279.

(53)

Lichtarge, O.; Bourne, H. R.; Cohen, F. E. An Evolutionary Trace Method Defines Binding Surfaces Common to Protein Families. J. Mol. Biol. 1996, 257, 342–358.

(54)

Moreira, I. S.; Fernandes, P. A.; Ramos, M. J. Hot Spots--a Review of the Protein-Protein Interface Determinant Amino-Acid Residues. Proteins 2007, 68, 803–812.

(55)

Ruffoni, A.; Ferri, N.; Bernini, S. K.; Ricci, C.; Corsini, A.; Maffucci, I.; Clerici, F.; Contini, A. 2-Amino-3-(Phenylsulfanyl)norbornane-2-Carboxylate: An Appealing Scaffold for the Design of Rac1–Tiam1 Protein–Protein Interaction Inhibitors. J. Med. Chem. 2014, 57, 2953–2962.

(56)

Maffucci, I.; Clayden, J.; Contini, A. Origin of Helical Screw Sense Selectivity Induced by Chiral Constrained Cα-Tetrasubstituted α-Amino Acids in Aib-Based Peptides. J. Phys. Chem. B 2015, 119, 14003–14013.

(57)

Maffucci, I.; Pellegrino, S.; Clayden, J.; Contini, A. Mechanism of Stabilization of Helix Secondary Structure by Constrained Cα-Tetrasubstituted α-Amino Acids. J. Phys. Chem. B 2015, 119, 1350–1361.

(58)

Vreven, T.; Hwang, H.; Pierce, B. G.; Weng, Z. Prediction of Protein-Protein Binding Free Energies. Protein Sci. 2012, 21, 396–404.

(59)

Frigerio, F.; Coda, A.; Pugliese, L.; Lionetti, C.; Menegatti, E.; Amiconi, G.; Schnebli, H.

ACS Paragon Plus Environment

31

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 32 of 46

P.; Ascenzi, P.; Bolognesi, M. Crystal and Molecular Structure of the Bovine αChymotrypsin-Eglin c Complex at 2.0 Å Resolution. J. Mol. Biol. 1992, 225, 107–123. (60)

Hou, Z.; Bernstein, D. A.; Fox, C. A.; Keck, J. L. Structural Basis of the Sir1-Origin Recognition Complex Interaction in Transcriptional Silencing. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 8489–8494.

(61)

Song, H. K.; Suh, S. W. Kunitz-Type Soybean Trypsin Inhibitor Revisited: Refined Structure of Its Complex with Porcine Trypsin Reveals an Insight into the Interaction between a Homologous Inhibitor from Erythrina Caffra and Tissue-Type Plasminogen Activator. J. Mol. Biol. 1998, 275, 347–363.

(62)

Chrencik, J. E.; Brooun, A.; Kraus, M. L.; Recht, M. I.; Kolatkar, A. R.; Han, G. W.; Seifert, J. M.; Widmer, H.; Auer, M.; Kuhn, P. Structural and Biophysical Characterization of the EphB4-EphrinB2 Protein-Protein Interaction and Receptor Specificity. J. Biol. Chem. 2006, 281, 28185–28192.

(63)

Sevcik, J.; Urbanikova, L.; Dauter, Z.; Wilson, K. S. Recognition of RNase Sa by the Inhibitor Barstar: Structure of the Complex at 1.7 A Resolution. Acta Crystallogr.,Sect.D 1998, 54, 954–963.

(64)

Simader, H.; Hothorn, M.; Kohler, C.; Basquin, J.; Simos, G.; Suck, D. Structural Basis of Yeast Aminoacyl-tRNA Synthetase Complex Formation Revealed by Crystal Structures of Two Binary Sub-Complexes. Nucleic Acids Res. 2006, 34, 3968–3979.

(65)

Wiegand, G.; Epp, O.; Huber, R. The Crystal Structure of Porcine Pancreatic α-Amylase in Complex with the Microbial Inhibitor Tendamistat. J. Mol. Biol. 1995, 247, 99–110.

(66)

Peschard, P.; Kozlov, G.; Lin, T.; Mirza, I. A.; Berghuis, A. M.; Lipkowitz, S.; Park, M.; Gehring, K. Structural Basis for Ubiquitin-Mediated Dimerization and Activation of the Ubiquitin Protein Ligase Cbl-B. Mol. Cell 2007, 27, 474–485.

(67)

Kühlmann, U. C.; Pommer, A. J.; Moore, G. R.; James, R.; Kleanthous, C. Specificity in Protein-Protein Interactions: The Structural Basis for Dual Recognition in Endonuclease Colicin-Immunity Protein Complexes. J. Mol. Biol. 2000, 301, 1163–1178.

(68)

Wang, S. X.; Pandey, K. C.; Scharfstein, J.; Whisstock, J.; Huang, R. K.; Jacobelli, J.; Fletterick, R. J.; Rosenthal, P. J.; Abrahamson, M.; Brinen, L. S.; Rossi, A.; Sali, A.; McKerrow, J. H. The Structure of Chagasin in Complex with a Cysteine Protease Clarifies the Binding Mode and Evolution of an Inhibitor Family. Structure 2007, 15, 535–543.

(69)

Tsunemi, M.; Matsuura, Y.; Sakakibara, S.; Katsube, Y. Crystal Structure of an ElastaseSpecific Inhibitor Elafin Complexed with Porcine Pancreatic Elastase Determined at 1.9 Å Resolution. Biochemistry 1996, 35, 11570–11576.

(70)

Takeuchi, Y.; Satow, Y.; Nakamura, K. T.; Mitsui, Y. Refined Crystal Structure of the Complex of Subtilisin BPN’ and Streptomyces Subtilisin Inhibitor at 1.8 Å Resolution. J. Mol. Biol. 1991, 221, 309–325.

(71)

Hurley, J. H.; Faber, H. R.; Worthylake, D.; Meadow, N. D.; Roseman, S.; Pettigrew, D. W.; Remington, S. J. Structure of the Regulatory Complex of Escherichia Coli IIIGlc with Glycerol Kinase. Science 1993, 259, 673–677.

(72)

McPhalen, C. A.; James, M. N. Structural Comparison of Two Serine Proteinase-Protein Inhibitor Complexes: Eglin-c-Subtilisin Carlsberg and CI-2-Subtilisin Novo. Biochemistry

ACS Paragon Plus Environment

32

Page 33 of 46

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

1988, 27, 6582–6598. (73)

Bewley, M. C.; Springer, K.; Zhang, Y. B.; Freimuth, P.; Flanagan, J. M. Structural Analysis of the Mechanism of Adenovirus Binding to Its Human Cellular Receptor, CAR. Science 1999, 286, 1579–1583.

(74)

Paesen, G. C.; Siebold, C.; Harlos, K.; Peacey, M. F.; Nuttall, P. A.; Stuart, D. I. A Tick Protein with a Modified Kunitz Fold Inhibits Human Tryptase. J. Mol. Biol. 2007, 368, 1172–1186.

(75)

Horn, J. R.; Ramaswamy, S.; Murphy, K. P. Structure and Energetics of Protein-Protein Interactions: The Role of Conformational Heterogeneity in OMTKY3 Binding to Serine Proteases. J. Mol. Biol. 2003, 331, 497–508.

(76)

Cho, S.; Swaminathan, C. P.; Kerzic, M. C.; Guan, R.; Yang, J.; Kieke, M. C.; Andersen, P. S.; Krantz, D. M.; Mariuzza, R. A.; Eric, S. J. Manipulating the Coupled Folding and Binding Process Drives Affinity Maturation in a Protein-Protein Complex. To be Publ.

(77)

Wang, S. X.; Pandey, K. C.; Somoza, J. R.; Sijwali, P. S.; Kortemme, T.; Brinen, L. S.; Fletterick, R. J.; Rosenthal, P. J.; McKerrow, J. H. Structural Basis for Unique Mechanisms of Folding and Hemoglobin Binding by a Malarial Protease. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 11503–11508.

(78)

Read, R. J.; Fujinaga, M.; Sielecki, A. R.; James, M. N. G. Structure of the Complex of Streptomyces Griseus Protease B and the Third Domain of the Turkey Ovomucoid Inhibitor at 1.8-Å Resolution. Biochem 1983, 22, 4420–4433.

(79)

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.

(80)

Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-Chain Torsion Potentials for the Amber ff99SB Protein Force Field. Proteins Struct. Funct. Bioinforma. 2010, 78, 1950–1958.

(81)

Onufriev, A.; Bashford, D.; Case, D. A. Exploring Protein Native States and Large-Scale Conformational Changes with a Modified Generalized Born Model. Proteins Struct. Funct. Bioinforma. 2004, 55, 383–394.

(82)

Nguyen, H.; Roe, D. R.; Simmerling, C. Improved Generalized Born Solvent Model Parameters for Protein Simulations. J. Chem. Theory Comput. 2013, 9, 2020–2034.

(83)

Molecular Operating Environment (MOE), 2013.08 ed.; Chemical Computing Group Inc.: Montreal, 2013.

(84)

Case, D. A.; Babin, V.; Berryman, J. T.; Betz, R. M.; Cai, Q.; Cerutti, D. S.; T.E. Cheatham, I.; Darden, T. A.; Duke, R. E.; Gohlke, H.; Goetz, A. W.; Gusarov, S.; Homeyer, N.; Janowski, P.; Kaus, J.; Kolossváry, I.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Luchko, T.; Luo, R.; Madej, B.; Merz, K. M.; Paesani, F.; Roe, D. R.; Roitberg, A.; Sagui, C.; Salomon-Ferrer, R.; Seabra, G.; Simmerling, C. L.; Smith, W.; Swails, J.; Walker, R. C.; Wang, J.; Wolf, R. M.; Wu, X.; Kollman, P. A. AMBER 14; University of California, San Francisco, 2014.

(85)

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.

ACS Paragon Plus Environment

33

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 34 of 46

1983, 79, 926–935. (86)

Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimera--a Visualization System for Exploratory Research and Analysis. J. Comput. Chem. 2004, 25, 1605–1612.

(87)

Schrödinger, L. The PyMOL Molecular Graphics System, Version 1.8. 2015.

(88)

Vertrees, J. InterfaceResidues http://www.pymolwiki.org/index.php/InterfaceResidues.

(89)

Weis, A.; Katebzadeh, K.; Söderhjelm, P.; Nilsson, I.; Ryde, U. Ligand Affinities Predicted with the MM/PBSA Method: Dependence on the Simulation Method and the Force Field. J. Med. Chem. 2006, 49, 6596–6606.

(90)

Maffucci, I.; Contini, A. Tuning the Solvation Term in the MM-PBSA/GBSA Binding Affinity Predictions. In Frontiers in Computational Chemistry; Ul-Haq, J., Madura, J. D., Eds.; Bentham Science Publishers, 2015; Vol. 1, pp 82–120.

(91)

Cheng, A. C.; Coleman, R. G.; Smyth, K. T.; Cao, Q.; Soulard, P.; Caffrey, D. R.; Salzberg, A. C.; Huang, E. S. Structure-Based Maximal Affinity Model Predicts SmallMolecule Druggability. Nat. Biotechnol. 2007, 25, 71–75.

(92)

Nguyen, H.; Pérez, A.; Bermeo, S.; Simmerling, C. Refinement of Generalized Born Implicit Solvation Parameters for Nucleic Acids and Their Complexes with Proteins. J. Chem. Theory Comput. 2015, 11, 3714–3728.

(93)

Maffucci, I.; Contini, A. An Updated Test of AMBER Force Fields and Implicit Solvent Models in Predicting the Secondary Structure of Helical, β-Hairpin and Intrinsically Disordered Peptides. J. Chem. Theory Comput. 2016, 12, 714–727.

(94)

da Silva, E. F.; Svendsen, H. F.; Merz, K. M. Explicitly Representing the Solvation Shell in Continuum Solvent Calculations. J. Phys. Chem. A 2009, 113, 6404–6409.

(95)

Vangone, A.; Bonvin, A. M. Contacts-Based Prediction of Binding Affinity in ProteinProtein Complexes. Elife 2015, 4, e07454.

-

PyMOLWiki

ACS Paragon Plus Environment

34

Page 35 of 46

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

TOC Graphics

ACS Paragon Plus Environment

35

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

Variation of r2 on dependency of the number of explicit water molecules included in the calculation (Nwat). Nwat-MMGBSA analyses were performed on the 4th ns of the ff14SB MD simulation, using GB-Neck2 as the implicit solvent model. Figure 1 122x71mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 36 of 46

Page 37 of 46

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

Snapshot of 2OUL complex (ff14SB, GB-Neck2, Nwat = 20) obtained by defining in the cpptraj input mask: A) All or polar residues with a dASA cutoff = 0.50 Å2, b) Polar residues with dASA cutoff = 0.50 or 0.75 Å2. Water molecules with red oxygen are those identified by both methods compared in A or B; water molecules with yellow or blue oxygen are those peculiar of one selection method, as explained in the legend. Figure 2 74x30mm (300 x 300 DPI)

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 between experimental ∆G_bind and predicted binding energies obtained from the Nwat-MMGBSA analysis of the 4th ns of the ff14SB MD trajectory with GB-Neck2 and Nwat = 0 (A) and 30 (B). Figure 3 60x21mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 38 of 46

Page 39 of 46

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

Water density plots obtained by grid analysis of the 3BZD (left) and 2UUY (right) complexes visualized with Chimera (step 1, level 15). Figure 4 69x27mm (300 x 300 DPI)

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

Dependency of r2 from Nwat in Nwat-MMGBSA analyses performed on the 4th (blue bars) and 12th (grey bars) ns of the ff14SB MD simulation using GB-Neck2 as implicit solvent model. The criteria for the selection of the interfacial residues are shown in the legend. Figure 5 74x43mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 40 of 46

Page 41 of 46

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

Correlation between experimental ∆G_bind and predicted binding energies obtained from the Nwat-MMGBSA analysis of the 12th ns of the ff14SB MD simulations with Nwat = 0 (A) and 30 (B). Figure 6 59x20mm (300 x 300 DPI)

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

Water density plots obtained by grid analysis of the 1AVX complex visualized with Chimera (step 1, level 15). The residues involved in the additional H-bonds found at the 12th ns of MD simulations are labelled. Figure 7 49x27mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 42 of 46

Page 43 of 46

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

Dependency of r2 from Nwat in Nwat-MMGBSA analyses performed on the 4th ns of the ff14SB (blue bars) and ff99SBildn (grey bars) MD simulations using GB-Neck2 as implicit solvent model. The criteria for the selection of the interfacial residues are showed in the legend. Figure 8 66x35mm (300 x 300 DPI)

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 between experimental ∆G_bind and predicted binding energies obtained from the Nwat-MMGBSA analysis of the 4th ns of the ff99SBildn MD simulations with GB-Neck2 and Nwat = 0 (A) or 30 (B). Figure 9 60x20mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 44 of 46

Page 45 of 46

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

Dependency of r2 from Nwat in Nwat-MMGBSA analyses performed on the 4th ns of MD simulations (top: ff14SB; bottom: ff99SBildn) using GB-Neck2 (blue bars) or GB-OBC(II) (grey bars) as the implicit solvent model. Figure 10 95x103mm (300 x 300 DPI)

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 between experimental ∆G_bind and predicted binding energies obtained from the Nwat-MMGBSA analysis of the 4th ns of the ff14SB MD simulations with Nwat = 0 (A) and 30 (B) and by using the GBOBC(II) as implicit solvent model. Figure 11 59x20mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 46 of 46