Evaluation of the utility of homology models in high throughput docking.

High throughput docking (HTD) is routinely used for in silico screening of compound libraries with the aim to find novel leads in a drug discovery program. In the absence of an experimentally determined structure, a homology model can be used instead. Here we present an assessment of the utility of homology models in HTD by docking 300,000 anticipated inactive compounds along with 642 known actives into the binding site of the insulin-like growth factor 1 receptor (IGF-1R) kinase constructed by homology modeling. Twenty-one different templates were selected and the enrichment curves obtained by the homology models were compared to those obtained by three IGF-1R crystal structures. The results show a wide range of enrichments from random to as good as two of the three IGF-1R crystal structures. Nevertheless, if we consider the enrichment obtained at 2% of the database screened as a performance criterion, the best crystal structure outperforms the best homology model. Surprisingly, the sequence identity of the template to the target is not a good descriptor to predict the enrichment obtained by a homology model. The three homology models that yield the worst enrichment have the smallest binding-site volume. Based on our results, we propose ensemble docking to perform HTD with homology models.

Assessing scoring functions for protein-ligand interactions.

An assessment of nine scoring functions commonly applied in docking using a set of 189 protein-ligand complexes is presented. The scoring functions include the CHARMm potential, the scoring function DrugScore, the scoring function used in AutoDock, the three scoring functions implemented in DOCK, as well as three scoring functions implemented in the CScore module in SYBYL (PMF, Gold, ChemScore). We evaluated the abilities of these scoring functions to recognize near-native configurations among a set of decoys and to rank binding affinities. Binding site decoys were generated by molecular dynamics with restraints. To investigate whether the scoring functions can also be applied for binding site detection, decoys on the protein surface were generated. The influence of the assignment of protonation states was probed by either assigning "standard" protonation states to binding site residues or adjusting protonation states according to experimental evidence. The role of solvation models in conjunction with CHARMm was explored in detail. These include a distance-dependent dielectric function, a generalized Born model, and the Poisson equation. We evaluated the effect of using a rigid receptor on the outcome of docking by generating all-pairs decoys ("cross-decoys") for six trypsin and seven HIV-1 protease complexes. The scoring functions perform well to discriminate near-native from misdocked conformations, with CHARMm, DOCK-energy, DrugScore, ChemScore, and AutoDock yielding recognition rates of around 80%. Significant degradation in performance is observed in going from decoy to cross-decoy recognition for CHARMm in the case of HIV-1 protease, whereas DrugScore and ChemScore, as well as CHARMm in the case of trypsin, show only small deterioration. In contrast, the prediction of binding affinities remains problematic for all of the scoring functions. ChemScore gives the highest correlation value with R(2) = 0.51 for the set of 189 complexes and R(2) = 0.43 for the set of 116 complexes that does not contain any of the complexes used to calibrate this scoring function. Neither a more accurate treatment of solvation nor a more sophisticated charge model for zinc improves the quality of the results. Improved modeling of the protonation states, however, leads to a better prediction of binding affinities in the case of the generalized Born and the Poisson continuum models used in conjunction with the CHARMm force field.

Weak temperature dependence of the free energy surface and folding pathways of structured peptides.

The thermodynamics and energetics of a 20-residue synthetic peptide with a stable three-stranded antiparallel beta-sheet fold are investigated by implicit solvent molecular dynamics (MD) at 330 K (slightly above the melting temperature in the model) and compared with previous simulation results at 360 K. At both temperature values, the peptide folds reversibly to the NMR solution conformation, irrespective of the starting conformation. The sampling of the conformational space (2.3 micros and 25 folding events at 330 K, and 3 micros and 50 folding events at 360 K) is sufficient to obtain a thermodynamic description of minima and transition states on the free energy surface, which is determined near equilibrium by counting populations. The free energy surface, plotted as a function of two-order parameters that monitor formation of either of the beta-hairpins, is similar at both temperature values. The statistically predominant folding pathway and its frequency (about two-thirds of the folding events) are the same at 330 K and 360 K. Furthermore, the main unfolding route is the reverse of the predominant folding pathway. The effective energy and its electrostatic and van der Waals contributions show a downhill profile at both temperatures, implying that the free energy barrier is of entropic origin and corresponds to the freezing of about two-thirds of the chain into a beta-hairpin conformation. The average folding rate is nearly the same at 330 K and 360 K, while the unfolding rate is about four times slower at 330 K than at 360 K. Taken together with previous MD analysis of alpha-helices and beta-hairpins, the present simulation results indicate that the free energy surface and folding mechanism of structured peptides have a weak temperature dependence.

Evaluation of a fast implicit solvent model for molecular dynamics simulations.

A solvation term based on the solvent accessible surface area (SASA) is combined with the CHARMM polar hydrogen force field for the efficient simulation of peptides and small proteins in aqueous solution. Only two atomic solvation parameters are used: one is negative for favoring the direct solvation of polar groups and the other positive for taking into account the hydrophobic effect on apolar groups. To approximate the water screening effects on the intrasolute electrostatic interactions, a distance-dependent dielectric function is used and ionic side chains are neutralized. The use of an analytical approximation of the SASA renders the model extremely efficient (i.e., only about 50% slower than in vacuo simulations). The limitations and range of applicability of the SASA model are assessed by simulations of proteins and structured peptides. For the latter, the present study and results reported elsewhere show that with the SASA model it is possible to sample a significant amount of folding/unfolding transitions, which permit the study of the thermodynamics and kinetics of folding at an atomic level of detail.

Calculation of protein ionization equilibria with conformational sampling: pK(a) of a model leucine zipper, GCN4 and barnase.

The use of conformational ensembles provided by nuclear magnetic resonance (NMR) experiments or generated by molecular dynamics (MD) simulations has been regarded as a useful approach to account for protein motions in the context of pK(a) calculations, yet the idea has been tested occasionally. This is the first report of systematic comparison of pK(a) estimates computed from long multiple MD simulations and NMR ensembles. As model systems, a synthetic leucine zipper, the naturally occurring coiled coil GCN4, and barnase were used. A variety of conformational averaging and titration curve-averaging techniques, or combination thereof, was adopted and/or modified to investigate the effect of extensive global conformational sampling on the accuracy of pK(a) calculations. Clustering of coordinates is proposed as an approach to reduce the vast diversity of MD ensembles to a few structures representative of the average electrostatic properties of the system in solution. Remarkable improvement of the accuracy of pK(a) predictions was achieved by the use of multiple MD simulations. By using multiple trajectories the absolute error in pK(a) predictions for the model leucine zipper was reduced to as low as approximately 0.25 pK(a) units. The validity, advantages, and limitations of explicit conformational sampling by MD, compared with the use of an average structure and a high internal protein dielectric value as means to improve the accuracy of pK(a) calculations, are discussed.