Toward an optimal docking and free energy calculation scheme in ligand design with application to COX-1 inhibitors.

Cyclooxygenase-1 (COX-1) is one of the main targets of most pain-relieving pharmaceuticals. Although the enzyme is well characterized, it is known to be a difficult target for automated molecular docking and scoring. We collected from the literature a structurally diverse set of 45 nonsteroidal anti-inflammatory drugs (NSAIDs) and COX-2-selective inhibitors (coxibs) with a wide range of binding affinities for COX-1. The binding of this data set to a homology model of human COX-1 was analyzed with different combinations of molecular docking algorithms, scoring functions, and the linear interaction energy (LIE) method for estimating binding affinities. It is found that the computational protocols for estimation of binding affinities are extremely sensitive to the initial orientations of the ligands in the binding pocket. To overcome this limitation, we propose a systematic exploration of docking poses using the LIE calculations as a postscoring function. This scheme yields predictions in excellent agreement with experiment, with a mean unsigned error of 0.9 kcal/mol for binding free energies and structures of high quality. A significant improvement of the results is also seen when averaging over experimental data from several independent measurements.

Computational prediction of alanine scanning and ligand binding energetics in G-protein coupled receptors.

Site-directed mutagenesis combined with binding affinity measurements is widely used to probe the nature of ligand interactions with GPCRs. Such experiments, as well as structure-activity relationships for series of ligands, are usually interpreted with computationally derived models of ligand binding modes. However, systematic approaches for accurate calculations of the corresponding binding free energies are still lacking. Here, we report a computational strategy to quantitatively predict the effects of alanine scanning and ligand modifications based on molecular dynamics free energy simulations. A smooth stepwise scheme for free energy perturbation calculations is derived and applied to a series of thirteen alanine mutations of the human neuropeptide Y1 receptor and series of eight analogous antagonists. The robustness and accuracy of the method enables univocal interpretation of existing mutagenesis and binding data. We show how these calculations can be used to validate structural models and demonstrate their ability to discriminate against suboptimal ones.

Mutagenesis and computational modeling of human G-protein-coupled receptor Y2 for neuropeptide Y and peptide YY.

Neuropeptide Y and peptide YY receptor type 2 (Y2) is involved in appetite regulation and several other physiological processes. We have investigated the structure of the human Y2 receptor. Computational modeling of receptor-agonist interactions was used as a guide to design a series of receptor mutants, followed by binding assays using full-length and truncated peptide agonists and the Y2-specific antagonist BIIE0246. Our model suggested a hydrogen bond network among highly conserved residues Thr2.61, Gln3.32, and His7.39, which could play roles in ligand binding and/or receptor structure. In addition, the C-terminus of the peptide could make contact with residues Tyr5.38 and Leu6.51. Mutagenesis of all these positions, followed by binding assays, provides experimental support for our computational model: most of the mutants for the residues forming the proposed hydrogen bond network displayed reduced peptide agonist affinities as well as reduced hPYY3-36 potency in a functional assay. The Ala and Leu mutants of Gln3.32 and His7.39 disrupted membrane expression of the receptor. Combined with the modeling, the experimental results support roles for these hydrogen bond network residues in peptide binding as well as receptor architecture. The reduced agonist affinity for mutants of Tyr5.38 and Leu6.51 supports their role in a binding pocket surrounding the invariant tyrosine at position 36 of the peptide ligands. The results for antagonist BIIE0246 suggest several differences in interactions compared to those of the peptides. Our results lead to a new structural model for NPY family receptors and peptide binding.

Computer simulations of structure-activity relationships for HERG channel blockers.

The hERG potassium channel is of major pharmaceutical importance, and its blockade by various compounds, potentially causing serious cardiac side effects, is a major problem in drug development. Despite the large amounts of existing biochemical data on blockade of hERG by drugs and druglike compounds, relatively little is known regarding the structural basis of binding of blockers to the channel. Here, we have used a recently developed homology model of hERG to conduct molecular docking experiments with a series of channel blockers, followed by molecular dynamics simulations of the complexes and evaluation of binding free energies with the linear interaction energy method. The calculations yield a remarkably good agreement with experimental binding affinities and allow for a rationalization of three-dimensional structure-activity relationships in terms of a number of key interactions. Two main interaction regions of the channel are thus identified with implications for further mutagenesis experiments and design of new compounds.

Mutagenesis of human neuropeptide Y/peptide YY receptor Y2 reveals additional differences to Y1 in interactions with highly conserved ligand positions.

Neuropeptide Y (NPY) and peptide YY (PYY) share approximately 70% of their 36 amino acids and bind to the same three human receptor subtypes, Y1, Y2 and Y5, even though these receptors only share approximately 30% sequence identity. Based on our previous investigation of human Y1 we describe here a mutagenesis study of three corresponding positions in human Y2, i.e. Tyr2.64, Val6.58 and Tyr7.31. Pharmacological characterization was performed with the four peptide agonists PYY, NPY, PYY(3-36) and NPY(13-36) as well as the non-peptide antagonist BIIE0246. Results from mutants where Tyr2.64 has been substituted by Ala suggest that Tyr2.64 is involved in the interaction with all investigated ligands whereas position Tyr7.31 seems to be more important for interaction with the truncated peptide PYY(3-36) than with intact NPY. Surprisingly, substitution of Tyr7.31 with His, the corresponding residue in Y1, resulted in total loss of binding of iodinated porcine PYY. The third position, Val6.58, did not influence binding of any of the ligands. These findings differ from those obtained for Y1 where Ala substitution resulted in lost or changed binding for each of the three positions. Although Tyr2.64 and Tyr7.31 in Y2 are involved in ligand binding, their interactions with the peptide ligands seem to be different from the corresponding positions in Y1. This suggests that the receptor-ligand interactions have changed during evolution after Y1 and Y2 arose from a common ancestral receptor.

