Abstract
Halogenated ligands are nowadays commonly designed in order to increase their potency against protein targets. Although novel computational methods of evaluating the affinity of such halogenated inhibitors have emerged, they still lack the sufficient accuracy, which is especially noticeable in the case of empirical scoring functions, being the method of choice in the drug design process. Here, we evaluated a series of halogenated inhibitors of phosphodiesterase type 5 with ab initio methods, revealing the physical nature of ligand binding and determining the components of interaction energy that are essential for proper inhibitor ranking. In particular, a nonempirical scoring model combining long-range contributions to the interaction energy provided a significant correlation with experimental binding potency, outperforming a number of commonly used empirical scoring functions. Considering the low computational cost associated with remarkable predictive abilities of the aforementioned model, it could be used for rapid assessment of the ligand affinity in the process of rational design of novel halogenated compounds.
Similar content being viewed by others
Introduction
It has been recognized that halogen compounds (iodine, bromine, chlorine or even, to some extent, fluorine) interact with atoms possessing a lone electron pair (Lewis base, e.g., O, N, S atoms or π-electron systems) through so called halogen bonding, where the halogen atom acts as an acceptor (Lewis acid) [1,2,3]. Halogen bonding is related to the anisotropy of the electron density and the emergence of increased electrostatic potential, i.e., σ-holes [4, 5]. This electrostatically driven, directional, intermolecular interaction [6] has been successfully exploited in the drug development process [7,8,9,10], and is often found in biomolecular complexes [2, 11, 12]. Since commonly used empirical or knowledge-based docking and scoring approaches do not always properly account for the anisotropy of electron charge distribution of halogen residues, the quantum mechanical (QM) methods are usually required, when it comes to the modeling of halogen bonded systems [12,13,14,15]. Importantly, scoring of ligands bearing halogen atoms with empirical scoring functions is often not reliable [14, 16], although a progress in this field has been achieved, and novel empirical scoring functions dedicated to halogenated compounds have appeared [17,18,19].
As the QM-based methods are too computationally demanding to be applied in the drug design process [3, 17], reliable and fast scoring approaches are needed. Our non-empirical \(E_{EL,MTP}^{(10)}+E_{Das}\) function [20], already validated with various protein-ligand complexes [21,22,23], might appear useful. Taking advantage of the long-range interaction energy terms including electrostatic atomic multipole expansion (\(E_{EL,MTP}^{(10)}\)) and approximate dispersion function (EDas), it constitutes a low cost approach that can be used in various biomolecular systems. The damped dispersion expression, EDas, has already been proven to successfully describe non-covalent interactions with atom-atom potentials fitted to reproduce the results of high-level quantum chemical calculations [24, 25]. Importantly, EDas function appears to be an alternative to the commonly applied DFT-based damped dispersion corrections that might not be able to recover dispersion energies properly due to recently demonstrated unphysical damping occurring at small intermonomer distances [26]. However, available EDas parameters have not covered halogen atoms [24, 25], restricting application of \(E_{EL,MTP}^{(10)}+E_{Das}\) function to nonhalogenated compounds.
Novel EDas parameters derived recentlyFootnote 1 have enabled application of the \(E_{EL,MTP}^{(10)}+E_{Das}\) model to scoring of halogen-substituted ligands. Preliminary tests of \(E_{EL,MTP}^{(10)}+E_{Das}\) function augmented by halogen parameters demonstrated its favourable performance for the series of phosphodiesterase type 5 (PDE5) inhibitors. PDE5 enzyme catalyzes hydrolysis of cyclic guanosine monophosphate (cGMP), an intracellular second messenger molecule involved in multiple signaling pathways. Decreased degradation of cGMP due to PDE5 inhibition enhances its effects, leading to, e.g., relaxation of vascular smooth muscle tissue. PDE5 inhibitors (PDE5Is) are used in the treatment of pulmonary hypertension or erectile dysfunction (ED) [27]. ED therapy with PDE5 inhibitors, e.g., sildenafil (Viagra, Pfizer Inc., first generation PDE5Is) or avanafil (Stendra, Metuchen Pharmaceuticals, LLC., second generation PDE5Is) is preferred as the first-line medication [28]. Moreover, the number of various PDE5Is in dietary supplements has increased from 1 in 2003 to 69 in 2016 [28]. Recently, PDE5Is have been proposed as medicines for heart failure [29]. PDE5 has also been linked to tumor development regulation [30]. Due to clinical potential of PDE5 inhibition, development of higher affinity PDE5Is is of special interest. Understanding of the molecular basis of PDE5Is binding could significantly aid the rational design of novel inhibitors.
\(E_{EL,MTP}^{(10)}+E_{Das}\) function complemented by novel halogen EDas parameters has been validated using the set of PDE5 inhibitors, that involved monocyclic pyrimidinones with halogen bonding introduced to strengthen the binding affinity [31]. Here, we report nonempirical investigation of these compounds, providing detailed analysis of the physical basis of their interactions with PDE5. Furthermore, our \(E_{EL,MTP}^{(10)}+E_{Das}\) approach is compared with a number of empirical scoring functions, including those recommended for halogen-bearing ligands (e.g., XBSF implemented in AutoDock VinaXB program [14]). Thermodynamical basis of binding of PDE5 inhibitors studied here have been characterized experimentally by Ren et al. [32]. The experimental results discussed therein were accompanied by theoretical analysis covering only a limited representation of the binding site, leaving a room for further and more complete computational examination, which we have therefore exploited. In particular, calculations for only 3 out of 5 compounds were performed by Ren et al. [32], whereas herein we report computational results for all inhibitors.
Our recent findings [23, 33] suggest that interaction energy calculations should be accompanied by assessment of the solvation effects, since the considerable differences in the solvation free energy might affect the inhibitory activity ranking produced by herein applied theoretical models. As the latter account only for the enthalpic contribution to the binding free energy, their applicability requires consistency of free energy of solvation (ΔGsolv) among the analyzed set of inhibitors. It could be especially important in the case of halogen-bearing ligands [34]. Accordingly, we complement our analysis of halogenated PDE5 inhibitors by the report on the related solvation effects.
Methods
Preparation of protein-inhibitor complexes
A set of 5 PDE5 inhibitors developed by Xu et al. [31] is examined with ab initio and empirical methods to analyze the effect of the halogen substitution on the binding affinity. Crystal structures of PDE5-inhibitor complexes are available via the following PDB codes: 4OEX, 4OEW [32], 3SIE, 3SHY, 3SHZ [31] (resolutions of 2.14, 2.44, 1.93, 2.65, and 2.45 Å, respectively). The inhibitor compounds share a common pyrimidinone scaffold. The only difference between these structures is the substitution at the 5-position of the pyrimidinone ring with a halogen or hydrogen atom (Fig. 1).
Optimal hydrogen bonding was determined with Propka [35,36,37,38], implemented in Maestro. Due to insufficient resolution of the most of the crystal structures (e.g., 3SHY complex with 2.65Å resolution), all protein-inhibitor structures were solvated with TIP3 water model [39] and optimized in Charmm program [40] (version c36b1). Both Charmm General Force Field v. 2b7 [41] and Charmm22 All-Hydrogen Force Field [42,43,44] parameter files were used. Missing parameters for inhibitor structures were generated with CGenFF interface at http://cgenff.paramchem.org [41, 45,46,47] (interface version 1.0.0). All amino acid residues further than 10 Å from each inhibitor were kept frozen throughout 1000 steps of steepest descent followed by conjugate gradient optimization until root mean square gradient of 0.01kcal ⋅mol− 1 ⋅ Å was reached.
To investigate the influence of the substitution on the activity of the inhibitors, the inhibitors were truncated, as the common scaffold was positioned similarly in the case of all complexes. Therefore, only the pyrimidinone substituted at 5-position (marked in red in Fig. 1) was taken into account during the calculations. All amino acid residues within 5 Å of the substituent at pyrimidinone 5-position were selected to serve as PDE5 binding site model. The latter included 11 amino acid residues (Tyr612, Asp764, Leu765, Ala767, Ile768, Gln775, Ile778, Ala779, Val782, Gln817, Phe820; see Fig. 2). Dangling bonds resulting from truncation of the inhibitor structures and cutting the amino acid residues from the protein structure were filled with hydrogen atoms minimized in Schrödinger Maestro [48] program using Opls 2005 force field [49].
Interaction energy calculations
Interaction energy between PDE5 binding site representation and each inhibitor was calculated within Hybrid Variation-Perturbation Theory (HVPT) [50, 51] decomposition scheme as the sum of binding energy values obtained for amino acid residue-inhibitor dimers. Each amino acid residue included in particular dimers was considered separately except for residues linked by a peptide bond, i.e., Asp764 and Leu765, Ala767 and Ile768, as well as Ile778 and Ala779 residues. Interaction energy calculations were carried out with GAMESS program [52] using polarized triple zeta valence basis set of Ahlrichs et al. (def2TZVP) [53, 54] and the corresponding effective core potential (ECP) on iodine atoms. Counterpoise correction was applied [55].
In HVPT, the Møller-Plesset second-order interaction energy (EMP2) is partitioned into the multipole electrostatic (\(E_{EL,MTP}^{(10)}\)), penetration (\(E_{EL,PEN}^{(10)}\)), exchange (\(E_{EX}^{(10)}\)), delocalization (\(E_{DEL}^{(R0)}\)) and correlation (\(E_{CORR}^{(2)}\)) terms:
Interaction energy components introduced above can be grouped into long- and short-range energy terms that vary with the intermolecular distance R as R−n or exp(−aR), respectively. The electrostatic multipole component of the binding energy, \(E_{EL,MTP}^{(10)}\), was assessed with the Cumulative Atomic Multipole Moment (CAMM) expansion [56]. The latter is truncated at R− 4 term and computed using the correlated wave-function. The first-order electrostatic energy, \(E_{EL}^{(10)}\), is equal to the sum of the penetration term, \(E_{EL,PEN}^{(10)}\), and the \(E_{EL,MTP}^{(10)}\) energy. Similarly, the first-order Heitler-London energy, E(10), is obtained by adding the first-order electrostatic energy to the exchange component, \(E_{EX}^{(10)}\). Subsequently, the higher order delocalization energy, \(E_{DEL}^{(R0)}\), which comprises classical induction and charge transfer terms, is the difference between the counterpoise-corrected self-consistent field variational energy, ESCF, and the first-order Heitler-London energy, E(10). Finally, the correlation term, \(E_{CORR}^{(2)}\), accounting primarily for the dispersion, exchange-dispersion, and intramolecular correlation contributions, is obtained as the difference of the second-order Møller-Plesset interaction energy, EMP2, and converged SCF energy, ESCF. The \(E_{EL,MTP}^{(10)}+E_{Das}\) energy values were calculated with a novel revision of the dispersion function, EDas.2
Empirical scoring was performed with several scoring functions. GoldScore, ChemScore, ChemPLP, and Astex Statistical Potential (ASP) functions, implemented in Gold 5.6.3 [57], were employed, with a spherical grid comprising amino acid residues within 10 Å radius from the point of origin, defined as the pyrimidinone 5-position. Further, ChemPLP and PLP scoring functions from Plants [58] docking program were used with a 10 Å radius sphere centered on the ligand. Since ChemPLP function is present in both Gold and Plants, it was assigned a subscript indicating its origin, namely ChemPLPG (Gold) and ChemPLPP (Plants). Despite sharing the same name, both functions are defined differently: ChemPLPG combines ChemScore hydrogen bonding term with linear potentials for modeling the van der Waals and repulsive terms [59], while ChemPLPP is based on the piecewise linear potential (PLP) scoring function [58]. Importantly, ChemPLPG score is given in arbitrary units with higher score indicating more potent inhibitor, whereas ChemPLPP yields the scoring results in kcal ⋅mol− 1 units (lower score indicates more potent inhibitor). GlideSP (standard precision) [60] and GlideXP (extra precision) [61] scoring functions from Schrödinger Glide program (revision 2018-1) [62] were also applied with a 10 Å grid centered on the ligand. XBSF scoring function, implemented in AutoDock VinaXB [14] and designed especially for halogenated ligands, was also employed for the comparison. Scoring was carried out with 10 Å cubic grid centered on the inhibitor.
In all docking programs full structures of protein-inhibitor complexes, optimized as described in Preparation of protein-inhibitor complexes section, were used for scoring of the existing poses of truncated inhibitor structures. In particular, no docking was performed, as the compounds re-docking would affect the results, precluding the comparison of the performance of both empirical and non-empirical scoring. __________ 2Unpublished results (2018)
Assessment of solvation effects
The solvation free energy (ΔGsolv) for each inhibitor was calculated in Gaussian09 with Polarizable Continuum Model (PCM) at the MP2/def2TZVP level of theory as single-point calculations performed for the inhibitor structures used throughout the computational protocol of interaction energy calculation. Integral equation formalism variant (IEFPCM) [63,64,65] was applied along with ExternalIteration [66, 67], DoVacuum, and SMD [68] options.
Evaluation of the results
The performance of each scoring model was estimated with Pearson correlation coefficient (R) calculated with respect to the experimentally determined inhibitory activity reported by Ren et al. [32]. Dissociation constant Kd was chosen as a reference, since all thermodynamical data provided by Ren et al. [32] were obtained with the recombinant catalytic domain of human PDE5, whereas IC50 values referred to the full-length of rabbit PDE5. As a further performance measure, the statistical predictor Npred, representing the success rate of prediction of relative affinities, was applied. Npred is calculated among all pairs of inhibitors as the percentage of concordant pairs with relative stability of the same sign as in the reference experimentally determined binding potency [69]. To enable comparison with the non-empirical interaction energy results, assigning lower binding energy values to more potent inhibitors, the scoring functions with higher score implying the greater binding potency were assigned the opposite values of the calculated correlation coefficient.
Results and discussion
Non-empirical models of PDE5 inhibitors binding
PDE5 binding site model considered herein includes 11 amino acid residues surrounding the inhibitor in the vicinity of 5 Å within the atom occupying pyrimidinone 5-position (Fig. 2). Asp764 is the only charged residue within this set. Except for Tyr612, Gln775, and Gln817, the remaining neutral residues are nonpolar. Total binding energies calculated by summing up the values characterizing individual amino acid residue-inhibitor dimers are given in Table 1. The corresponding EMP2 pairwise interaction energy values are provided in the Supplementary Material.
Negative values of total interaction energy calculated at the reference MP2 level of theory, EMP2, indicate the favorable interaction with PDE5 binding site. The prevalent contribution to binding energy at MP2 level of theory appears to be due to electrostatic interaction, \(E_{EL}^{(10)}\) (Table 1). The positive values of the first-order Heitler-London energy, E(10), imply the overestimated repulsion, as the short-range exchange contribution, \(E_{EX}^{(10)}\), is particularly sensitive to any structural inaccuracies [70]. Accounting for the subsequent delocalization term, \(E_{DEL}^{(R0)}\), restores the stabilizing nature of the resulting ESCF interaction energy (Table 1), however, the proper ranking of inhibitors is established only after including the correlation contribution, \(E_{CORR}^{(2)}\). Except for 5-Br and 5-Cl inhibitors, featuring the same value of EMP2 binding energy, the ranking of PDE5 inhibitors resulting from interaction energy calculated at MP2 level of theory (Table 1) is in agreement with experimental data [32].
To quantify the performance of particular interaction energy models, the relationship between the binding energy calculated at the subsequent levels of theory and the experimentally determined dissociation constants [32] was expressed by means of the Pearson correlation coefficient (R) and the success rate of prediction of relative inhibitory potency (Npred; see Table 1). The experimental binding affinities are reproduced only when dispersion is accounted for, i.e., at the EMP2 (R = − 0.95,Npred = 100%) or \(E_{EL,MTP}^{(10)}+E_{Das}\) (R = − 0.93,Npred = 90%) levels of theory. In fact, the predictive ability of the \(E_{EL,MTP}^{(10)}+E_{Das}\) model is essentially due to the dispersion contribution, EDas. As demonstrated by low values of R correlation coefficients, neither multipole electrostatic nor the first-order electrostatic interaction energy alone provide reliable predictions of the binding affinity in terms of the correlation with the experiment (R = − 0.45 and − 0.54 for \(E_{EL,MTP}^{(10)}\) and \(E_{EL}^{(10)}\), respectively; Table 1). Npred values associated with these levels of theory (70% in both cases) are lower compared to the models accounting for the dispersion interactions, mostly due to the underestimation of the 5-Cl inhibitor affinity.
Binding energies calculated at E(10) level of theory anti-correlate with the experimental inhibitory activity, i.e., higher repulsion is associated with stronger inhibitors (Table 1). Due to exponential dependence of \(E_{EX}^{(10)}\) contribution (accounted for at the E(10) level of theory) on the interatomic distance, any minor structural deficiencies resulting from, e.g., insufficiently accurate starting geometries, significantly affect this particular interaction energy term [20, 70]. The extent to which \(E_{EX}^{(10)}\) term is altered seems to be influenced by the overall binding strength, as the more potent ligands tend to experience higher repulsion originating from overestimated \(E_{EX}^{(10)}\) contribution. In consequence, delocalization energy, \(E_{DEL}^{(R0)}\), that contributes to the subsequent level of theory, namely ESCF, is not able to recover the proper inhibitory activity trend, as ESCF remains anti-correlated with respect to the experimental ligand affinity (Table 1). Finally, the inverse ranking of PDE5 inhibitors is overcome with \(E_{CORR}^{(2)}\) contribution and the resulting EMP2 energy provides a reasonable estimate of the inhibitory activity. The important role of the correlation term appears to be in agreement with the results of Riley et al. [71] demonstrating that halogen bonding involves the interplay between electrostatic and dispersion forces. As already emphasized, the performance of the approximate \(E_{EL,MTP}^{(10)}+E_{Das}\) model, accounting for long-range interaction energy terms only, is comparable to the predictive capabilities of EMP2 binding energy. However, the computational cost of \(E_{EL,MTP}^{(10)}+E_{Das}\) model is equivalent to force field calculations [20], as it scales with the square number of atoms, O(A2), while computational scaling of EMP2 energy is described by the fifth power of the number of atomic orbitals, O(N5).
The higher the atomic number of the substituting halogen (Cl, Br, I), the stronger the interaction with the protein, as described by nonempirical \(E_{EL,MTP}^{(10)}+E_{Das}\) model. This is in agreement with other available studies of halogen series, where it has been shown that the strength of the halogen bond increases with the increasing polarizability of the halogen atoms [72]. In particular, the increased halogen atom polarizability has been related to the magnitude of σ-hole potential [73, 74]. At the MP2 level of theory, the interaction energy values of both 5-Cl and 5-Br inhibitors are roughly the same, probably due to the incorporation of the short-range contributions, especially \(E_{EX}^{(10)}\). As a result, accounting for delocalization and correlation terms improves the overall performance of the EMP2 energy, but the binding affinities of the aforementioned compounds are still more similar than expected.
Solvation effects
In our recent work on protein-protein interaction (PPI) inhibitors [23, 33], ΔGsolv has emerged as an important factor determining whether significant correlation between relative interaction energy values and experimental results could be obtained. Estimation of binding affinity with interaction energy-based models is justified only if the remaining contributions to binding free energy are essentially constant across the inhibitor series. Relatively high standard deviation (SD) values of solvation free energy characterizing a set of ligands indicate that no predictions could be made based on the interaction energy alone. Accordingly, any correlation with the experimental binding affinity other than coincidental should not be expected. To overcome this issue, the inhibitor set could be limited to compounds featuring similar ΔGsolv values.
As for the analyzed PDE5 inhibitors, the correlation between experimental affinity and EMP2 or \(E_{EL,MTP}^{(10)}+E_{Das}\) binding energy already is significant. Nevertheless, it seemed important to verify solvation effects here as well, especially given the findings of Fanfrlik et al. [34], where binding energy of halogenated ligands yielded no correlation with experiment unless solvation energy was accounted for. The relatively low SD values of ΔGsolv calculated for PDE5 inhibitors with SMD model (Table 2) further support the predictive abilities of interaction energy-based models tested herein. The coincidence of low ΔGsolv standard deviation and the significant R correlation coefficient values for PDE5 inhibitors is in line with our previous analyses (see Refs. [23] and [33]), showing that in order to obtain correlation for the interaction energy, the associated solvent effects estimated with SMD model need to be relatively similar within the ligand set.
As reported by Ren et al. [32], entropy contribution related to 5-F inhibitor is substantially higher than for the remaining compounds (such observations are also described in the literature [75, 76]). In fact,the entropy term characterizing 5-F compound exceeds the corresponding enthalpic term, while the latter constitutes the main binding free energy contribution for the remaining PDE5 inhibitors within this set [32]. Interestingly, this is not reflected in our solvation energy calculations, where 5-F compound is associated with the lowest ΔGsolv value compared to the remaining inhibitors (see Table 2). Moreover, low SD values calculated for the analyzed PDE5Is indicate similar solvation effects across the series of ligands. Possibly, the methods for computing the solvation free energy of ligands are only approximate, yielding ΔGsolv results that suffer from significant errors [77,78,79], and, as a result, they do not represent the experimental results accurately enough.
Insights into the nature of the PDE5Is’ binding
EMP2 binding energy values obtained for individual amino acid residues or residue pairs (Asp764-Leu765, Ala767-Ile768, and Ile778-Ala779) are shown in Fig. 3. The highest interaction energy values are due to Gln817 residue (− 13.4 to − 14.6kcal ⋅mol− 1, Table S1 in the Supplementary Material). Phe820 residue also contributes significantly to the overall binding energy (− 4.7 to − 6.4kcal ⋅mol− 1). Interestingly, the interaction with Tyr612 residue, that seems to be halogen bonded to the inhibitor molecules, is of moderate strength (− 0.7 to − 2.4kcal ⋅mol− 1). Apparently, interactions with the remaining closely positioned amino acid residues (Gln817, Phe820, Val782, Asp764, and Leu765; Fig. 3) contribute together to the overall binding strength. In this case, focusing on a particular interaction (i.e., halogen/hydrogen bonding with Tyr612), while neglecting the contributions arising from the other nearby residues during novel inhibitors design should be reconsidered. In particular, our results do not support the decreased binding strength in 5-I-Tyr612 complex, accompanied by enhanced interaction with a buried water molecule, as observed by Ren et al. [32]. Further discussion on these differences is provided in the Supplementary Material.
As mentioned in the previous section, Ren et al. [32] reported substantial difference in the experimentally derived enthalpic/entropic contribution to the binding of 5-F inhibitor, as compared to the remaining compounds. It was suggested that the low enthalpic term of 5-F compound could be due to the electron withdrawal from the pyrimidinone ring upon the fluorine substitution, which probably results in reduced interactions with PDE5 amino acid residues, e.g., Phe820 or Val782 [32]. This was not confirmed by computational results reported herein (see Fig. 3 and Table S1 in the Supplementary Material), as binding of 5-F due to Phe820 and Val782 residues is of a similar magnitude as in the case of 5-H ligand.
Another possible reason pointed out by Ren et al. [32] for the decreased enthalpic contribution of 5-F inhibitor involved weakening of the bond between Tyr612 and Gln817 residues and 5-F compound. As shown in Fig. 3, 5-F–Tyr612 interaction amounting to − 0.7kcal ⋅mol− 1 is indeed weaker than for the remaining compounds, featuring EMP2 interaction energy within the range of − 1.3 to − 2.4kcal ⋅mol− 1 (see Table S1 in the Supplementary Material). However, 5-F–Gln817 EMP2 contribution is comparable to the corresponding value characterizing 5-H compound (− 13.4 and − 13.5kcal ⋅mol− 1, respectively). Another weaker interaction of the 5-F compound is the one with the Asp764-Leu765 residue pair (see Table S1 in Supplementary Material for details). It seems that the latter, together with the Tyr612 interaction, contribute to the somewhat smaller enthalpic term of 5-F inhibitor. Weakening of the interaction upon fluorine substitution might be due to prevalent negative electrostatic potential characterizing fluorine atoms covalently bonded to aromatic systems [80]. Nevertheless, the difference in the enthalpic/entropic contribution of 5-F compound is not as pronounced as one could expect given the experimental data.
Empirical scoring
Scoring functions applied for the evaluation of PDE5Is included GoldScore, ChemPLP (ChemPLPG), ChemScore, and Astex Statistical Potential (ASP) from Gold 5.6.3, ChemPLP (ChemPLPP), and PLP scoring functions from Plants, XBSF from VinaXB, as well as GlideSP and GlideXP from Glide. Correlation coefficients associated with all these methods are compared in Fig. 4. The values of each empirical score obtained for PDE5 inhibitors, along with the corresponding correlation coefficients and Npred values, are provided in Table S2 in the Supplementary Material.
Out of the nine scoring functions tested herein, only GoldScore function was capable of ranking the PDE5 inhibitors in sufficient agreement with the experimental inhibitory activity (R = − 0.96,Npred = 90%, see Fig. 4). In most cases, the affinity of 5-Cl and/or 5-F inhibitor is severely overestimated (see Table S2 in Supplementary Material). This might indicate that the majority of the currently used scoring functions lack the ability to properly describe halogen interactions, and there still is a room for improvement. Overall, the performance of scoring functions appear to be different depending on the system studied [23] and there is no apparent choice of an appropriate scoring approach to be used with a certain protein-ligand system. However, our nonempirical \(E_{EL,MTP}^{(10)}+E_{Das}\) model has been shown to possess the satisfactory predictive abilities for a variety of enzyme-inhibitor complexes [20, 21] including the PPI inhibitors [22, 23, 33]. Our current results seem to extend the range of application of \(E_{EL,MTP}^{(10)}+E_{Das}\) model to halogen-bearing compounds.
Summary
In this work, a congeneric series of halogenated PDE5 inhibitors, sharing a common pyrimidinone scaffold, was evaluated with the ab initio HVPT energy decomposition scheme, followed by the analysis with an approximate \(E_{EL,MTP}^{(10)}+E_{Das}\) model, including only multipole electrostatic and dispersion interaction energy terms. The latter method of the assessment of the inhibitory activity was compared to a number of empirical scoring approaches.
The interaction energy calculated for such a model accurately reproduces the experimentally derived affinities, as the correlation coefficients with respect to experimental inhibitory activity are equal to − 0.95 and − 0.93 for the reference EMP2 and nonempirical \(E_{EL,MTP}^{(10)}+E_{Das}\) binding energies, respectively. It must be noted that the experimental results are reproduced by the nonempirical models only when dispersion forces are accounted for. In the case of the computationally affordable \(E_{EL,MTP}^{(10)}+E_{Das}\) model, the damped dispersion expression EDas was included in the revised version, containing the halogen atom parameters. As for the empirical scoring, only GoldScore performance is adequate (R correlation coefficient of − 0.96) and comparable to either an approximate \(E_{EL,MTP}^{(10)}+E_{Das}\) or the most robust EMP2 level of theory.
The performance achieved for the \(E_{EL,MTP}^{(10)}+E_{Das}\) model with the revised EDas parameters suggests that it could be successfully used in systems that involve halogenated species. It is especially important in the view of mostly unsatisfactory results obtained with the empirical scoring functions, which demonstrates that reliable description of such systems remains a challenge.
Notes
Unpublished results (2018)
References
Metrangolo P, Resnati G (2001) Halogen bonding: a paradigm in supramolecular chemistry. Chem Eur J 7(12):2511–2519
Politzer P, Murray J (2013) Halogen bonding: an interim discussion. Chem Phys Chem 14(2):278–294
Kolar M, Hobza P (2016) Computer modeling of halogen bonds and other sigma-hole interactions. Chem Rev 116(9):5155–5187
Politzer P, Murray J (2018) Sigma-holes and pi-holes: similarities and differences. J Comp Chem 39(9):464–471
Bartashevich E, Tsirelson V (2018) A comparative view on the potential acting on an electron in a molecule and the electrostatic potential through the typical halogen bonds. J Comp Chem 39(10):573–580
Politzer P, Murray JS, Clark T (2010) Halogen bonding: an electrostatically-driven highly directional noncovalent interaction. Phys Chem Chem Phys 12:7748–7757
Wilcken R, Zimmermann M, Lange A, Zahn S, Boeckler F (2012) Using halogen bonds to address the protein backbone: a systematic evaluation. J Comput Aid Mol Des 26(8):935–945
Lu Y, Liu Y, Xu Z, Li H, Liu H, Zhu W (2012) Halogen bonding for rational drug design and new drug discovery. Expert Opin Drug Discov 7(5):375–383
Xu Z, Yang Z, Liu Y, Lu Y, Chen K, Zhu W (2014) Halogen bond: Its role beyond drug-target binding affinity for drug discovery and development. J Chem Inf Model 54(1):69–78
Kolar M, Tabarrini O (2017) Halogen bonding in nucleic acid complexes. J Med Chem 60(21):8681–8690
Sirimulla S, Bailey J, Vegesna R, Narayan M (2013) Halogen interactions in protein-ligand complexes: implications of halogen bonding for rational drug design. J Chem Inf Model 53(11):2781–2791
Kurczab R (2017) The evaluation of QM/MM-driven molecular docking combined with MM/GBSA calculations as a halogen-bond scoring strategy. Acta Cryst B 73:188–194
Ho P (2017) Halogen bonding in medicinal chemistry: from observation to prediction. Future Med Chem 9(7):637–640
Koebel M, Schmadeke G, Posner R, Sirimulla, S (2016) Autodock VinaXB: implementation of XBSF, new empirical halogen bond scoring function, into AutoDock Vina. J Cheminf 8:27
Wolters L, Schyman P, Pavan M, Jorgensen W, Bickelhaupt F, Kozuch S (2014) The many faces of halogen bonding: a review of theoretical models and methods. Wiley Interdiscip Rev Comput Mol Sci 4(6):523–540
Wang G, Zhu W (2016) Molecular docking for drug discovery and development: a widely used approach but far from perfect. Future Med Chem 8(14)
Ford M, Ho P (2016) Computational tools to model halogen bonds in medicinal chemistry. J Med Chem 59(5):1655–1670
Titov O, Shulga D, Palyulin V, Zefirov N (2016) Quadrupole correction for halogen bonding description in virtual screening and molecular docking. Dokl Chem 471:338–342
Lin F, MacKerell A (2018) Polarizable empirical force field for halogen-containing compounds based on the classical Drude oscillator. J Chem Theory Comput 14(2):1083–1098
Giedroyć-Piasecka W, Dyguda-Kazimierowicz E, Beker W, Mor M, Lodola A, Sokalski WA (2014) Physical nature of fatty acid amide hydrolase interactions with its inhibitors: Testing a simple nonempirical scoring model. J Phys Chem B 118(51):14727–14736
Jedwabny W, Panecka-Hofman J, Dyguda-Kazimierowicz E, Wade RC, Sokalski WA (2017a) Application of a simple quantum chemical approach to ligand fragment scoring for Trypanosoma brucei pteridine reductase 1 inhibition. J Comput Aid Mol Des 31:715–728
Jedwabny W, Kłossowski S, Purohit T, Cierpicki T, Grembecka J, Dyguda-Kazimierowicz E (2017b) Theoretical models of inhibitory activity for inhibitors of protein-protein interactions: targeting menin-Mixed Lineage Leukemia with small molecules. Med Chem Commun 8:2216–2227
Jedwabny W, Lodola A, Dyguda-Kazimierowicz E, Sokalski WA (2018) Theoretical model of EphA2 – Ephrin A1 inhibition. Molecules 23:1688–1707
Podeszwa R, Pernal K, Patkowski K, Szalewicz K (2010) Extension of the Hartree-Fock plus dispersion method by first-order correlation effects. J Phys Chem Lett 1(2):550–555
Pernal K, Podeszwa R, Patkowski K, Szalewicz K (2009) Dispersionless density functional theory. Phys Rev Lett 103(26): 263201
Shahbaz M, Szalewicz K (2018) Do semilocal density-functional approximations recover dispersion energies at small intermonomer separations?. Phys Rev Lett 121(11). Article Number: 113402
Pauls M, Moynihan B, Barrick T, Kruuse C, Madigan J, Hainsworth A, Isaacs J (2018) The effect of phosphodiesterase-5 inhibitors on cerebral blood flow in humans: a systematic review. J Cereb Blood Flow Metab 38(2):189–203
Kee C, Ge X, Gilard V, Malet-Martino M, Low M (2018) A review of synthetic phosphodiesterase type 5 inhibitors (PDE-5i) found as adulterants in dietary supplements. J Pharm Biomed Anal 147:250–277
Korkmaz-Icoz S, Radovits T, Szabo G (2018) Targeting phosphodiesterase 5 as a therapeutic option against myocardial ischaemia/reperfusion injury and for treating heart failure. Br J Pharmacol 175(2):223–231
Barone I, Giordano C, Bonofiglio D, Ando S, Catalano S (2017) Phosphodiesterase type 5 and cancers: progress and challenges. Oncotarget 8(58):99179–99202
Xu Z, Liu Z, Chen T, Chen T, Wang Z, Tian G, Shi J, Wang X, Lu Y, Yan X, Wang G, Jiang H, Chen K, Wang S, Xu Y, Shen J, Zhu W (2011) Utilization of halogen bond in lead optimization: a case study of rational design of potent phosphodiesterase type 5 (PDE5) inhibitors. J Med Chem 54 (15):5607–5611
Ren J, He Y, Chen W, Chen T, Wang G, Wang Z, Xu Z, Luo X, Zhu W, Jiang H, Shen J, Xu Y (2014) Thermodynamic and structural characterization of halogen bonding in protein-ligand interactions: a case study of PDE5 and its inhibitors. J Med Chem 57(8):3588–3593
Jedwabny W, Cierpicki T, Grembecka J, Dyguda-Kazimierowicz E (2018) Validation of approximate nonempirical scoring model for menin-mixed lineage leukemia inhibitors. Theor Chem Acc 137(11). Article Number: 148
Fanfrlik J, Kolar M, Kamlar M, Hurny D, Ruiz F, Cousido-Siah A, Mitschler A, Rezac J, Munusamy E, Lepsik M, Matejicek P, Vesely J, Podjarny A, Hobza P (2013) Modulation of aldose reductase inhibition by halogen bond tuning. ACS Chem Biol 8(11):2484–2492
Li H, Robertson A, Jensen J (2005) Very fast empirical prediction and interpretation of protein pKa values. Proteins 61:704–721
Bas D, Rogers D, Jensen J (2008) Very fast prediction and rationalization of pKa values for protein-ligand complexes. Proteins 73:765–783
Olsson M, Søndergard C, Rostkowski M, Jensen J (2011) PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions. J Chem Theory Comput 7:525–537
Søndergard C, Olsson M, Rostkowski M, Jensen J (2011) Improved treatment of ligands and coupling effects in empirical calculation and rationalization of pKa values. J Chem Theory Comput 7:2284–2295
Mahoney M, Jorgensen W (2000) A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions. J Chem Phys 112(20):8910–8922
Brooks B, Bruccoleri R, Olafson B, States D, Swaminathan S, Karplus M (1983) CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. J Comput Chem 4(2):187–217
Vanommeslaeghe K, Hatcher E, Acharya C, Kundu S, Zhong S, Shim J, Darian E, Guvench O, Lopes P, Vorobyov I, MacKerell A Jr (2010) CHARMM General force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J Comput Chem 31(4):671–690
MacKerell A, Bashford D, Bellott M, Dunbrack R, Evanseck J, Field M, Fischer S, Gao J, Guo H, Ha S, Joseph-McCarthy D, Kuchnir L, Kuczera K, Lau F, Mattos C, Michnick S, Ngo T, Nguyen D, Prodhom B, Reiher W, Roux B, Schlenkrich M, Smith J, Stote R, Straub J, Watanabe M, Wiorkiewicz-Kuczera J, Yin D, Karplus M (1998) All-atom empirical potential for molecular modeling and dynamics studies of proteins. J Phys Chem B 102(18):3586–3616
Mackerell A, Feig M, Brooks C (2004) Extending the treatment of backbone energetics in protein force fields: Limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. J Comput Chem 25(11):1400–1415
Stote R, Karplus M (1995) Zinc binding in proteins and solution: a simple but accurate nonbonded representation. Proteins 23(1):12–31
Vanommeslaeghe K, MacKerell A Jr (2012) Automation of the CHARMM general force field (CGenFF) I: bond perception and atom typing. J Chem Inf Model 52(12):3144–3154
Vanommeslaeghe K, Raman EP, MacKerell A Jr (2012) Automation of the CHARMM general force field (CGenFF) II: assignment of bonded parameters and partial atomic charges. J Chem Inf Model 52(12):3155–3168
Yu W, He X, Vanommeslaeghe K, MacKerell AD Jr (2012) Extension of the CHARMM general force field to sulfonyl-containing compounds and its utility in biomolecular simulations. J Comput Chem 33(31):2451–2468
Maestro version 9.3 (2012) Schrödinger, LLC, New York
Banks J, Beard H, Cao Y, Cho A, Damm W, Farid R, Felts A, Halgren T, Mainz D, Maple J, Murphy R, Philipp D, Repasky M, Zhang L, Berne B, Friesner R, Gallicchio E, Levy R (2005) Integrated modeling program, applied chemical theory (IMPACT). J Comp Chem 26:1752–1780
Sokalski WA, Roszak S, Pecul K (1988) An efficient procedure for decomposition of the SCF interaction energy into components with reduced basis set dependence. Chem Phys Lett 153:153– 159
Góra R, Sokalski WA, Leszczyński J, Pett V (2005) The nature of interactions in the ionic crystal of 3-pentenenitrile, 2-nitro-5-oxo, ion(-1) sodium. J Phys Chem B 109:2027–2033
Schmidt MW, Baldridge KK, Boatz JA, Elbert ST, Gordon MS, Jensen JH, Koseki S, Matsunaga N, Nguyen KA, Su SJ, Windus TL, Dupuis M, Montgomery JA (1993) General atomic and molecular electronic structure system. J Comput Chem 14(11):1347–1363
Weigend F (2006) Accurate coulomb-fitting basis sets for H to Rn. Phys Chem Chem Phys 8(9):1057–1065
Weigend F, Ahlrichs R (2005) Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: design and assessment of accuracy. Phys Chem Chem Phys 7(18):3297–3305
Boys S, Bernardi F (1970) The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors. Mol Phys 19:553–566
Sokalski W, Poirier R (1983) Cumulative atomic multipole representation of the molecular charge distribution and its basis set dependence. Chem Phys Lett 98(1):86–92
Jones G, Willett P, Glen R, Leach A, Taylor R (1997) Development and validation of a genetic algorithm for flexible docking. J Mol Biol 267:727–748
Korb O, Stutzle T, Exner T (2009) Empirical scoring functions for advanced protein-ligand docking with PLANTS. J Chem Inf Model 49(1):84–96
(2018) The Cambridge Crystallographic Data Centre (CCDC). Available online: https://www.ccdc.cam.ac.uk/support-and-resources/support/case/?caseid=5d1a2fc0-c93a-49c3-a8e2-f95c472dcff0 (accessed on 11.10.2018)
Friesner R, Banks J, Murphy R, Halgren T, Klicic J, Mainz D, Repasky M, Knoll E, Shelley M, Perry J, Shaw D, Francis P, Shenkin P (2004) Glide: a new approach for rapid, accurate docking and scoring. 1. Method and assessment of docking accuracy. J Med Chem 47(7):1739–1749
Friesner RA, Murphy RB, Repasky MP, Frye LL, Greenwood JR, Halgren TA, Sanschagrin PC, Mainz DT (2006) Extra precision glide: docking and scoring incorporating a model of hydrophobic enclosure for protein-ligand complexes. J Med Chem 49(21):6177–6196
Schrödinger LLC (2018) Schrödinger Release 2018-1, Glide, New York
Tomasi J, Mennucci B, Cances E (1999) The IEF version of the PCM solvation method: an overview of a new method addressed to study molecular solutes at the QM ab initio level. J Mol Struct Theochem 464(1-3):211–226
Pascualahuir J, Silla E, Tunon I (1994) GEPOL: An improved description of molecular surfaces. III. a new algorithm for the computation of a solventexcluding surface. J Comp Chem 15(10):1127–1138
Tomasi J, Mennucci B, Cammi R (2005) Quantum mechanical continuum solvation models. Chem Rev 105(8):2999–3093
Improta R, Scalmani G, Frisch MJ, Barone V (2007) Toward effective and reliable fluorescence energies in solution by a new state specific polarizable continuum model time dependent density functional theory approach. J Chem Phys 127(7). Article Number: 074504
Improta R, Barone V, Scalmani G, Frisch MJ (2006) A state-specific polarizable continuum model time dependent density functional theory method for excited state calculations in solution. J Chem Phys 125(5). Article Number: 054103
Marenich AV, Cramer CJ, Truhlar DG (2009) Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions. J Phys Chem B 113(18):6378–6396
Langner K, Beker W, Sokalski W (2012) Robust predictive power of the electrostatic term at shortened intermolecular distances. J Phys Chem Lett 3:2785–2789
Grzywa R, Dyguda-Kazimierowicz E, Sieńczyk M, Feliks M, Sokalski WA, Oleksyszyn J (2007) The molecular basis of urokinase inhibition: from the nonempirical analysis of intermolecular interaction to the prediction of binding affinity. J Mol Model 13:677–683
Riley K, Murray J, Fanfrlik J, Rezac J, Sola R, Concha M, Ramos F, Politzer P (2013) Halogen bond tunability II: The varying roles of electrostatic and dispersion contributions to attraction in halogen bonds. J Mol Model 19(11):4651–4659
Auffinger P, Hays F, Westhof E, Ho P (2004) Halogen bonds in biological molecules. Proc Natl Acad Sci USA 101(48):16789–16794
Murray J, Macaveiu L, Politzer P (2014) Factors affecting the strengths of sigma-hole electrostatic potentials. J Comput Sci 5(4):590–596
Bundhun A, Ramasami P, Murray J, Politzer P (2013) Trends in sigma-hole strengths and interactions of F3MX molecules (M = C, Si, Ge and X = F, Cl, Br, I). J Mol Model 19(7):2739–2746
Oda M, Saito M, Tsumuraya T, Fujii L (2010) Contribution of the trifluoroacetyl group in the thermodynamics of antigen-antibody binding. J Mol Recognit 23(3):263–270
Krintel C, Francotte P, Pickering D, Juknaite L, Pohlsgaard J, Olsen L, Frydenvang K, Goffin E, Pirotte B, Kastrup J (2016) Enthalpy-entropy compensation in the binding of modulators at ionotropic glutamate receptor GluA2. Biophys J 110(11):2397–2406
Ángyán JG (1992) Common theoretical framework for quantum chemical solvent effect theories. J Math Chem 10(1):93–137
Ho JM, Klamt A, Coote ML (2010) Comment on the correct use of continuum solvent models. J Phys Chem A 114(51):13442–13444
Yokogawa D (2018) Toward accurate solvation free energy calculation with the reference interaction site model self-consistent field: introduction of a new bridge function. J Chem Theory Comput 14:3272–3278
Riley K, Murray J, Fanfrlik J, Rezac J, Sola R, Concha M, Ramos F, Politzer P (2011) Halogen bond tunability I: The effects of aromatic fluorine substitution on the strengths of halogen-bonding interactions involving chlorine, bromine, and iodine. J Mol Model 17(12):3309–3318
Acknowledgements
We thank Professor W. Andrzej Sokalski for thoughtful reading and comments. The project was financed by the National Science Centre Poland (grant number 2016/21/N/ST4/00516). Calculations were carried out using resources provided by Wrocław Centre for Networking and Supercomputing, Poznań Supercomputing and Networking Center, and Warsaw Interdisciplinary Centre for Mathematical and Computational Modelling. We also thank Wrocław University of Science and Technology for support.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher’s note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This paper belongs to Topical Collection 8th conference on Modeling & Design of Molecular Materials (MDMM 2018)
Electronic supplementary material
Below is the link to the electronic supplementary material.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Jedwabny, W., Dyguda-Kazimierowicz, E. Revisiting the halogen bonding between phosphodiesterase type 5 and its inhibitors. J Mol Model 25, 29 (2019). https://doi.org/10.1007/s00894-018-3897-z
Received:
Accepted:
Published:
DOI: https://doi.org/10.1007/s00894-018-3897-z