Comparative Analysis of Electrostatic Models for Ligand Docking

The precise modeling of molecular interactions remains an important goal among molecular modeling techniques. Some of the challenges in the field include the precise definition of a Hamiltonian for biomolecular systems, together with precise parameters derived from Molecular Mechanics Force Fields, for example. The problem is even more challenging when interaction energies from different species are computed, such as the interaction energy involving a ligand and a protein, given that small differences must be computed from large energies. Here we evaluated the effects of the electrostatic model for ligand binding energy evaluation in the context of ligand docking. For this purpose, a classical Coulomb potential with distance-dependent dielectrics was compared with a Poisson-Boltzmann (PB) model for electrostatic potential computation, based on DelPhi calculations. We found that, although the electrostatic energies were highly correlated for the Coulomb and PB models, the ligand pose and the enrichment of actual ligands against decoy compounds, were improved when binding energies were computed using PB as compared to the Coulomb model. We observed that the electrostatic energies computed with the Coulomb model were, on average, ten times larger than the energies computed with the PB model, suggesting a strong overestimation of the polar interactions in the Coulomb model. We also found that a slightly smoothed Lennard-Jones potential combined with the PB model resulted in a good compromise between ligand sampling and energetic scoring.

Introduction

The quantitative description of molecular interactions, at an atomic level, remains an important challenge even in current days of Petascale computing. Some of the difficulties found in this field include: (i) the energetic description of biomolecular systems; (ii) the fact that binding energies are small differences taken from large energies, resulting in large uncertainties; and (iii) the limited sampling for some calculations. Taken together, these obstacles are exactly the challenge of scoring solutions in the docking problem (Halperin et al., 2002).

The second problem, due to the small differences taken from bigger numbers, can be alleviated with accurate calculations and appropriate sampling. In the context of single point calculations, such as in ligand docking, this challenge remains as an important issue and is handled in some applications with a posterior analysis of ligand candidates using molecular dynamics (MD) or Monte Carlo (MC) simulations to generate an ensemble of thermally accessible configurations of the system and binding energy calculations. In this context, the MM-GBSA or MM-PBSA approaches became very popular (Graves et al., 2008; Zhang et al., 2014; Genheden and Ryde, 2015).

The energetic description of a biomolecular system is tackled in many docking approaches using molecular mechanics force fields (Luty et al., 1995; Schulz-Gasch and Stahl, 2004; Cisneros et al., 2014), where the intermolecular interaction energies are typically computed as a sum of polar interactions, modeled as a Coulomb potential, and van der Waals interactions, modeled through a Lennard-Jones potential (Wang and Lin, 2013). Additional terms can be added to model the influence of the solvent, for example (Muniz and Nascimento, 2017).

Modeling polar interactions using a Coulomb potential introduces some potentially important issues. First, polarization is not considered. Although this effect might be important, a quantum description of the system would be required for appropriate treatment of the dynamics in the electron density within the active site, increasing the computational costs of the calculation. Polarization could also be taken into account by the use of polarizable force fields. However, the computational cost associated with these calculations limits their use in the context of the docking of large compound databases (Illingworth et al., 2008). Second, the dielectric medium of a protein might not be exactly a constant medium, since the protein surface faces the solvent while its core might be closer to a highly hydrophobic medium. So, a representation of the electrostatic potential (and energies) as a function of a varying continuum dielectrics might be necessary, such as the treatment given by the Poisson-Boltzmann (PB) equation (Honig and Nicholls, 1995; Oron et al., 2003; Li et al., 2012, 2013).

Interestingly, Luty et al. observed that, for 20 poses of benzamidine within 8 Å of trypsin binding site, the electrostatic interaction energy computed with PB and using a simple Coulomb model assuming ε = f(r), i.e., the dielectric constant ε is a linear function of the interatomic distance r, showed a high correlation (r² = 0.96) (Luty et al., 1995). In contrary, Gilson and Honig observed that this simple distance-dependent dielectric model overestimates electrostatic interactions (also observed by Luty et al.) and concluded that this model does not seem to be a realistic way of treating polar interactions in biomolecular systems (Gilson and Honig, 1988).

In late '80s Honig et al. developed the DelPhi program (Gilson et al., 1988), that numerically solves the PB equation for macromolecular structures, of any shape, given atomic coordinates, atomic van der Waals parameters, and atomic charges. The calculation of electrostatic potentials within the current versions of Delphi (Li et al., 2012, 2013; Jia et al., 2017) is fast, taking a few seconds in typical workstation computers for a small size protein. Although it might not be fast enough to be used in MD simulations, it is very competitive for docking studies, where the receptor is kept as rigid, in many strategies, and the interaction potentials can be pre-computed in grids and stored for the actual docking calculations (Meng et al., 1992; Luty et al., 1995). Worth of note, the PB calculation is under constant improvement. Recently, Li et al. showed that Gaussian-based smoothed dielectric function could better reproduce the assignment of PKa's for protein residues (Li et al., 2013). The same approach was also applied to the ion distribution (Jia et al., 2017).

