Predicting absolute ligand binding free energies to a simple model site

A central challenge in structure-based ligand design is the accurate prediction of binding free energies. Here we apply alchemical free energy calculations in explicit solvent to predict ligand binding in a model cavity in T4 lysozyme. Even in this simple site, there are challenges. We made systematic improvements, beginning with single poses from docking, then including multiple poses, additional protein conformational changes, and using an improved charge model. Computed absolute binding free energies had an RMS error of 1.9 kcal/mol relative to previously determined experimental values. In blind prospective tests, the methods correctly discriminated between several true ligands and decoys in a set of putative binders identified by docking. In these prospective tests, the RMS error in predicted binding free energies relative to those subsequently determined experimentally was only 0.6 kcal/mol. X-ray crystal structures of the new ligands bound in the cavity corresponded closely to predictions from the free energy calculations, but sometimes differed from those predicted by docking. Finally, we examined the impact of holding the protein rigid, as in docking, with a view to learning how approximations made in docking affect accuracy and how they may be improved.     

Protein surface roughness accounts for binding free energy of Plasmepsin II‐ligand complexes

The calculation of absolute binding affinities for protein‐inhibitor complexes remains as one of the main challenges in computational structure‐based ligand design. The present work explored the calculations of surface fractal dimension (as a measure of surface roughness) and the relationship with experimental binding free energies of Plasmepsin II complexes. Plasmepsin II is an attractive target for novel therapeutic compounds to treat malaria. However, the structural flexibility of this enzyme is a drawback when searching for specific inhibitors. Concerning that, we performed separate explicitly solvated molecular dynamics simulations using the available high‐resolution crystal structures of different Plasmepsin II complexes. Molecular dynamics simulations allowed a better approximation to systems dynamics and, therefore, a more reliable estimation of surface roughness. This constitutes a novel approximation in order to obtain more realistic values of fractal dimension, because previous works considered only x‐ray structures. Binding site fractal dimension was calculated considering the ensemble of structures generated at different simulation times. A linear relationship between binding site fractal dimension and experimental binding free energies of the complexes was observed within 20 ns. Previous studies of the subject did not uncover this relationship. Regression model, coined FD model, was built to estimate binding free energies from binding site fractal dimension values. Leave‐one‐out cross‐validation showed that our model reproduced accurately the absolute binding free energies for our training set (R² = 0.76; <|error|> =0.55 kcal/mol; SD_error = 0.19 kcal/mol). The fact that such a simple model may be applied raises some questions that are addressed in the article.     

Fragmentation-tree density representation for crystallographic modelling of bound ligands

The identification and modelling of ligands into macromolecular models is important for understanding molecule's function and for designing inhibitors to modulate its activities. We describe new algorithms for the automated building of ligands into electron density maps in crystal structure determination. Location of the ligand-binding site is achieved by matching numerical shape features describing the ligand to those of density clusters using a "fragmentation-tree" density representation. The ligand molecule is built using two distinct algorithms exploiting free atoms with inter-atomic connectivity and Metropolis-based optimisation of the conformational state of the ligand, producing an ensemble of structures from which the final model is derived. The method was validated on several thousand entries from the Protein Data Bank. In the majority of cases, the ligand-binding site could be correctly located and the ligand model built with a coordinate accuracy of better than 1 ?. We anticipate that the method will be of routine use to anyone modelling ligands, lead compounds or even compound fragments as part of protein functional analyses or drug design efforts.     

Accurate determination of the binding free energy for KcsA-charybdotoxin complex from the potential of mean force calculations with restraints

Free energy calculations for protein-ligand dissociation have been tested and validated for small ligands (50 atoms or less), but there has been a paucity of studies for larger, peptide-size ligands due to computational limitations. Previously we have studied the energetics of dissociation in a potassium channel-charybdotoxin complex by using umbrella sampling molecular-dynamics simulations, and established the need for carefully chosen coordinates and restraints to maintain the physiological ligand conformation. Here we address the ligand integrity problem further by constructing additional potential of mean forces for dissociation of charybdotoxin using restraints. We show that the large discrepancies in binding free energy arising from simulation artifacts can be avoided by using appropriate restraints on the ligand, which enables determination of the binding free energy within the chemical accuracy. We make several suggestions for optimal choices of harmonic potential parameters and restraints to be used in binding studies of large ligands.     

