Comparison of end-point continuum-solvation methods for the calculation of protein-ligand binding free energies

We have compared the predictions of ligand‐binding affinities from several methods based on end‐point molecular dynamics simulations and continuum solvation, that is, methods related to MM/PBSA (molecular mechanics combined with Poisson–Boltzmann and surface area solvation). Two continuum‐solvation models were considered, viz., the Poisson–Boltzmann (PB) and generalised Born (GB) approaches. The nonelectrostatic energies were also obtained in two different ways, viz., either from the sum of the bonded, van der Waals, nonpolar solvation energies, and entropy terms (as in MM/PBSA), or from the scaled protein–ligand van der Waals interaction energy (as in the linear interaction energy approach, LIE). Three different approaches to calculate electrostatic energies were tested, viz., the sum of electrostatic interaction energies and polar solvation energies, obtained either from a single simulation of the complex or from three independent simulations of the complex, the free protein, and the free ligand, or the linear‐response approximation (LRA). Moreover, we investigated the effect of scaling the electrostatic interactions by an effective internal dielectric constant of the protein (?_int). All these methods were tested on the binding of seven biotin analogues to avidin and nine 3‐amidinobenzyl‐1H‐indole‐2‐carboxamide inhibitors to factor Xa. For avidin, the best results were obtained with a combination of the LIE nonelectrostatic energies with the MM+GB electrostatic energies from a single simulation, using ?_int = 4. For fXa, standard MM/GBSA, based on one simulation and using ?_int = 4–10 gave the best result. The optimum internal dielectric constant seems to be slightly higher with PB than with GB solvation. © Proteins 2012; © 2012 Wiley Periodicals, Inc.     

Discovery of Rab1 binding sites using an ensemble of clustering methods

Targeting non‐native‐ligand binding sites for potential investigative and therapeutic applications is an attractive strategy in proteins that share common native ligands, as in Rab1 protein. Rab1 is a subfamily member of Rab proteins, which are members of Ras GTPase superfamily. All Ras GTPase superfamily members bind to native ligands GTP and GDP, that switch on and off the proteins, respectively. Rab1 is physiologically essential for autophagy and transport between endoplasmic reticulum and Golgi apparatus. Pathologically, Rab1 is implicated in human cancers, a neurodegenerative disease, cardiomyopathy, and bacteria‐caused infectious diseases. We have performed structural analyses on Rab1 protein using a unique ensemble of clustering methods, including multi‐step principal component analysis, non‐negative matrix factorization, and independent component analysis, to better identify representative Rab1 proteins than the application of a single clustering method alone does. We then used the identified representative Rab1 structures, resolved in multiple ligand states, to map their known and novel binding sites. We report here at least a novel binding site on Rab1, involving Rab1‐specific residues that could be further explored for the rational design and development of investigative probes and/or therapeutic small molecules against the Rab1 protein. Proteins 2017; 85:859–871. © 2016 Wiley Periodicals, Inc.     

Prediction of the binding energy for small molecules,peptides and proteins

A fast and reliable evaluation of the binding energy from a single conformation of a molecular complex is an important practical task. Knowledge‐based scoring schemes may not be sufficiently general and transferable, while molecular dynamics or Monte Carlo calculations with explicit solvent are too computationally expensive for many applications. Recently, several empirical schemes using finite difference Poisson–Boltzmann electrostatics to predict energies for particular types of complexes were proposed. Here, an improved empirical binding energy function has been derived and validated on three different types of complexes: protein–small ligand, protein–peptide and protein–protein. The function uses the boundary element algorithm to evaluate the electrostatic solvation energy. We show that a single set of parameters can predict the relative binding energies of the heterogeneous validation set of complexes with 2.5 kcal/mol accuracy. We also demonstrate that global optimization of the ligand and of the flexible side‐chains of the receptor improves the accuracy of the evaluation. Copyright © 1999 John Wiley & Sons, Ltd.     

