Molecular docking is a computational method that investigates the intermolecular complexes formed between two or more constituent molecules. It comprises the process of generating a model of a complex based on the known three-dimensional structures of its components, i.e., the target (protein, or nucleic acids) and the ligand (a peptide, a protein, a small organic molecule), free or bound to other species (Rognan, 2013). The docking procedure consists in the search for near-native ligand conformations and orientations (usually referred to as docking poses) with respect to a target protein, where the structure of the latter is known or modeled. Fast approximate mathematical expressions (so called scoring functions) are used to rank the docking poses based on estimates of the goodness of the conformations obtained for each putative binder and of the binding affinity estimate of the two interacting partners. Pioneered during the early 1980s (Kuntz et al., 1982), molecular docking is still a field of intensive research, as it represents a fundamental component in many drug discovery programs (Meng et al., 2011) and a primary tool for the virtual screening of large chemical libraries (Kitchen et al., 2004). The typical system considered in docking calculations includes the ligand, the receptor, and the solvent molecules. Because of the enormous number of degrees of freedom associated with the solvent, it is usually neglected in the calculations, or implicitly accounted for in the scoring functions. Despite some valuable improvements in the accuracy and efficiency of the molecular docking algorithms, there are still considerable drawbacks and limitations to face. Among these, the reliability of the scoring functions is probably one of the aspects deserving more attention, since discriminating native pose and obtaining a fair correlation between docking scores and experimental activity data remain difficult tasks. These limitations are responsible for the occurrence of false-positive and false-negative hits in the ranked lists resulting from the screenings performed with standard docking methods. Over the years, since the pioneering work of Kuntz et al. (1982), several scoring functions have been developed (Gilson and Zhou, 2007; Huang et al., 2010; Sarti et al., 2013), based on several terms. Despite empirical scoring functions are still widely used in drug discovery since they are faster and relatively accurate, first-principle methods for ranking decoy structures and for predicting affinity should be considered the first desirable choice in docking scoring stage.

Hence, it is a general opinion that molecular docking results may benefit from post-processing with more accurate tools, able to provide higher accuracy in energy scoring of the putative docked poses. Among several docking approaches, HADDOCK is one of the few computational docking programs that follow a data-driven strategy, using experimental data (generated either via NMR experiments, mutagenesis, or mass spectrometry) as pivotal information to generate docking poses (Dominguez et al., 2003). Moreover, the program allows the receptor to undergo small conformational changes upon association with the ligand, a feature that has been deemed as crucial in simulating the binding process (Spiliotopoulos and Caflisch, 2014). HADDOCK has been applied successfully to a plethora of biomolecular systems (Dominguez et al., 2003). Its reliability is highlighted by the excellent evaluations in the CAPRI experiments (van Dijk et al., 2005) and by the fact that more than 60 structures solved via HADDOCK docking have been deposited in the Protein Data Bank (Berman et al., 2000). Moreover, continuous efforts are devoted to integrate HADDOCK with experimental methodologies (Hennig et al., 2012) and other computational techniques (Kastritis et al., 2014) in order to improve its built-in scoring function. Among several other scoring approaches, Molecular Mechanics Poisson–Boltzmann Surface Area, MM-PBSA, is routinely used to evaluate the strength of the complex formation between protein and ligands. MM-PBSA represents a good trade-off between calculation efficiency and accuracy in binding energy calculations and it has been profitably exploited in virtual design since it allows a ranking in different docking runs and between different ligands (Graves et al., 2008; Venken et al., 2011; Yang et al., 2011; Barakat et al., 2012; Genheden and Ryde, 2015). Although MM-PBSA is one of the most used approximate methods for the estimate of binding free energies, it also presents weaknesses that should not be overlooked. In particular, a source of error can be represented by the entropy contribution, which is often neglected when relative binding free energies of similar molecules are computed. Furthermore, the quality of results depends on different computational factors, including the conformational sampled space, the force field, internal dielectric constant, and the set of atomic radii (Weis et al., 2006). Additional MM-PBSA limitations in the estimation of binding free energies are for highly polar molecules such as DNA and RNA (Kongsted et al., 2009), buried ligands (Singh and Warshel, 2010), and in presence of explicit water molecules that might contribute to the binding free energy (Homeyer and Gohlke, 2012). Therefore, several attempts have been made to improve accuracy and predictivity of the MM-PBSA method acting on the solvation term, including polar and non-polar terms. Expedients such as using different PB solvers (Feig et al., 2004), tuning the grid mesh (Harris et al., 2013), and/or the internal dielectric constant (Singh and Warshel, 2010; Hou et al., 2011; Genheden and Ryde, 2015), including crystallographic and/or specific water molecules (Treesuwan and Hannongbua, 2009; Liu et al., 2013; Maffucci and Contini, 2013), have allowed to successfully use MM-PBSA in the binding free energy calculations (Wang et al., 2001). However, there is still room for improvement to make MM-PBSA more efficient and reliable in the binding free energy calculations in different respects.