Discrete restraint-based protein modeling and the Calpha-trace problem

We present a novel de novo method to generate protein models from sparse, discretized restraints on the conformation of the main chain and side chain atoms. We focus on Calpha-trace generation, the problem of constructing an accurate and complete model from approximate knowledge of the positions of the Calpha atoms and, in some cases, the side chain centroids. Spatial restraints on the Calpha atoms and side chain centroids are supplemented by constraints on main chain geometry, phi/xi angles, rotameric side chain conformations, and inter-atomic separations derived from analyses of known protein structures. A novel conformational search algorithm, combining features of tree-search and genetic algorithms, generates models consistent with these restraints by propensity-weighted dihedral angle sampling. Models with ideal geometry, good phi/xi angles, and no inter-atomic overlaps are produced with 0.8 A main chain and, with side chain centroid restraints, 1.0 A all-atom root-mean-square deviation (RMSD) from the crystal structure over a diverse set of target proteins. The mean model derived from 50 independently generated models is closer to the crystal structure than any individual model, with 0.5 A main chain RMSD under only Calpha restraints and 0.7 A all-atom RMSD under both Calpha and centroid restraints. The method is insensitive to randomly distributed errors of up to 4 A in the Calpha restraints. The conformational search algorithm is efficient, with computational cost increasing linearly with protein size. Issues relating to decoy set generation, experimental structure determination, efficiency of conformational sampling, and homology modeling are discussed.     

Ligand-induced conformational changes: improved predictions of ligand binding conformations and affinities

A computational docking strategy using multiple conformations of the target protein is discussed and evaluated. A series of low molecular weight, competitive, nonpeptide protein tyrosine phosphatase inhibitors are considered for which the x-ray crystallographic structures in complex with protein tyrosine phosphatase 1B (PTP1B) are known. To obtain a quantitative measure of the impact of conformational changes induced by the inhibitors, these were docked to the active site region of various structures of PTP1B using the docking program FlexX. Firstly, the inhibitors were docked to a PTP1B crystal structure cocrystallized with a hexapeptide. The estimated binding energies for various docking modes as well as the RMS differences between the docked compounds and the crystallographic structure were calculated. In this scenario the estimated binding energies were not predictive inasmuch as docking modes with low estimated binding energies corresponded to relatively large RMS differences when aligned with the corresponding crystal structure. Secondly, the inhibitors were docked to their parent protein structures in which they were cocrystallized. In this case, there was a good correlation between low predicted binding energy and a correct docking mode. Thirdly, to improve the predictability of the docking procedure in the general case, where only a single target protein structure is known, we evaluate an approach which takes possible protein side-chain conformational changes into account. Here, side chains exposed to the active site were considered in their allowed rotamer conformations and protein models containing all possible combinations of side-chain rotamers were generated. To evaluate which of these modeled active sites is the most likely binding site conformation for a certain inhibitor, the inhibitors were docked against all active site models. The receptor rotamer model corresponding to the lowest estimated binding energy is taken as the top candidate. Using this protocol, correct inhibitor binding modes could successfully be discriminated from proposed incorrect binding modes. Moreover, the ranking of the estimated ligand binding energies was in good agreement with experimentally observed binding affinities.     

High resolution X-ray structures of mouse major urinary protein nasal isoform in complex with pheromones