A model binding site for testing scoring functions in molecular docking

Prediction of interaction energies between ligands and their receptors remains a major challenge for structure-based inhibitor discovery. Much effort has been devoted to developing scoring schemes that can successfully rank the affinities of a diverse set of possible ligands to a binding site for which the structure is known. To test these scoring functions, well-characterized experimental systems can be very useful. Here, mutation-created binding sites in T4 lysozyme were used to investigate how the quality of atomic charges and solvation energies affects molecular docking. Atomic charges and solvation energies were calculated for 172,118 molecules in the Available Chemicals Directory using a semi-empirical quantum mechanical approach by the program AMSOL. The database was first screened against the apolar cavity site created by the mutation Leu99Ala (L99A). Compared to the electronegativity-based charges that are widely used, the new charges and desolvation energies improved ranking of known apolar ligands, and better distinguished them from more polar isosteres that are not observed to bind. To investigate whether the new charges had predictive value, the non-polar residue Met102, which forms part of the binding site, was changed to the polar residue glutamine. The structure of the resulting Leu99Ala and Met102Gln double mutant of T4 lysozyme (L99A/M102Q) was determined and the docking calculation was repeated for the new site. Seven representative polar molecules that preferentially docked to the polar versus the apolar binding site were tested experimentally. All seven bind to the polar cavity (L99A/M102Q) but do not detectably bind to the apolar cavity (L99A). Five ligand-bound structures of L99A/M102Q were determined by X-ray crystallography. Docking predictions corresponded to the crystallographic results to within 0.4A RMSD. Improved treatment of partial atomic charges and desolvation energies in database docking appears feasible and leads to better distinction of true ligands. Simple model binding sites, such as L99A and its more polar variants, may find broad use in the development and testing of docking algorithms.     

Calculation of the total electrostatic energy of a macromolecular system: solvation energies, binding energies, and conformational analysis

In this report we describe an accurate numerical method for calculating the total electrostatic energy of molecules of arbitrary shape and charge distribution, accounting for both Coulombic and solvent polarization terms. In addition to the solvation energies of individual molecules, the method can be used to calculate the electrostatic energy associated with conformational changes in proteins as well as changes in solvation energy that accompany the binding of charged substrates. The validity of the method is examined by calculating the hydration energies of acetate, methyl ammonium, ammonium, and methanol. The method is then used to study the relationship between the depth of a charge within a protein and its interaction with the solvent. Calculations of the relative electrostatic energies of crystal and misfolded conformations of Themiste dyscritum hemerythrin and the VL domain of an antibody are also presented. The results indicate that electrostatic charge-solvent interactions strongly favor the crystal structures. More generally, it is found that charge-solvent interactions, which are frequently neglected in protein structure analysis, can make large contributions to the total energy of a macromolecular system.     

Development, validation, and application of adapted PEOE charges to estimate pKa values of functional groups in protein-ligand complexes

For routine pK(a) calculations of protein-ligand complexes in drug design, the PEOE method to compute partial charges was modified. The new method is applicable to a large scope of proteins and ligands. The adapted charges were parameterized using experimental free energies of solvation of amino acids and small organic ligands. For a data set of 80 small organic molecules, a correlation coefficient of r(2) = 0.78 between calculated and experimental solvation free energies was obtained. Continuum electrostatics pK(a) calculations based on the Poisson-Boltzmann equation were carried out on a validation set of nine proteins for which 132 experimental pK(a) values are known. In total, an overall RMSD of 0.88 log units between calculated and experimentally determined data is achieved. In particular, the predictions of significantly shifted pK(a) values are satisfactory, and reasonable estimates of protonation states in the active sites of lysozyme and xylanase could be obtained. Application of the charge-assignment and pK(a)-calculation procedure to protein-ligand complexes provides clear structural interpretations of experimentally observed changes of protonation states of functional groups upon complex formation. This information is essential for the interpretation of thermodynamic data of protein-ligand complex formation and provides the basis for the reliable factorization of the free energy of binding in enthalpic and entropic contributions. The modified charge-assignment procedure forms the basis for future automated pK(a) calculations of protein-ligand complexes.     