Here we focus our investigation on the field of protein-peptide interactions (PPIs), which are gaining large interest in the biological and pharmaceutical research (Scott et al., 2016). In fact, the inhibition of PPIs is of paramount importance in drug discovery and development. The main problem in the PPIs simulation is that the protein-peptide interface is large, shallow, and involving several contacts characterized by being weak, transient and non-specific. Therefore, not all the PPI interface contributes equally to the strength of the binding between the partners. PPIs are rather mediated by hot spots, small regions that give the largest contribution to the binding.

In the present study, we investigated the effectiveness of MM-PBSA on evaluating protein-peptide docked complexes using HADDOCK software in order to consider the possibility of exploiting this relatively fast approach as an additional scoring function. We show in the results section the performance of MM-PBSA as a scoring function for the 19 systems (Results—Section MM-PBSA As Scoring Function for Protein-Peptide Docking) and as binding affinity predictor for the systems where reliable experimental binding affinity data were available (Results—Section Correlation between Experimental Binding Free Energies and Scores).

MM-PBSA is an end-point method devised to estimate binding free energy (ΔG_comp) as the difference of the free energy of the complex and those of the unbound receptor and peptide (Massova and Kollman, 1999). Normally, it is performed from a set of snapshots obtained from Molecular Dynamics simulation (Hou et al., 2011). This method is significantly less computationally demanding than alternatives such as free energy perturbation (FEP) calculations and therefore it represents a possible alternative to FEP for virtual screening of large chemical libraries. It relies on the use of implicit solvent (for the PB part) and it requires energy calculations only on the endpoint (bound/unbound) states whereas other approaches require energy calculation along a reaction coordinate. MM-PBSA has already been used as a scoring function in the past with various outcomes (Kuhn et al., 2005; Thompson et al., 2008; Zhou et al., 2009; Genheden and Ryde, 2015) but to the best of our knowledge this is the first time it has been used for a set of PPIs obtained from docking calculations. Previous studies have shown that MM-PBSA is efficient to identify the correct binding poses and rank small molecules for a specific target (Thompson et al., 2008; Hou et al., 2011; Zhu et al., 2013). However, there is no systematic evaluation of the performance of MM-PBSA in identifying the correct docking poses in the protein-peptide context. As mentioned, the free energy of binding was calculated as the difference in free energy between the product state and the reactants state, that is, between the energy of the protein-peptide complex and the sum of the energies of the protein and the ligand in their unbound forms.

We investigated 19 protein-peptide complexes for which structural and thermodynamic data (binding free energy values ΔG_bind) were available (Table 1). The final MM-PBSA values are calculated as the sum of two molecular mechanics terms (namely Coulomb and Lennard-Jones), which are calculated by HADDOCK, and two solvation terms, including polar and non-polar solvation contributions, which here were calculated using two different Poisson–Boltzmann equation solvers, APBS (Baker et al., 2001) and DelPhi (Rocchia et al., 2001, 2002) in combination with the NanoShaper program (Decherchi and Rocchia, 2013). We decide to use two different solvers to minimize the MM-PBSA aforementioned weakness and interestingly the binding free energy values from the two solvers were in good agreement in all case studies, indicating that the consistency of the approach. (Figure S1).