In mice, the major urinary proteins (MUP) play a key role in pheromonal communication by binding and transporting semiochemicals. MUP‐IV is the only isoform known to be expressed in the vomeronasal mucosa. In comparison with the MUP isoforms that are abundantly excreted in the urine, MUP‐IV is highly specific for the male mouse pheromone 2‐sec‐butyl‐4,5‐dihydrothiazole (SBT). To examine the structural basis of this ligand preference, we determined the X‐ray crystal structure of MUP‐IV bound to three mouse pheromones: SBT, 2,5‐dimethylpyrazine, and 2‐heptanone. We also obtained the structure of MUP‐IV with 2‐ethylhexanol bound in the cavity. These four structures show that relative to the major excreted MUP isoforms, three amino acid substitutions within the binding calyx impact ligand coordination. The F103 for A along with F54 for L result in a smaller cavity, potentially creating a more closely packed environment for the ligand. The E118 for G substitution introduces a charged group into a hydrophobic environment. The sidechain of E118 is observed to hydrogen bond to polar groups on all four ligands with nearly the same geometry as seen for the water‐mediated hydrogen bond network in the MUP‐I and MUP‐II crystal structures. These differences in cavity size and interactions between the protein and ligand are likely to contribute to the observed specificity of MUP‐IV.     

Flexible docking using tabu search and an empirical estimate of binding affinity

This article describes the implementation of a new docking approach. The method uses a Tabu search methodology to dock flexibly ligand molecules into rigid receptor structures. It uses an empirical objective function with a small number of physically based terms derived from fitting experimental binding affinities for crystallographic complexes. This means that docking energies produced by the searching algorithm provide direct estimates of the binding affinities of the ligands. The method has been tested on 50 ligand-receptor complexes for which the experimental binding affinity and binding geometry are known. All water molecules are removed from the structures and ligand molecules are minimized in vacuo before docking. The lowest energy geometry produced by the docking protocol is within 1.5 Å root-mean square of the experimental binding mode for 86% of the complexes. The lowest energies produced by the docking are in fair agreement with the known free energies of binding for the ligands. Proteins 33:367–382, 1998. © 1998 Wiley-Liss, Inc.     

First principles predictions of the structure and function of g-protein-coupled receptors: validation for bovine rhodopsin

G-protein-coupled receptors (GPCRs) are involved in cell communication processes and with mediating such senses as vision, smell, taste, and pain. They constitute a prominent superfamily of drug targets, but an atomic-level structure is available for only one GPCR, bovine rhodopsin, making it difficult to use structure-based methods to design receptor-specific drugs. We have developed the MembStruk first principles computational method for predicting the three-dimensional structure of GPCRs. In this article we validate the MembStruk procedure by comparing its predictions with the high-resolution crystal structure of bovine rhodopsin. The crystal structure of bovine rhodopsin has the second extracellular (EC-II) loop closed over the transmembrane regions by making a disulfide linkage between Cys-110 and Cys-187, but we speculate that opening this loop may play a role in the activation process of the receptor through the cysteine linkage with helix 3. Consequently we predicted two structures for bovine rhodopsin from the primary sequence (with no input from the crystal structure)-one with the EC-II loop closed as in the crystal structure, and the other with the EC-II loop open. The MembStruk-predicted structure of bovine rhodopsin with the closed EC-II loop deviates from the crystal by 2.84 A coordinate root mean-square (CRMS) in the transmembrane region main-chain atoms. The predicted three-dimensional structures for other GPCRs can be validated only by predicting binding sites and energies for various ligands. For such predictions we developed the HierDock first principles computational method. We validate HierDock by predicting the binding site of 11-cis-retinal in the crystal structure of bovine rhodopsin. Scanning the whole protein without using any prior knowledge of the binding site, we find that the best scoring conformation in rhodopsin is 1.1 A CRMS from the crystal structure for the ligand atoms. This predicted conformation has the carbonyl O only 2.82 A from the N of Lys-296. Making this Schiff base bond and minimizing leads to a final conformation only 0.62 A CRMS from the crystal structure. We also used HierDock to predict the binding site of 11-cis-retinal in the MembStruk-predicted structure of bovine rhodopsin (closed loop). Scanning the whole protein structure leads to a structure in which the carbonyl O is only 2.85 A from the N of Lys-296. Making this Schiff base bond and minimizing leads to a final conformation only 2.92 A CRMS from the crystal structure. The good agreement of the ab initio-predicted protein structures and ligand binding site with experiment validates the use of the MembStruk and HierDock first principles' methods. Since these methods are generic and applicable to any GPCR, they should be useful in predicting the structures of other GPCRs and the binding site of ligands to these proteins.     