Here, we compared the results of docking enrichments and pose reproduction within the same algorithm when using Coulomb electrostatics with a distance-dependent dielectric model (i.e., ε = r) and using a PB electrostatic potential pre-computed using DelPhi (Li et al., 2012). Concurrently, we evaluated the influence of Lennard-Jones soft-core potential on docking efficacy with both PB and Coulomb models. We found that the PB electrostatic model resulted in modest improvement in pose reproduction and enrichment. However, when this model was combined with a smoothed van der Waals potential, an important improvement of pose reproduction and enrichment was observed, suggesting that fine-tuning of these terms is necessary.

Materials and Methods

Docking Calculations

For all docking calculations reported in this work, the software LiBELa (dos Santos Muniz and Nascimento, 2015) was used. LiBELa (Ligand Binding Energy Landscape) uses a combination of ligand- and receptor-based strategies for docking. For this purpose, the algorithm requires a reference ligand, that indicates the initial binding mode. The docking procedure starts with a superposition of the search ligand onto the reference ligand by using a ligand-based approach, as previously described (Vaz de Lima and Nascimento, 2013). Briefly, LiBELa describes the volume of each i ligand atom as a Gaussian function (Vaz de Lima and Nascimento, 2013):

$\begin{array}{l}\rho (r)=p_i{exp}^{\left\{\left[-\pi {\left(\frac{3p_i}{4\pi {\sigma }_i^3}\right)}^{\frac{2}{3}}\right]{(r-r_i)}^2\right\}}\end{array}$

where p_i is the Gaussian amplitude, defined as $2\sqrt{2}$, r_i is atomic coordinate for atom i and σ_i is the van der Waals radius for the same atom. Using this Gaussian-based description of shape, an overlay volume for two molecules, A and B, can be defined as (Vaz de Lima and Nascimento, 2013):

$\begin{array}{l}{V}_{AB}=w{\displaystyle \sum _{i\in A}}{\displaystyle \sum _{j\in B}}{\displaystyle \int }dr{\rho }_i(r){\rho }_j(r)\end{array}$

Similar terms are added to V_AB to compute for the superposition of atoms with positive charge and negative charge with weights defined by w (here, set to 1.0 for both terms). Thus, by maximizing the overlay volume V_AB in Cartesian space, an initial optimized placement of the search ligand is obtained. Afterward, this initial binding mode is re-optimized to find a minimum in the binding energy using a global optimization algorithm. In this step, a typical force field-based definition of binding energy is used as the objective function:

$\begin{array}{l}{E}_{bind}={\displaystyle \sum _i^{rec}}{\displaystyle \sum _j^{lig}}\frac{q_i{q}_j}{ϵr_{ij}}+\frac{A_{ij}}{r_{ij}^{12}}-\frac{B_{ij}}{r_{ij}^6}\end{array}$

Where q is atomic charge, r_ij is the interatomic distance between atoms i and j and A_ij and B_ij are the Lennard-Jones parameters for the atom pair ij, computed by the geometric mean approximation. Here, $A_i=2{\delta }_i{(2r_0)}^{12}$ and $B_i=2{\delta }_i{(2r_0)}^6$, where r₀ is the atomic radius and δ is the well depth parameter necessary for computing van der Waals interactions according to AMBER force field. Both parameters are taken from AMBER FF14SB force field (Maier et al., 2015). A final similarity index can be computed using a Hodgkin's similarity index (Hodgkin and Richards, 1987) defined as:

$\begin{array}{l}SI=\frac{2V_{AB}}{V_{AA}+V_{BB}}\end{array}$

To speed up the calculations, the receptor interaction potential is pre-computed and stored in grids. In this point, LiBELA can compute a typical Coulomb electrostatic potential:

$\begin{array}{l}\varphi (r)={\displaystyle \sum _i^{rec}}\frac{q_i}{\epsilon r_i}\end{array}$

where the dielectric constant was set to the interatomic distance r, i.e., ϵ = r_ij (Luty et al., 1995). Alternatively, LiBELa can parse a DelPhi electrostatic map with the electrostatic potential ϕ_DelPhi instead and compute binding energies using this stronger electrostatic model. For the calculations shown in this work, DelPhi 8.4 was used (Li et al., 2012, 2013). Typically, a computation box of 30 × 30 × 30 Å with a spacing of 0.4 Å (gsize 75 and scale 2.5 Å, in Delphi parameters), with interior dielectrics of 2.06 and exterior dielectrics of 78.5, and salt concentration set to 145 mM. The same grid spacing was used in calculation employing the Coulomb model.