In our calculations, we observed that the computed binding free energies were larger than those obtained by experiments (Table 2), an overestimation that has been already observed in other systems. This behavior has been often ascribed to the omission of the entropic contribution, which is an approximation typical of these calculations (Gilson and Zhou, 2007; Spiliotopoulos et al., 2012).

The constituents of each complex (i.e., protein and peptide) were separated and re-docked using HADDOCK. The HADDOCK scores and MM-PBSA binding free energies corresponding to 200 poses were calculated in each system with the HADDOCK's clustering-based approach. Similarly to HADDOCK score calculation, where coulombic interactions are scaled to a fifth, we also dampened MM-PBSA, i.e., the MM-PBSA energies were calculated multiplying both coulombic and polar solvation terms by 0.2. We will therefore compare three different scoring functions, including HADDOCK built-in, MM-PBSA, and dampened MM-PBSA, which we call dMM-PBSA. The two following sections show the performances of each scoring function in discriminating between the correctly and incorrectly docked peptide poses (Section MM-PBSA As Scoring Function for Protein-Peptide Docking) and in correlating with the experimental ΔG_bind through the whole dataset (Section Correlation between Experimental Binding Free Energies and Scores).

In this work, we made use of the MM-PBSA technique in docking scoring and in affinity prediction of protein-peptide complexes. We also compared the results with the HADDOCK built-in scoring function. Overall, HADDOCK and dMM-PBSA, a dampened MM-PBSA version, behaved similarly in ranking the near-native poses in the top 3 clusters, improving over the standard MM-PBSA version. The introduction of weights for the different MM-PBSA terms is not unprecedented in the literature (Zhou et al., 2009), but this approach has never been applied to PPIs. Notably, despite the fact that different experimental conditions (where the main difference regarded the type of buffer and the ionic strength, whereas as both pH and temperature were comparable) and techniques (Table 1) were used to determine the dissociation constants, we observed a good correlation between experimental and computational ΔG of binding using HADDOCK and dMM-PBSA scoring functions, 0.63 and 0.66 respectively. Interestingly, lack of modulation of the solvation MM-PBSA terms resulted in worse correlation between the experimental and simulated figures (r = 0.49).

Our findings were then validated on two additional systems, one with known structure and binding affinity and one for which only the ΔG_bind has been reported. Both HADDOCK and dMM-PBSA methods perform remarkably well in ranking the two additional protein-peptide complexes, and lead to good correlation with experimentally measured ΔG_bind.