Solution-state molecular structure of apo and oleate-liganded liver fatty acid-binding protein

Rat liver fatty acid-binding protein (LFABP) is distinctive among intracellular lipid-binding proteins (iLBPs): more than one molecule of long-chain fatty acid and a variety of diverse ligands can be bound within its large cavity, and in vitro lipid transfer to model membranes follows a mechanism that is diffusion-controlled rather than mediated by protein-membrane collisions. Because the apoprotein has proven resistant to crystallization, nuclear magnetic resonance spectroscopy offers a unique route to functionally informative comparisons of molecular structure and dynamics for LFABP in free (apo) and liganded (holo) forms. We report herein the solution-state structures determined for apo-LFABP at pH 6.0 and for holoprotein liganded to two oleates at pH 7.0, as well as the structure of the complex including locations of the ligands. 1H, 13C, and 15N resonance assignments revealed very similar types and locations of secondary structural elements for apo- and holo-LFABP as judged from chemical shift indices. The solution-state tertiary structures of the proteins were derived with the CNS/ARIA computational protocol, using distance and angular restraints based on 1H-1H nuclear Overhauser effects (NOEs), hydrogen-bonding networks, 3J(HNHA) coupling constants, intermolecular NOEs, and residual dipolar (NH) couplings. The holo-LFABP solution-state conformation is in substantial agreement with a previously reported X-ray structure [Thompson, J., Winter, N., Terwey, D., Bratt, J., and Banaszak, L. (1997) The crystal structure of the liver fatty acid-binding protein. A complex with two bound oleates, J. Biol. Chem. 272, 7140-7150], including the typical beta-barrel capped by a helix-turn-helix portal. In the solution state, the internally bound oleate has the expected U-shaped conformation and is tethered electrostatically, but the extended portal ligand can adopt a range of conformations based on the computationally refined structures, in contrast to the single conformation observed in the crystal structure. The apo-LFABP also has a well-defined beta-barrel structural motif typical of other members of the iLBP protein family, but the portal region that is thought to facilitate ligand entry and exit exhibits conformational variability and an unusual "open cap" orientation with respect to the barrel. These structural results allow us to propose a model in which ligand binding to LFABP occurs through conformational fluctuations that adjust the helix-turn-helix motif to open or close the top of the beta-barrel, and solvent accessibility to the protein cavity favors diffusion-controlled ligand transport.     

Docking multiple conformations of a flexible ligand into a protein binding site using NMR restraints

A method is described for docking a large, flexible ligand using intra-ligand conformational restraints from exchange-transferred NOE (etNOE) data. Numerous conformations of the ligand are generated in isolation, and a subset of representative conformations is selected. A crude model of the protein-ligand complex is used as a template for overlaying the selected ligand structures, and each complex is conformationally relaxed by molecular mechanics to optimize the interaction. Finally, the complexes were assessed for structural quality. Alternative approaches are described for the three steps of the method: generation of the initial docking template; selection of a subset of ligand conformations; and conformational sampling of the complex. The template is generated either by manual docking using interactive graphics or by a computational grid-based search of the binding site. A subset of conformations from the total number of peptides calculated in isolation is selected based on either low energy and satisfaction of the etNOE restraints, or a cluster analysis of the full set. To optimize the interactions in the complex, either a restrained Monte Carlo-energy minimization (MCM) protocol or a restrained simulated annealing (SA) protocol were used. This work produced 53 initial complexes of which 8 were assessed in detail. With the etNOE conformational restraints, all of the approaches provide reasonable models. The grid-based approach to generate an initial docking template allows a large volume to be sampled, and as a result, two distinct binding modes were identified for a fifteen-residue peptide binding to an enzyme active site.     

X-ray structures of small ligand-FKBP complexes provide an estimate for hydrophobic interaction energies