Model assembly study of the ligand binding by p-hydroxybenzoate hydroxylase: Correlation between the calculated binding energies and the experimental dissociation constants

The energies of binding of seven ligands by p-hydroxybenzoate hydroxylase (PHBH) were calculated theoretically. Direct enzyme–ligand interaction energies were calculated using the ab initio quantum mechanical model assembly of the active site at the 3-21G level. Solvation energies of the ligands needed in the evaluation of the binding energies were calculated with the semiempirical AM1–SM2 method and the long-range electrostatic interaction energies between the ligands and the protein matrix classically using the static charge distributions of the ligands and the protein. Energies for proton-transfer between the ligands OH or SH substituent at position 4 and the active-site tyrosine within the ab initio model assemblies were calculated and compared to the corresponding pK_as in aqueous solution. Excluding 3,4-dihydroxybenzoate, the natural product of PHBH, a linear relationship between the calculated binding energies and the experimental binding free energies was found with a correlation coefficient of 0.90. Contributions of the direct enzyme–ligand interaction energies, solvation energies and the long-range electrostatic interaction energies to the calculated binding energies were analyzed. The proton-transfer energies of the ligands with substituents ortho to the ionized OH were found to be perturbed less in the model calculations than the energies of their meta isomers as deduced from the corresponding pK_as. © 1995 Wiley-Liss, Inc.     

Conformational energy range of ligands in protein crystal structures: The difficult quest for accurate understanding

In this review, we address a fundamental question: What is the range of conformational energies seen in ligands in protein‐ligand crystal structures? This value is important biophysically, for better understanding the protein‐ligand binding process; and practically, for providing a parameter to be used in many computational drug design methods such as docking and pharmacophore searches. We synthesize a selection of previously reported conflicting results from computational studies of this issue and conclude that high ligand conformational energies really are present in some crystal structures. The main source of disagreement between different analyses appears to be due to divergent treatments of electrostatics and solvation. At the same time, however, for many ligands, a high conformational energy is in error, due to either crystal structure inaccuracies or incorrect determination of the reference state. Aside from simple chemistry mistakes, we argue that crystal structure error may mainly be because of the heuristic weighting of ligand stereochemical restraints relative to the fit of the structure to the electron density. This problem cannot be fixed with improvements to electron density fitting or with simple ligand geometry checks, though better metrics are needed for evaluating ligand and binding site chemistry in addition to geometry during structure refinement. The ultimate solution for accurately determining ligand conformational energies lies in ultrahigh‐resolution crystal structures that can be refined without restraints.     

Field‐based comparison of ligand and coactivator binding sites of nuclear receptors

A structure‐based comparison of the ligand‐binding domains of 35 nuclear receptors from five different subfamilies is presented. Their ligand and coactivator binding sites are characterized using knowledge‐based contact preference fields for hydrophobic and hydrophilic interactions implemented in the MOE modeling environment. Additionally, for polar knowledge‐based field points the preference for negative or positive electrostatic interactions is estimated using the Poisson‐Boltzmann equation. These molecular‐interaction fields are used to cluster the nuclear receptor family based on similarities of their binding sites. By analyzing the similarities and differences of hydrophobic and polar fields in binding pockets of related receptors it is possible to identify conserved interactions in ligand and coactivator binding pockets, which support e.g. design of specific ligands during lead optimization or virtual screening as docking filter. Examples of remarkable similarities between ligand binding sites of members from phylogenetically different nuclear receptor families (RXR, RAR, HNF4, NR5) and differences between closely related subtypes (LXR, RAR, TR) are discussed in more detail. Significant similarities and differences of coactivator binding sites are shown for NR3Cs, LXRs and PPARs. © 2009 Wiley Periodicals, Inc. Biopolymers 91: 884–894, 2009. This article was originally published online as an accepted preprint. The “Published Online” date corresponds to the preprint version. You can request a copy of the preprint by emailing the Biopolymers editorial office at biopolymers@wiley.com     