In order to assess the presence of possible systematic errors in the binding free energy calculations, we used two different PB solvers, namely APBS (Baker et al., 2001) and DelPhi (Rocchia et al., 2001, 2002). No differences in term of correlation with experimental data were found using the two different solvers, except that in absolute ΔG_comp values (Figure S1) and in the calculation speed (see Section Materials and Methods). APBS provided a ΔG_comp larger than DelPhi. This could arise from different aspects, including the different approaches used to describe the dielectric interface, the approach used to estimate the reaction field energy, and/or how the different solvers treat cavities that are internal to the solute. The accuracy of the PB equation solution has been reported to be sensitive to the grid size, in favor of smaller grid spacing (Sørensen et al., 2015) However, decreasing the grid spacing increases the computational resources needed to perform the calculation, both in terms of physical memory, and the computational time required. We choose a grid size of 0.5 Å for both solvers since it represents a good trade-off between speed and accuracy. With lower grid resolution, APBS would have encountered problems in convergence in ΔG_comp calculation (see Figure 3.2 in Sørensen et al., 2015). Finally, the choice of the interior dielectric (ε_int) value in the PB calculation is not trivial since in the literature its value can be found spanning between 1 and 20. Higher ε_int (4–20) aims to effectively mimic polarization and local rearrangement effects, as well as transient penetration of water molecules into the solute interior. They should in principle be preferred when the PB calculation is performed on individual structures. On the other hand, lower ε_int (2–3) is preferred in order to mainly account for electronic polarization and it is commonly used on ensembles of structures, which explicitly account for conformational flexibility. In light of this, there is still no consensus on the most appropriate ε_int value. We decided to use ε_int = 2 in all cases, while being aware that scaling the polar contribution in dMM-PBSA is similar to considering a dielectric screening of the solute medium, accounting for polarization and rearrangement due to the reaction to the existing fields (Figure S2). This is also in agreement with previous studies where the rank-ordering performance for MM-PBSA improves with increasing dielectric constant (Wang et al., 2013). Simultaneously, this supports the relative importance of the non-polar components and results in better performance. In fact, the van der Waals and BSA terms, which are directly related to the MM-PBSA non-polar terms, correlated well with the experimental data and they provided the main contribution to the final score. These factors might explain the better performance of the dampened MM-PBSA. Finally, we analyzed the calculated data from HADDOCK and MM-PBSA in order to evaluate the reproducibility of the results in terms of scoring and final correlation with the experimental data. First of all, for a given complex, with the corresponding MM energy terms provided by HADDOCK, the only variations in the final results could arise from the PB calculations since it depends on many parameters, including the grid spacing, the atomic radii, and the dielectric interior. For this reason, we used two different PB solvers, APBS and DelPhi, using identical parameters. In Figures S1, S2 it is shown the good agreement in the MM-PBSA and dMM-PBSA calculations with the two PB solvers, indicating that the calculations are quite robust. Second, we calculated the correlations between the experimental and calculated binding free energy using different MM-PBSA values from the cluster_BEST4 and not the average value. The final correlations R were in a range of 0.39–0.49 for MM-PBSA and 0.60–0.67 for dMM-PBSA, indicating that in principle the MM-PBSA could be performed on a single pose with less computational effort. Finally, the calculated standard error for each cluster_BEST4 reported in Table 2 represent the MM-PBSA limit when we try to compare the binding affinities of different complexes, mainly when dealing with docking since the conformational sampling is poor with respect to performing MM-PBSA from molecular dynamics simulations (Spiliotopoulos et al., 2012). Nevertheless, in case of AIRE-PHD1 mutants, the MM-PBSA and dMM-PBSA uncertainty over the NMR structure bundle (20 structures) reduces, leading to an average standard error of 34 and 7 kJ/mol, respectively, indicating that dMM-PBSA is fairly reliable in predicting protein-peptide interface alanine mutations. Moreover, the CAS dMM-PBSA of AIRE-PHD1 error range values are fairly in agreement with the CAS MM-PBSA error range (ca. 4 kJ/mol) shown in Spiliotopoulos et al., in which the sampling was carried out by molecular dynamics simulations (Spiliotopoulos et al., 2012). This result indicates that dMM-PBSA carried out on a small number of structures (e.g., 20) behaves similarly to MM-PBSA from molecular dynamics, in which the sampling is more extended. In fact, the modulation of the polar terms in MM-PBSA is probably taking into account the possible local rearrangement similarly to what is done via a higher internal dielectric, indicating that dMM-PBSA represents a simple and promising approach in evaluating alanine mutations from single structure, either from docking or from X-Ray or from NMR.

We believe that in parallel with the recently developed optimization of the HADDOCK score for PPIs inhibitors (Kastritis et al., 2014), the combination of HADDOCK score and modified MM-PBSA binding free energy might lay the groundwork for novel approaches to study in silico PPIs inhibitors in a quick and automatic fashion. The advantage of using MM-PBSA as scoring function is threefold. First, MM-PBSA is versatile. It provides estimates of the equilibrium averages over the solvent degrees of freedom, permitting the post-processing of solute representative snapshots from docking poses. Since MM-PBSA estimates binding free energy, it represents a valuable alternative since it is in principle transferable between different docking runs and can be used to score both intra-ligand and inter-ligand poses, saving individual validation for each system under study. This makes MM-PBSA more suitable for novel problems with limited experimental data as we demonstrated in the case study of NPHP1-SH3. Second, MM-PBSA can be used to quantify the thermodynamical strength of the putative poses. In particular, our dMM-PBSA is a reliable scoring function in the protein-peptide field showing good a correlation with experimental data. Therefore, the advantage to use dMM-PBSA with respect to MM-PBSA relies on the possibility to modulate the polar terms without rerunning the Poisson–Boltzmann calculations at different internal dielectric values. Finally, MM-PBSA could better allow to disclose atomistic details of protein-peptide binding, supporting the rational design of bioactive compounds. In fact, post-processing task such as CAS approach has been very recently applied to MD simulations for successfully evaluating the importance of key residues in the protein-peptide binding complex (Spiliotopoulos et al., 2012). Therefore, MM-PBSA can be used to calculate the ΔΔG, defined as ΔG_wt − ΔG_mutALA, on the protein-peptide best docked poses in order to identify residues for which mutation to alanine strongly attenuates binding. The latter behavior occurred in our short CAS study of AIRE-PHD1 mutants (Figure S3), by which an acceptable correlation of R = 0.86/0.88 between the experimental and the calculated binding energies of the mutants demonstrated the possibility to use MM-PBSA as a promising tool at low computational cost to evaluate the hot-spots in the PPIs field with respect to the HADDOCK score.