We also tested the effect of a smoothed Lennard-Jones potential by applying the same strategy as suggested by Verkhivker et al. (1999). Here, the binding energy is evaluated as (dos Santos Muniz and Nascimento, 2015; Muniz and Nascimento, 2017):

$\begin{array}{l}{E}_{bind}={\displaystyle \sum _j^{lig}}{q}_j\varphi \left(r\right){\displaystyle \sum _i^{rec}}{\displaystyle \sum _j^{lig}}\frac{A_{ij}}{{\left(r_{ij}^6+{\delta }_{VDW}^6\right)}^2}-\frac{B_{ij}}{r_{ij}^6+{\delta }_{VDW}^6}\end{array}$

The smoothing term δ_VDW was systematically varied in the interval 0.5 to 2.5 Å with a step of 0.5 Å to evaluate the effect of the Lennard-Joned soft-core potential in pose reproduction and enrichment when combined with a Coulomb electrostatic potential (ϕ_Coulomb) or a PB electrostatic potential (ϕ_Delphi).

Docking Pose Reproduction

Self-Docking Test

For docking pose reproduction, we used three data sets. The dataset SB2012 (Mukherjee et al., 2010) includes 1,043 crystal structures of protein-ligand complexes, distributed as SYBYL MOL2 files. Here and all over this text, the “receptor” is defined as a protein where an organic small molecular, the “ligand.,” binds. In these files, the atomic charges are already defined using AMBER forcefield for receptor and AM1-BCC (Jakalian et al., 2000, 2002) for ligands. The dataset files were used as provided, with no further optimizations or modifications of atomic coordinates. Here, a docking calculation was set using each ligand-receptor pair, using the own ligand as the reference ligand in LiBELa.

Cross-Docking Test

For a cross-docking experiment, the Astex dataset was used (Verdonk et al., 2008). In this dataset, 58 structures with analogous complexes are provided. From this dataset, 54 targets were used together with 860 ligands in total. The targets (receptors) were prepared using DockPrep tool as available in UCSF Chimera (Pettersen et al., 2004) using AMBER FF14SB atomic charges. For the ligands, AM1-BCC atomic charges were attributed using ANTECHAMBER (Wang et al., 2006) and SYBYL atom types were assigned using the same tool. In this experiment, each ligand was docked on different (non-native) crystal structures of its own target. Afterward, the root mean square deviation (RMSD) was computed using the native ligand structure as a reference.

Enrichment Tests

In order to evaluate the ability to enrich actual ligands against decoys, i.e., compounds with similar physicochemical properties but not expected to bind to a given target, the DUD38 subset of DUD-E database, which contains 38 targets from the original DUD dataset (Huang et al., 2006), but rebuild with the same protocol as used in DUD-Enhanced (DUD-E) (Mysinger et al., 2012). This subset includes the PDB files for the receptors and over 630,000 compounds, among binders and decoys, with an average decoy-to-ligand ratio of 33. The compounds were used as provided (as SYBYL MOL2 files) with atomic charges defined following the default ZINC protocol (Irwin and Shoichet, 2005; Irwin et al., 2012). The receptor files were prepared using the DockPrep tool available in UCSF Chimera (Pettersen et al., 2004). In this tool, atomic charges are attributed to receptor atoms following AMBER FF14SB parameters. Finally, the prepared receptor is saved as a SYBYL MOL2 file type.

The target-specific ligands and decoys were docked to each target using LiBELa default parameters and using either a Coulomb electrostatic model or a pre-computed Delphi electrostatic potential. The Delphi calculations were carried out in two steps. In the first step, a calculation is set where the protein represents 50% of the calculation box. In a second step, a focused calculation was carried out using a grid of 0.4 Å for a 30 × 30 × 30 Å calculation box centered in the center of mass of the reference ligand. The energies computed after docking calculations were used to rank the docked molecules and ROC curves were computed with locally developed python scripts. The enrichment was quantified using the Adjusted LogAUC metric (Mysinger and Shoichet, 2010). This metric is similar to the well-known AUC but is computed for a semi-logarithmic plot of the ROC curve spanning three decades in the horizontal axis. The computed area is then corrected to remove the area expected for a random enrichment (14.5%).

Results