Limited validity of group additivity for the folding energetics of the peptide group

The principle of group additivity is a standard feature of analyses of the energetics of protein folding, but it is known that it may not always be valid for the polar peptide group. The neighboring residue effect shows that group additivity is not strictly valid for a heteropeptide. We show here that group additivity fails seriously for peptide groups close to either peptide end, even for a homopeptide that has blocked end groups with no formal charges involved. The failure of group additivity is caused by the electrostatic character of the solvation of peptide polar groups and is illustrated with values of the electrostatic solvation free energy (ESF) calculated by DelPhi. Solvation free energies and enthalpies are known experimentally for monoamides and are often used to model the solvation of peptide groups, but ESF results show that monoamide values are very different from those of peptide groups. A main cause of the difference is that peptide solvation depends on the dipole-dipole interactions made between adjacent peptide groups, which vary with peptide conformation. Ligands that interact with the peptide backbone by an electrostatic mechanism could show a similar peptide end effect, and hydrogen exchange results from the literature confirm that exchange rates are position-dependent close to peptide ends.     

A comprehensive examination of the contributions to the binding entropy of protein–ligand complexes

One of the most important requirements in computer‐aided drug design is the ability to reliably evaluate the binding free energies. However, the process of ligand binding is very complex because of the intricacy of the interrelated processes that are difficult to predict and quantify. In fact, the deeper understanding of the origin of the observed binding free energies requires the ability to decompose these free energies to their contributions from different interactions. Furthermore, it is important to evaluate the relative entropic and enthalpic contributions to the overall free energy. Such an evaluation is useful for assessing temperature effects and exploring specialized options in enzyme design. Unfortunately, calculations of binding entropies have been much more challenging than calculations of binding free energies. This work is probably the first to present microscopic evaluation of all of the relevant components to the binding entropy, namely configurational, polar solvation, and hydrophobic entropies. All of these contributions are evaluated by the restraint release approach. The calculated results shed an interesting light on major compensation effects in both the solvation and hydrophobic effect and, despite some overestimate, can provide very useful insight. This study also helps in analyzing some problems with the widely used molecular mechanics/Poisson‐Boltzmann surface area approach. Proteins 2010. © 2010 Wiley‐Liss, Inc.     

De novo ligand design to an ensemble of protein structures

We describe a combinatorial method for de novo ligand design to an ensemble of receptor structures. Receptor conformations, protonation states, and structural water molecules are considered consistently within the framework of de novo ligand design. The method relies on Monte Carlo optimization to search the space of ligand structures, conformations, and rigid-body movements as well as receptor models. The method is applied to an ensemble of HIV protease and human collagenase receptor models. Ligand structures generated de novo exhibit the correct hydrogen-bonding pattern in the core of the active site, with hydrophobic groups extending into the receptor S1 and S1' pocket space. Furthermore, it is shown that known ligands are recovered in the correct binding mode and in the native, most tightly binding receptor model.     

An automated method for dynamic ligand design

An automated method for the dynamic ligand design (DLD) for a binding site of known structure is described. The method can be used for the creation of de novo ligands and for the modification of existing ligands. The binding site is saturated with atoms (sp³ carbon atoms in the present implementation) that form molecules under the influence of a potential function that joins atoms to each other with the correct stereochemistry. The resulting molecules are linked to precomputed functional group minimum energy positions in the binding site. The generalized potential function allows atoms to sample a continuous parameter space that includes the Cartesian coordinates and their occupancy and type, e.g., the method allows change of an sp³ carbon into an sp² carbon or oxygen. A parameter space formulated in this way can then be sampled and optimized by a variety of methods. In this work, molecules are generated by use of a Monte Carlo simulated annealing algorithm. The DLD method is illustrated by its application to the binding site of FK506 binding protein (FKBP), an immunophilin. De novo ligands are designed and modification of the immunosuppressant drug FK506 are suggested. The results demonstrate that the dynamic ligand design approach can automatically construct ligands which complement both the shape and charge distribution of the binding site. © 1995 Wiley-Liss, Inc.     