A new crystal form of native FK506 binding protein (FKBP) has been obtained which has proved useful in ligand binding studies. Three different small molecule ligand complexes and the native enzyme have been determined at higher resolution than 2.0 A. Dissociation constants of the related small molecule ligands vary from 20 mM for dimethylsulphoxide to 200 microM for tetrahydrothiophene 1-oxide. Comparison of the four available crystal structures shows that the protein structures are identical to within experimental error, but there are differences in the water structure in the active site. Analysis of the calculated buried surface areas of these related ligands provides an estimated van der Waals contribution to the binding energy of -0.5 kJ/A(2) for non-polar interactions between ligand and protein.     

Contributions to the binding free energy of ligands to avidin and streptavidin

The free energy of binding of a ligand to a macromolecule is here formally decomposed into the (effective) energy of interaction, reorganization energy of the ligand and the macromolecule, conformational entropy change of the ligand and the macromolecule, and translational and rotational entropy loss of the ligand. Molecular dynamics simulations with implicit solvation are used to evaluate these contributions in the binding of biotin, biotin analogs, and two peptides to avidin and streptavidin. We find that the largest contribution opposing binding is the protein reorganization energy, which is calculated to be from 10 to 30 kcal/mol for the ligands considered here. The ligand reorganization energy is also significant for flexible ligands. The translational/rotational entropy is 4.5-6 kcal/mol at 1 M standard state and room temperature. The calculated binding free energies are in the correct range, but the large statistical uncertainty in the protein reorganization energy precludes precise predictions. For some complexes, the simulations show multiple binding modes, different from the one observed in the crystal structure. This finding is probably due to deficiencies in the force field but may also reflect considerable ligand flexibility.     

Non-Boltzmann thermodynamic integration (NBTI) for macromolecular systems: relative free energy of binding of trypsin to benzamidine and benzylamine

The relative free energies of binding of trypsin to two amine inhibitors, benzamidine (BZD) and benzylamine (BZA), were calculated using non-Boltzmann thermodynamic integration (NBTI). Comparison of the simulations with the crystal structures of both complexes, trypsin-BZD and trypsin-BZA, shows that NBTI simulations better sample conformational space relative to thermodynamic integration (TI) simulations. The relative binding free energy calculated using NBTI was much closer to the experimentally determined value than that obtained using TI. The error in the TI simulation was found to be primarily due to incorrect sampling of BZA's conformation in the binding pocket. In contrast, NBTI produces a smooth mutation from BZD to BZA using a surrogate potential, resulting in a much closer agreement between the inhibitors' conformations and the omit electron density maps. This superior agreement between experiment and simulation, of both relative binding free energy differences and conformational sampling, demonstrates NBTI's usefulness for free energy calculations in macromolecular simulations.     

Predicting Ligand Binding Affinity with Alchemical Free Energy Methods in a Polar Model Binding Site

We present a combined experimental and modeling study of organic ligand molecules binding to a slightly polar engineered cavity site in T4 lysozyme (L99A/M102Q). For modeling, we computed alchemical absolute binding free energies. These were blind tests performed prospectively on 13 diverse, previously untested candidate ligand molecules. We predicted that eight compounds would bind to the cavity and five would not; 11 of 13 predictions were correct at this level. The RMS error to the measurable absolute binding energies was 1.8 kcal/mol. In addition, we computed “relative” binding free energies for six phenol derivatives starting from two known ligands: phenol and catechol. The average RMS error in the relative free energy prediction was 2.5 kcal/mol (phenol) and 1.1 kcal/mol (catechol). To understand these results at atomic resolution, we obtained x-ray co-complex structures for nine of the diverse ligands and for all six phenol analogs. The average RMSD of the predicted pose to the experiment was 2.0 Å (diverse set), 1.8 Å (phenol-derived predictions), and 1.2 Å (catechol-derived predictions). We found that predicting accurate affinities and rank-orderings required near-native starting orientations of the ligand in the binding site. Unanticipated binding modes, multiple ligand binding, and protein conformational change all proved challenging for the free energy methods. We believe that these results can help guide future improvements in physics-based absolute binding free energy methods.     