The calculations of the electrostatic potentials with DelPhi are very fast, typically taking <5 s in an Intel Xeon E5645 (2.40 GHz) processor running in a single thread. This is much faster than the calculation of the interaction potential grids in LiBELa, which took about 5.4 min averaging over the 38 targets of the DUD38 dataset. The computational efficiency of the electrostatic calculations with DelPhi makes it tempting to use this more robust model in docking calculations. However, what is the actual role of PB-based calculation on protein-ligand interactions in the context of ligand docking? In order to address this question, we set up a comparative analysis of ligand binding poses and ligand enrichments using the Coulomb electrostatic model or PB as the electrostatic model to assess the effect of the model in sampling and enrichment, respectively.

Effect on Sampling

The simplest experiment one can think of to evaluate the sampling effect on ligand docking is to assess the ability of a model to reproduce ligand poses from crystal structures. In this context, we compared the RMSDs observed for LiBELa using either Coulomb or PB as the electrostatic model.

After the redocking of 1,029 ligands on their respective receptors, the RMSD for all atoms, including hydrogen atoms, was computed in comparison with the original (experimental) structures. Averaging over the entire dataset an RMSD of 1.215 Å was observed for the Coulomb model, while for the PB model an average RMSD of 1.129 Å was achieved. The median RMSD for these models were 0.535 and 0.598 Å with a standard deviation slightly increased for Coulomb as compared to PB (1.689, compared to 1.419 Å). For both models the fraction of targets with RMSD values found below the typical cutoff value of 3.0 Å was close to 90%, as indicated in Table 1.

TABLE 1

Table 1. Summary of the self-docking experiment using the SB2012 dataset.

When a smoothed Lennard-Jones potential was combined with the electrostatic models under evaluation in this work, we found very interesting differences. For a small smoothing parameter δ = 0.5 Å, the differences between the electrostatic models are small, similarly to what is observed in the AMBER Lennard-Jones model. However, as δ becomes larger, the differences between the Coulomb model and the PB model become more evident. When δ is set to 2.0 Å, the average RMSD found for the Coulomb model was 5.643 Å (median 5.566 Å), while the average RMSD for the PB model was 1.185 Å, with median in 0.611 Å. So, it appears that the combination of the PB model with a soft-core VDW potential still leads to good results in pose reproduction while the Coulomb model rapidly seems to dominate the binding energy resulting in meaningless ligand poses.

Another interesting observation comes from the comparison between the polar term in the interaction energies. An analysis for 1,029 protein-ligand complexes reveals a good correlation between the electrostatic interaction energies computed using a Coulomb model with distance-dependent dielectrics, i.e., ε = r_ij, and electrostatic interaction energies computed using the Poisson-Boltzmann model. As shown in Figure 1, there is a good correlation between the computed energy terms (r = 0.7 for N = 1,029), as also observed previously by Luty et al. (1995). Additionally, one can observe that the electrostatic interaction energies computed by the Coulomb model are about 10 times more favorable, on average than those computed using the PB model, indicating a typical overestimation of the interaction energies in this model. In the context of ligand docking, this overestimation may result in binding modes that are biased toward a few polar contacts that are too favorable as compared to the overall fitting of the ligand and receptor binding pockets.

FIGURE 1

Figure 1. Correlation of the electrostatic interaction energies computed with Coulomb model (horizontal axis) and PB (vertical axis). The line shows a linear regression of the obtained data (N = 1,029) and the regression coefficients are shown in the figure.

A more stringent test is to assess the ability of the model to reproduce experimentally determined poses in a cross-docking experiment, i.e., in a receptor structure different from the one used in the actual docking calculation. In brief, it involves docking ligand L_A in receptor structure R_B and comparing the docking pose to the pose observed when L_A was crystallized bound to receptor R_A. For this task, the Astex non-diverse dataset was used.

The employed dataset includes 603 diverse (non-native) complexes. The results obtained are summarized in Table 2. Again, when the typical AMBER Lennard-Jones potential is used, a slight improvement in the binding poses is observed, with average RMSD going from 4.48 Å, for the Coulomb model, to 3.95 Å in the PB model (median values of 3.27 and 2.94 Å, respectively). When a smoothed Lennard-Jones potential is used, on the other hand, the differences between both electrostatic models increase. For a smoothing parameter δ set to 0.5 Å, the average RMSD decreases from 4.65 in the Coulomb model to 3.90 in the PB model (median values 3.52 and 3.06 Å). And when δ is set to 2.0 Å, the average RMSD decreases from 7.60 to 3.71 Å, with median values decreasing from 7.52 to 2.54 Å (Table 2).

TABLE 2

Table 2. Summary of the cross-docking experiment using the Astex dataset.