Toward a consensus model of the HERG potassium channel.

Malfunction of hERG potassium channels, due to inherited mutations or inhibition by drugs, can cause long QT syndrome, which can lead to life-threatening arrhythmias. A three-dimensional structure of hERG is a prerequisite to understand the molecular basis of hERG malfunction. To achieve a consensus model, we carried out an extensive analysis of hERG models based on various alignments of helix S5. We analyzed seven models using a combination of conventional geometry/packing/normality validation methods as well as molecular dynamics simulations and molecular docking. A synthetic test set with the X-ray crystal structure of K(v)1.2 with artificially shifted S5 sequences modeled into the structure served as a reference case. We docked the known hERG inhibitors (+)-cisapride, (S)-terfenadine, and MK-499 into the hERG models and simulation snapshots. None of the single analyses unambiguously identified a preferred model, but the combination of all three revealed that there is only one model that fulfils all quality criteria. This model is confirmed by a recent mutation scanning experiment (P. Ju, G. Pages, R. P. Riek, P. C. Chen, A. M. Torres, P. S. Bansal, S. Kuyucak, P. W. Kuchel, J. I. Vandenberg, J. Biol. Chem. 2009, 284, 1000-1008). We expect the modeled structure to be useful as a basis both for computational studies of channel function and kinetics as well as the design of experiments.

Mutagenesis and homology modeling of the Tn21 integron integrase IntI1.

Horizontal DNA transfer between bacteria is widespread and a major cause of antibiotic resistance. For logistic reasons, single or combined genes are shuttled between vectors such as plasmids and bacterial chromosomes. Special elements termed integrons operate in such shuttling and are therefore vital for horizontal gene transfer. Shorter elements carrying genes, cassettes, are integrated in the integrons, or excised from them, by virtue of a recombination site, attC, positioned in the 3' end of each unit. It is a remarkable and possibly restricting elementary feature of attC that it must be single-stranded while the partner target site, attI, may be double-stranded. The integron integrases belong to the tyrosine recombinase family, and this work reports mutations of the integrase IntI1 from transposon Tn21, chosen within a well-conserved region characteristic of the integron integrases. The mutated proteins were tested for binding to a bottom strand of an attC substrate, by using an electrophoresis mobility shift assay. To aid in interpreting the results, a homology model was constructed on the basis of the crystal structure of integron integrase VchIntIA from Vibrio cholerae bound to its cognate attC substrate VCRbs. The local stability and hydrogen bonding network of key domains of the modeled structure were further examined using molecular dynamics simulations. The homology model allowed us to interpret the roles of several amino acid residues, four of which were clearly binding assay responsive upon mutagenesis. Notably, we also observed features indicating that IntI1 may be more prone to base-specific contacts with VCRbs than VchIntIA.

Predicting binding modes from free energy calculations.

To produce reliable predictions of bioactive conformations is a major challenge in the field of structure-based inhibitor design and is a requirement for accurate binding free energy predictions with structure-based methods. A series of HIV-1 reverse transcriptase inhibitors was cross-docked using a non-native crystal structure that resulted in two distinct clusters of possible conformations. One of these clusters was compatible with an existing crystal structure, whereas the other displayed a flipped heterocyclic group. Binding free energies, using the non-native crystal structure, calculated from several scoring functions, were similar for the two clusters, and no conclusion about the binding mode could be drawn from these results. The two clusters could be separated through rescoring with the linear interaction method (LIE) in combination with molecular dynamics simulations, which leads to a binding mode prediction in line with experimental crystallographic data. Further, the LIE model produces the best correlation between experimental and calculated binding free energies among the tested scoring methods.

Combining docking, molecular dynamics and the linear interaction energy method to predict binding modes and affinities for non-nucleoside inhibitors to HIV-1 reverse transcriptase.

Docking, scoring, molecular dynamics (MD), and the linear interaction energy (LIE) method are used here to predict binding modes and affinities for a set of 43 non-nucleoside inhibitors to HIV-1 reverse transcriptase. Starting from a crystallographic structure, the binding modes of 43 inhibitors are predicted using automated docking. The Goldscore scoring function and the LIE method are then used to determine the relative binding free energies for the inhibitors. The Goldscore scoring function does not reproduce the relative binding affinities for the inhibitors, while the standard parametrization of the LIE method reproduces the experimental binding free energies for 39 inhibitors with an R (2) = 0.70 and an unsigned average error of 0.8 kcal/mol. The present calculations provide a validation of the combination of docking, MD, and LIE as a powerful tool in structure-based drug design, and the methodology is easily scalable for attaining a higher throughput of compounds.