Ensemble docking of multiple protein structures: considering protein structural variations in molecular docking

One approach to incorporate protein flexibility in molecular docking is the use of an ensemble consisting of multiple protein structures. Sequentially docking each ligand into a large number of protein structures is computationally too expensive to allow large-scale database screening. It is challenging to achieve a good balance between docking accuracy and computational efficiency. In this work, we have developed a fast, novel docking algorithm utilizing multiple protein structures, referred to as ensemble docking, to account for protein structural variations. The algorithm can simultaneously dock a ligand into an ensemble of protein structures and automatically select an optimal protein structure that best fits the ligand by optimizing both ligand coordinates and the conformational variable m, where m represents the m-th structure in the protein ensemble. The docking algorithm was validated on 10 protein ensembles containing 105 crystal structures and 87 ligands in terms of binding mode and energy score predictions. A success rate of 93% was obtained with the criterion of root-mean-square deviation <2.5 A if the top five orientations for each ligand were considered, comparable to that of sequential docking in which scores for individual docking are merged into one list by re-ranking, and significantly better than that of single rigid-receptor docking (75% on average). Similar trends were also observed in binding score predictions and enrichment tests of virtual database screening. The ensemble docking algorithm is computationally efficient, with a computational time comparable to that for docking a ligand into a single protein structure. In contrast, the computational time for the sequential docking method increases linearly with the number of protein structures in the ensemble. The algorithm was further evaluated using a more realistic ensemble in which the corresponding bound protein structures of inhibitors were excluded. The results show that ensemble docking successfully predicts the binding modes of the inhibitors, and discriminates the inhibitors from a set of noninhibitors with similar chemical properties. Although multiple experimental structures were used in the present work, our algorithm can be easily applied to multiple protein conformations generated by computational methods, and helps improve the efficiency of other existing multiple protein structure(MPS)-based methods to accommodate protein flexibility.     

Evaluation of site-directed spin labeling for characterizing protein-ligand complexes using simulated restraints

Simulation studies have been performed to evaluate the utility of site-directed spin labeling for determining the structures of protein-ligand complexes, given a known protein structure. Two protein-ligand complexes were used as model systems for these studies: a 1.9-A-resolution x-ray structure of a dihydrofolate reductase mutant complexed with methotrexate, and a 1.5-A-resolution x-ray structure of the V-Src tyrosine kinase SH2 domain complexed with a five-residue phosphopeptide. Nitroxide spin labels were modeled at five dihydrofolate reductase residue positions and at four SH2 domain residue positions. For both systems, after energy minimization, conformational ensembles of the spin-labeled residues were generated by simulated annealing while holding the remainder of the protein-ligand complex fixed. Effective distances, simulating those that could be obtained from (1)H-NMR relaxation measurements, were calculated between ligand protons and the spin labels. These were converted to restraints with several different levels of precision. Restrained simulated annealing calculations were then performed with the aim of reproducing target ligand-binding modes. The effects of incorporating a few supplementary short-range (< or =5.0 A) distance restraints were also examined. For the dihydrofolate reductase-methotrexate complex, the ligand-binding mode was reproduced reasonably well using relatively tight spin-label restraints, but methotrexate was poorly localized using loose spin-label restraints. Short-range and spin-label restraints proved to be complementary. For the SH2 domain-phosphopeptide complex without the short-range restraints, the peptide did not localize to the correct depth in the binding groove; nevertheless, the orientation and internal conformation of the peptide was reproduced moderately well. Use of the spin-label restraints in conjunction with the short-range restraints resulted in relatively well defined structural ensembles. These results indicate that restraints derived from site-directed spin labeling can contribute significantly to defining the orientations and conformations of bound ligands. Accurate ligand localization appears to require either a few supplementary short-range distance restraints, or relatively tight spin-label restraints, with at least one spin label positioned so that some of the restraints draw the ligand into the binding pocket in the latter case.