Taken together, the results shown here indicate that the PB model for electrostatic computation result in better pose reproduction in the scenario of a typical AMBER FF binding energy calculation and, more significantly, in the scenario of a soft docking, i.e., when the Lennard-Jones potential is smoothed. Given the results obtained, we moved for the evaluation of the changes in the enrichment of actual binders against decoy compounds.

Ligand Enrichment

In order to assess the ability of the electrostatic models to recover actual ligands against decoy compounds, we choose the DUD38 dataset. In this dataset, 38 targets are given with a set of binder compounds and a set of decoy compounds. In this context, a decoy is defined as a compound that has similar physicochemical properties to the binders but is not expected to bind to the receptor. After docking all the binders and decoys, the compounds are ranked by their binding energy and a receiver-operating characteristic (ROC) curve is plotted. Finally, the enrichment is computed using Adjusted LogAUC metric, as previously proposed (Mysinger and Shoichet, 2010). Briefly, this metric describes the area under the curve of the ROC plot with the x-axis in the logarithm scale and spanning three decades. The area computed is corrected by subtracting the area expected for a random enrichment.

The results obtained are summarized in Table 3 and shown in the complete version in the Supplementary Material. From the data shown here, we note that, for the usual Lennard-Jones model used in AMBER force field, i.e., δ_VDW = 0.0, the electrostatic models performed almost similarly in terms of enrichment, with an average enrichment of 4.50 or 4.91 for Coulomb or PB, respectively, with a slight improvement of the enrichment with the PB model. Using the smoothed VDW potential with δ_VDW = 0.5 Å, similar enrichments are observed but with an improvement in the median enrichment for the PB model. Here, the average enrichments were 5.01 and 4.99 with median enrichments of 3.93 and 5.40 for Coulomb and PB models, respectively. Finally increasing the smoothing constant to δ_VDW = 2.0 Å, a maximum in the average/median logAUC is observed for the PB model (5.36 and 6.09 for average and median, respectively), while a marked decrease in the enrichment for the Coulomb model is observed.

TABLE 3

Table 3. Summary of the enrichment experiment using the DUD38 dataset.

For the sake of comparison, the same docking calculations using the DUD38 were set up using the Grid Score model of DOCK 6.7 (Brozell et al., 2012). The average and median logAUC observed for this model was 1.3 and −1.2, respectively (Table 3). Since logAUC corrects for the expected random enrichment, this metric can achieve negative results if results are worse than random. It is important to add that the Grid Score used here has a 6–12 Lennard-Jones potential with a Coulomb electrostatic model that uses a distance-depend dielectric function (ε = r_ij), similar to the model used in LiBELa.

Discussion

The enrichment data shown in Table 3 for the DUD38 dataset strongly suggests that the continuum electrostatic model can lead to important improvements in the ability to recover actual binders and separate them from decoy compounds. On the other hand, as we already noted from the data shown in Figure 1, the Coulomb interaction electrostatic energies are, on average, 10 times more favorable than interaction electrostatic energies computed with PB. Then, it makes sense that the balance between the electrostatic and van der Waals terms should be also fine-tuned. We assessed this balance by introducing a smoothed Lennard-Jones term to model the van der Waals interactions. A good balance seems to be achieved when the smoothing constant δ_VDW was set to 2.0 Å. With this calculation setup, a maximum in the enrichment is observed, without compromising the docking poses, according to the results of ligand enrichment with DUD38, self-docking with the SB2012 dataset (Table 1) and cross-docking with the Astex dataset (Table 2). A complete comparison of the effect of the smoothing parameter δ_VDW is shown in the Supplementary Material, where the ligand pose and ligand enrichment can be compared as a function of the smoothing parameter.

A second effect of the electrostatic treatment given to the docking calculations can be observed in the distribution of the net charges of the top-scored molecules in docking calculations. The analysis of the charge distribution for the target ace, shown as an example in Supplementary Figure 2, reveals that among the top-scored molecules when the Coulomb model was used, almost half of them have net charges −2 or −3 e, indicating a favoring of the non-specific electrostatic interactions to the total docking score. On the other hand, the PB model favors neutral molecules or molecules with net charge −1 e (Supplementary Material). No molecule with net charge −2 or −3 e is observed among the top-scored molecules, suggesting a much more specific scoring of the biomolecular interactions. As a piece of evidence of the correctness of the PB model, an inspection of the distribution of net charges among the actual binders in the DUD dataset for this target shows that 66% of the binders have net charge 0, 30% have charge −1 and 4% have charge −2, indicating that the PB model more closely reflects the molecular interactions observed in experimental conditions.