Spatial Separation of Molecular Conformers and Clusters

Gas-phase molecular physics and physical chemistry experiments commonly use supersonic expansions through pulsed valves for the production of cold molecular beams. However, these beams often contain multiple conformers and clusters, even at low rotational temperatures. We present an experimental methodology that allows the spatial separation of these constituent parts of a molecular beam expansion. Using an electric deflector the beam is separated by its mass-to-dipole moment ratio, analogous to a bender or an electric sector mass spectrometer spatially dispersing charged molecules on the basis of their mass-to-charge ratio. This deflector exploits the Stark effect in an inhomogeneous electric field and allows the separation of individual species of polar neutral molecules and clusters. It furthermore allows the selection of the coldest part of a molecular beam, as low-energy rotational quantum states generally experience the largest deflection. Different structural isomers (conformers) of a species can be separated due to the different arrangement of functional groups, which leads to distinct dipole moments. These are exploited by the electrostatic deflector for the production of a conformationally pure sample from a molecular beam. Similarly, specific cluster stoichiometries can be selected, as the mass and dipole moment of a given cluster depends on the degree of solvation around the parent molecule. This allows experiments on specific cluster sizes and structures, enabling the systematic study of solvation of neutral molecules.     

Blind Prediction of Charged Ligand Binding Affinities in a Model Binding Site

Predicting absolute protein–ligand binding affinities remains a frontier challenge in ligand discovery and design. This becomes more difficult when ionic interactions are involved because of the large opposing solvation and electrostatic attraction energies. In a blind test, we examined whether alchemical free-energy calculations could predict binding affinities of 14 charged and 5 neutral compounds previously untested as ligands for a cavity binding site in cytochrome c peroxidase. In this simplified site, polar and cationic ligands compete with solvent to interact with a buried aspartate. Predictions were tested by calorimetry, spectroscopy, and crystallography. Of the 15 compounds predicted to bind, 13 were experimentally confirmed, while 4 compounds were false negative predictions. Predictions had a root-mean-square error of 1.95 kcal/mol to the experimental affinities, and predicted poses had an average RMSD of 1.7 Å to the crystallographic poses. This test serves as a benchmark for these thermodynamically rigorous calculations at predicting binding affinities for charged compounds and gives insights into the existing sources of error, which are primarily electrostatic interactions inside proteins. Our experiments also provide a useful set of ionic binding affinities in a simplified system for testing new affinity prediction methods.     

Escherichia coli glutaminyl-tRNA synthetase is electrostatically optimized for binding of its cognate substrates

Natural evolution has resulted in protein molecules displaying a wide range of binding properties that include extremes of affinity and specificity. A detailed understanding of the principles underlying protein structure-function relationships, particularly with respect to binding properties, would greatly enhance molecular engineering and ligand design studies. Here, we have analyzed the interactions of an aminoacyl-tRNA synthetase for which strong evolutionary pressure has enforced high specificity for substrate binding and catalysis. Electrostatic interactions have been identified as one efficient mechanism for enhancing binding specificity; as such, the effects of charged and polar groups were the focus of this study. The binding of glutaminyl-tRNA synthetase from Escherichia coli to several ligands, including the natural substrates, was analyzed. The electrostatic complementarity of the enzyme to its ligands was assessed using measures derived from affinity optimization theory. The results were independent of the details of the calculational parameters, including the value used for the protein dielectric constant. Glutamine and ATP, two of the natural ligands, were found to be extremely complementary to their binding sites, particularly in regions seen to make electrostatic interactions in the structure. These data suggest that the optimization of electrostatic interactions has played an important role in guiding the evolution of this enzyme. The results also show that the enzyme is able to effectively select for high affinity and specificity for the same chemical moieties both in the context of smaller substrates, and in that of a larger reactive intermediate. The regions of greatest non-complementarity between the enzyme and ligands are the portions of the ligand that make few polar contacts with the binding site, as well as the sites of chemical reaction, where overly strong electrostatic binding interactions with the substrate could hinder catalysis. The results also suggest that the negative charge on the phosphorus center of glutaminyl-adenylate plays an important role in the tight binding of this intermediate, and thus that adenylate analogs that preserve the negative charge in this region may bind substantially tighter than analogs where this group is replaced with a neutral group, such as the sulfamoyl family, which can make similar hydrogen bonds but is uncharged.     