Finally, structural prediction of protein-peptide complexes remains challenging due to two major obstacles: peptides are highly flexible and they often interact weakly with their substrate, underlining their importance in signal transduction or regulation which often relies on transient processes. This leaves flexible docking as one of the few amenable computational techniques to model these complexes (Verkhivker et al., 2000; Hetényi and van der Spoel, 2002; Niv and Weinstein, 2005; Raveh et al., 2010, 2011; Trellet et al., 2013). In our study, we are considering the two interacting partners, protein and peptide, already in the bound conformation, which represents a strong assumption in term of binding mechanism and also the best scenario for a docking calculation, especially in the PPIs field. This undoubtedly increases the success rate in the native pose determination and it reduces the error in the binding affinity calculations since in our MM-PBSA approach we are neglecting most of the solute entropic contribution (i.e., under this assumption the estimate of ΔS = S_Complex − (S_protein − S_peptide) can be poor). Recent work have highlighted the efforts to improve HADDOCK protocol in the field of protein-peptide (Trellet et al., 2013) but still predicting large conformational changes remains a challenge as indicated by several failure to accurately predict cases where the protein undergoes large conformational changes upon binding (Trellet et al., 2013). In this case, even the more accurate and reliable scoring function and binding free energy methods will struggle in discriminating the correct binding mode due to both hard and soft docking failure (Verkhivker et al., 2000).

In conclusion, in contrast to other scoring functions and approximate binding free energy calculation methods such as the linear interaction energy (LIE) method, MM-PBSA contains less empirical parameters and, thus, it is more likely to be useful in determining the relative free energies of binding of quite different compounds and systems for which there is more limited experimental data, although a protein structure of the target is, of course, required.

DS, designed research, ran experiments, analyzed results, wrote the paper; PK designed research, ran experiments, analyzed results; AM designed research, ran experiments, analyzed results; AB designed research, analyzed results; GM designed research, analyzed results, wrote the paper; WR analyzed results, wrote the paper; AS, designed research, ran experiments, analyzed results, wrote the paper.

The research leading to these results has received funding under the Horizon 2020 Program, FET-Open: PROSEQO, Grant Agreement n. [687089]. AS acknowledges funding from AIRC (Associazione Italiana Ricerca sul Cancro) through grant MFAG11899. The development of HADDOCK is supported by grants from the Netherlands Organization for Scientific Research (NWO) (TOP-PUNT grant no. 718.015.001) and by European H2020 e-Infrastructure grants (EGI-Engage, grant no. 654142; INDIGO-DataCloud, grant no. 653549; West-Life grant no. 675858 and BioExcel grant no. 675728). The EGI infrastructure and DIRAC4EGI service with the dedicated support of CESNET-MetaCloud, INFN-PADOVA, NCG-INGRID-PT, RAL-LCG2, TW-NCHC, SURFsara and NIKHEF, and the additional support of the national GRID Initiatives of Belgium, France, Italy, Germany, the Netherlands, Poland, Portugal, Spain, UK, South Africa, Malaysia, Taiwan and the US Open Science Grid are acknowledged.