A recent development in the PB calculations introduced a Gaussian-based approach “to deliver a smooth dielectric function for the entire space domain” (Li et al., 2013). The authors showed that the Gaussian-based function resulted in better assignments of PKa's and also dielectric values for protein interior and protein-water interface in agreement with previous works (Li et al., 2013). Interestingly, a comparison of the ligand enrichment obtained after docking calculations using the Non-Gaussian dielectric model and the Gaussian dielectric model showed a decrease of about 15% in the enrichment of actual ligands against decoys in the DUD38 dataset. This decrease is observed for the AMBER Lennard-Jones model as well as for the smoothed Lennard-Jones model. This still preliminary observation highlights that even the robust PB model can still be improved in the context of ligand docking to result in even more reliable calculations and predictions of protein-ligand interactions.

What is the actual role of PB calculations in protein-ligand recognition? A simpler answer to this question would be a better treatment of intermolecular electrostatic interactions avoiding the overestimation of charge-charge interaction, as observed in the Coulomb model. However, as observed by Basu et al. (2012), there is a close association between electrostatic complementarity and shape complementarity in biological macromolecules. This close association may be also associated with the findings shown here, where tradeoff between a better electrostatic treatment for intermolecular interactions and adjustment of van der Waals interactions results in better recognition model, directly affecting sampling (better binding poses for ligands) as well as scoring (better enrichment of known ligands when compared to decoy compounds).

In conclusion, here we evaluated the effect of scoring docking calculations with a Coulomb model or with a Poisson-Boltzmann model for electrostatic energies. As the major findings, we observe that (i) PB model improves both docking pose and docking scoring capabilities. In the most stringent test for pose reproducibility, the cross-docking test, the best results were found for the combination of PB and softcore Lennard-Jones potential. Here we observed a reduction of the average RMSD from 3.95 to 3.90 Å and then to 3.71 Å, as the softcore constant δ goes from 0.0 to 0.5 Å and then to 2.0 Å. (ii) The best enrichment results for the challenging dataset DUD38 was also observed for the combination of PB and softcore Lennard-Jones potential. Again, taking the median values over the DUD38 dataset, we observed an increase of 15 or 30% if δ_VDW is set to 0.5 or 2.0, respectively. These results confirm previous observations that the Coulomb potential overestimate the electrostatic energy in protein-ligand interaction calculations and also show that a re-parametrization of the Lennard-Jones parameters might be necessary. Finally, given the quick calculations of PB-based electrostatic potentials in modern workstations, we don't see a reason for not using this model in protein-ligand calculations in the context of ligand docking.

Data Availability

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

AN conceived the project. GS conducted the experiments. GS and AN collected and analyzed the data. All authors wrote and approved the manuscript.

Funding

Financial support was provided by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) through grants 2017/18173-0, 2015/26722-8 and 2015/13684-0 and by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), through grants 303165/2018-9, 406936/2017-0 and 168769/2017-4.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors thank the DelPhi developers for making their tool freely available to the scientific community. We also thank Heloisa Muniz and Victor Nogueira for the discussions. We thank the financial support provided by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) through grants 2017/18173-0, 2015/26722-8 and 2015/13684-0 and by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), through grants 303165/2018-9 and 406936/2017-0.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2019.00052/full#supplementary-material

References

Basu, S., Bhattacharyya, D., and Banerjee, R. (2012). Self-complementarity within proteins: bridging the gap between binding and folding. Biophys. J. 102, 2605–2614. doi: 10.1016/j.bpj.2012.04.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Brozell, S. R., Mukherjee, S., Balius, T. E., Roe, D. R., Case, D. A., and Rizzo, R. C. (2012). Evaluation of DOCK 6 as a pose generation and database enrichment tool. J. Comput. Aided. Mol. Des. 26, 749–773. doi: 10.1007/s10822-012-9565-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Cisneros, G. A., Karttunen, M., Ren, P., and Sagui, C. (2014). Classical electrostatics for biomolecular simulations. Chem. Rev. 114, 779–814. doi: 10.1021/cr300461d

PubMed Abstract | CrossRef Full Text | Google Scholar

dos Santos Muniz, H., and Nascimento, A. S. (2015). Ligand- and receptor-based docking with LiBELa. J. Comput. Aided. Mol. Des. 29, 713–723. doi: 10.1007/s10822-015-9856-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Genheden, S., and Ryde, U. (2015). The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin. Drug Discov. 441, 1–13. doi: 10.1517/17460441.2015.1032936

CrossRef Full Text | Google Scholar