Fragment-based virtual screening approach and molecular dynamics simulation studies for identification of BACE1 inhibitor leads

Traditional structure-based virtual screening method to identify drug-like small molecules for BACE1 is so far unsuccessful. Location of BACE1, poor Blood Brain Barrier permeability and P-glycoprotein (Pgp) susceptibility of the inhibitors make it even more difficult. Fragment-based drug design method is suitable for efficient optimization of initial hit molecules for target like BACE1. We have developed a fragment-based virtual screening approach to identify/optimize the fragment molecules as a starting point. This method combines the shape, electrostatic, and pharmacophoric features of known fragment molecules, bound to protein conjugate crystal structure, and aims to identify both chemically and energetically feasible small fragment ligands that bind to BACE1 active site. The two top-ranked fragment hits were subjected for a 53 ns MD simulation. Principle component analysis and free energy landscape analysis reveal that the new ligands show the characteristic features of established BACE1 inhibitors. The potent method employed in this study may serve for the development of potential lead molecules for BACE1-directed Alzheimer’s disease therapeutics.     

Artificial protein cavities as specific ligand-binding templates: characterization of an engineered heterocyclic cation-binding site that preserves the evolved specificity of the parent protein.

Cavity complementation has been observed in many proteins, where an appropriate small molecule binds to a cavity-forming mutant. Here, the binding of compounds to the W191G cavity mutant of cytochrome c peroxidase is characterized by X-ray crystallography and binding thermodynamics. Unlike cavities created by removal of hydrophobic side-chains, the W191G cavity does not bind neutral or hydrophobic compounds, but displays a strong specificity for heterocyclic cations, consistent with the role of the protein to stabilize a tryptophan radical at this site. Ligand dissociation constants for the protonated cationic state ranged from 6 microM for 2-amino-5-methylthiazole to 1 mM for neutral ligands, and binding was associated with a large enthalpy-entropy compensation. X-ray structures show that each of 18 compounds with binding behavior bind specifically within the artificial cavity and not elsewhere in the protein. The compounds make multiple hydrogen bonds to the cavity walls using a subset of the interactions seen between the protein and solvent in the absence of ligand. For all ligands, every atom that is capable of making a hydrogen bond does so with either protein or solvent. The most often seen interaction is to Asp235, and most compounds bind with a specific orientation that is defined by their ability to interact with this residue. Four of the ligands do not have conventional hydrogen bonding atoms, but were nevertheless observed to orient their most polar CH bond towards Asp235. Two of the larger ligands induce disorder in a surface loop between Pro190 and Asn195 that has been identified as a mobile gate to cavity access. Despite the predominance of hydrogen bonding and electrostatic interactions, the small variation in observed binding free energies were not correlated readily with the strength, type or number of hydrogen bonds or with calculated electrostatic energies alone. Thus, as with naturally occurring binding sites, affinities to W191G are likely to be due to a subtle balance of polar, non-polar, and solvation terms. These studies demonstrate how cavity complementation and judicious choice of site can be used to produce a protein template with an unusual ligand-binding specificity.