Gilson, M. K., and Honig, B. H. (1988). Energetics of charge-charge interactions in proteins. Proteins Struct. Funct. Genet. 3, 32–52. doi: 10.1002/prot.340030104

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilson, M. K., Sharp, K. A., and Honig, B. H. (1988). Calculating the electrostatic potential of molecules in solution: method and error assessment. J. Comput. Chem. 9, 327–335. doi: 10.1002/jcc.540090407

CrossRef Full Text | Google Scholar

Graves, A. P., Shivakumar, D. M., Boyce, S. E., Jacobson, M. P., Case, D. A., and Shoichet, B. K. (2008). Rescoring docking hit lists for model cavity sites: predictions and experimental testing. J. Mol. Biol. 377, 914–934. doi: 10.1016/j.jmb.2008.01.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Halperin, I., Ma, B., Wolfson, H., and Nussinov, R. (2002). Principles of docking: an overview of search algorithms and a guide to scoring functions. Proteins Struct. Funct. Genet. 47, 409–443. doi: 10.1002/prot.10115

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodgkin, E. E., and Richards, W. G. (1987). Molecular similarity based on electrostatic potential and electric-field. Int. J. Quantum Chem. 32, 105–110.

Google Scholar

Honig, B., and Nicholls, A. (1995). Classical electrostatics in biology and chemistry. Science 268, 1144–1149. doi: 10.1126/science.7761829

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, N., Shoichet, B. K., and Irwin, J. J. (2006). Benchmarking sets for molecular docking. J. Med. Chem. 49, 6789–6801. doi: 10.1021/jm0608356

PubMed Abstract | CrossRef Full Text | Google Scholar

Illingworth, C. J., Morris, G. M., Parkes, K. E., Snell, C. R., and Reynolds, C. A. (2008). Assessing the role of polarization in docking. J. Phys. Chem. A 112, 12157–12163. doi: 10.1021/jp710169m

PubMed Abstract | CrossRef Full Text | Google Scholar

Irwin, J. J., and Shoichet, B. K. (2005). ZINC - a free database of commercially available compounds for virtual screening. J. Chem. Inf. Model. 45, 177–182. doi: 10.1021/ci049714+

PubMed Abstract | CrossRef Full Text | Google Scholar

Irwin, J. J., Sterling, T., Mysinger, M. M., Bolstad, E. S., and Coleman, R. G. (2012). ZINC: a free tool to discover chemistry for biology. J. Chem. Inf. Model. 52, 1757–1768. doi: 10.1021/ci3001277

PubMed Abstract | CrossRef Full Text | Google Scholar

Jakalian, A., Bush, B. L., Jack, D. B., and Bayly, C. I. (2000). Fast, efficient generation of high-quality atomic charges. AM1-BCC model: I. Method. J. Comput. Chem. 21, 132–146. doi: 10.1002/(SICI)1096-987X(20000130)21:2<132::AID-JCC5>3.0.CO;2-P

CrossRef Full Text | Google Scholar

Jakalian, A., Jack, D. B., and Bayly, C. I. (2002). Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation. J. Comput. Chem. 23, 1623–1641. doi: 10.1002/jcc.10128

PubMed Abstract | CrossRef Full Text | Google Scholar

Jia, Z., Li, L., Chakravorty, A., and Alexov, E. (2017). Treating ion distribution with Gaussian-based smooth dielectric function in DelPhi. J. Comput. Chem. 38, 1974–1979. doi: 10.1002/jcc.24831

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuhn, H. W. (2010). “The Hungarian method for the assignment problem,” in 50 Years of Integer Programming 1958-2008: From the Early Years to the State-of-the-Art, eds M. Jünger, M. Liebling, D. Naddef, G. L. Nemhauser, W. R. Pulleyblank, G. Reinelt, G. Rinaldi, L. A. Wolsey (Berlin; Heidelberg: Springer), 29–47. doi: 10.1007/978-3-540-68279-0_2

CrossRef Full Text | Google Scholar

Li, L., Li, C., Sarkar, S., Zhang, J., Witham, S., Zhang, Z., et al. (2012). DelPhi: a comprehensive suite for DelPhi software and associated resources. BMC Biophys. 5:9. doi: 10.1186/2046-1682-5-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Li, C., Zhang, Z., and Alexov, E. (2013). On the dielectric “constant” of proteins: smooth dielectric function for macromolecular modeling and its implementation in DelPhi. J. Chem. Theory Comput. 9, 2126–2136. doi: 10.1021/ct400065j

PubMed Abstract | CrossRef Full Text | Google Scholar

Luty, B. A., Wasserman, Z. R., Stouten, P. F. W., Hodge, C. N., Zacharias, M., and McCammon, J. A. (1995). A molecular mechanics/grid method for evaluation of ligand-receptor interactions. J. Comput. Chem. 16, 454–464. doi: 10.1002/jcc.540160409

CrossRef Full Text | Google Scholar

Maier, J. A., Martinez, C., Kasavajhala, K., Wickstrom, L., Hauser, K. E., and Simmerling, C. (2015). ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 11, 3696–3713. doi: 10.1021/acs.jctc.5b00255

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, E. C., Shoichet, B. K., and Kuntz, I. D. (1992). Automated docking with grid-based energy evaluation. J. Comput. Chem. 13, 505–524.

Google Scholar

Mukherjee, S., Balius, T. E., and Rizzo, R. C. (2010). Docking validation resources: protein family and ligand flexibility experiments. J. Chem. Inf. Model. 50, 1986–2000. doi: 10.1021/ci1001982

PubMed Abstract | CrossRef Full Text | Google Scholar

Muniz, H. S., and Nascimento, A. S. (2017). Towards a critical evaluation of an empirical and volume-based solvation function for ligand docking. PLoS ONE 12:e0174336. doi: 10.1371/journal.pone.0174336

PubMed Abstract | CrossRef Full Text | Google Scholar

Mysinger, M. M., Carchia, M., Irwin, J. J., and Shoichet, B. K. (2012). Directory of Useful Decoys, Enhanced (DUD-E): better ligands and decoys for better benchmarking. J. Med. Chem. 55, 6582–6594. doi: 10.1021/jm300687e

PubMed Abstract | CrossRef Full Text | Google Scholar

Mysinger, M. M., and Shoichet, B. K. (2010). Rapid context-dependent ligand desolvation in molecular docking. J. Chem. Inf. Model. 50, 1561–1573. doi: 10.1021/Ci100214a

PubMed Abstract | CrossRef Full Text | Google Scholar

Oron, A., Wolfson, H., Gunasekaran, K., and Nussinov, R. (2003). Using DelPhi to compute electrostatic potentials and assess their contribution to interactions. Curr. Protoc. Bioinformatics Chapter 2, 8.4.1–8.4.12. doi: 10.1002/0471250953.bi0804s02

CrossRef Full Text | Google Scholar

Pettersen, E. F., Goddard, T. D., Huang, C. C., Couch, G. S., Greenblatt, D. M., Meng, E. C., et al. (2004). UCSF chimera - a visualization system for exploratory research and analysis. J. Comput. Chem. 25, 1605–1612. doi: 10.1002/Jcc.20084

PubMed Abstract | CrossRef Full Text | Google Scholar

Schulz-Gasch, T., and Stahl, M. (2004). Scoring functions for protein-ligand interactions: a critical perspective. Drug Discov. Today Technol. 1, 231–239. doi: 10.1016/j.ddtec.2004.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Trott, O., and Olson, A. J. (2009). AutoDock vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J. Comput. Chem. 31, 455–461. doi: 10.1002/jcc.21334

PubMed Abstract | CrossRef Full Text | Google Scholar

Vaz de Lima, L. A., and Nascimento, A. S. (2013). MolShaCS: a free and open source tool for ligand similarity identification based on Gaussian descriptors. Eur. J. Med. Chem. 59, 296–303. doi: 10.1016/j.ejmech.2012.11.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Verdonk, M. L., Mortenson, P. N., Hall, R. J., Hartshorn, M. J., and Murray, C. W. (2008). Protein–ligand docking against non-native protein conformers. J. Chem. Inf. Model. 48, 2214–2225. doi: 10.1021/ci8002254

PubMed Abstract | CrossRef Full Text | Google Scholar

Verkhivker, G. M., Rejto, P. A., Bouzida, D., Arthurs, S., Colson, A. B., Freer, S. T., et al. (1999). Towards understanding the mechanisms of molecular recognition by computer simulations of ligand-protein interactions. J. Mol. Recognit. 12, 371–389. doi: 10.1002/(SICI)1099-1352(199911/12)12:6<371::AID-JMR479>3.0.CO;2-O

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Wang, W., Kollman, P. A., and Case, D. A. (2006). Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graph. Model. 25, 247–260. doi: 10.1016/j.jmgm.2005.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J. C., and Lin, J. H. (2013). Scoring functions for prediction of protein-ligand interactions. Curr. Pharm. Des. 19, 2174–2182. doi: 10.2174/1381612811319120005

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X., Wong, S. E., and Lightstone, F. C. (2014). Toward fully automated high performance computing drug discovery: a massively parallel virtual screening pipeline for docking and molecular mechanics/generalized born surface area rescoring to improve enrichment. J. Chem. Inf. Model. 54, 324–337. doi: 10.1021/ci4005145

PubMed Abstract | CrossRef Full Text | Google Scholar