Methodology for protein-ligand binding studies: application to a model for drug resistance, the HIV/FIV protease system

Template-Based Method for Conformation Generation and Scoring for Congeneric Series of Ligands

Physics-based prediction of protein-ligand binding affinities for a congeneric series of ligands in lead optimization requires their geometries as a first step. In this paper, we report a method that uses the 3D conformation of a lead compound in complex with a protein as a template to generate conformations of a series of related analog compounds. The method uses the Maximal Common Substructure (MCS) computed between lead and analog ligands to assign coordinates for the atoms shared between the ligands. For the differing atoms, a conformation generation procedure is implemented that results in a diversity of conformations. The generated conformations are sorted using a score based on the Molecular Mechanics and Generalized Born with Solvent Accessible Surface Area contribution (MM-GBSA) method. The accuracy of the generated conformations is tested retrospectively using a cross-validation approach applied to four data sets obtained from the Drug Design Data Resource (D3R) by measuring the RMSD of the top scored conformation with respect to the crystallographic pose. The scoring ability of the method is independently assessed using data for the same protein targets to test the rank ordering ability and separating active and inactive ligands. We tested the effect of protein flexibility during structural optimization and scoring approaches with and without strain energies. Retrospective validation on data sets comprising 4 targets shows that the method outperforms random selection for all targets and outperforms a molecular weight-based null model in 3 out of 4 targets in separating active and inactive compounds. Therefore, the presented method is expected to be of utility in lead optimization for rapidly screening analog ligands and generating initial conformations for use in more detailed physics-based binding affinity prediction methods.

Effects from metal ion in tumor endothelial marker 8 and anthrax protective antigen: BioLayer Interferometry experiment and molecular dynamics simulation study

One of the anthrax receptors, tumor endothelial marker 8 (TEM8), is reported to be a potential anticancer target due to its over-expression during tumor angiogenesis. To extend our BioLayer Interferometry study in PA-TEM8 binding, we present a computational approach to reveal the role of an integral metal ion on receptor structure and binding thermodynamics. We estimated the interaction energy between PA and TEM8 using computer simulation. Consistent with experimental study, computational results indicate the metal ion in TEM8 contributes significantly to the binding affinity, and PA-TEM8 binding is more favorable in the presence of Mg²⁺ than Ca²⁺ . Further, computational analysis suggests that the differences in PA-TEM8 binding affinity are comparable to the closely related integrin proteins. The conformation change, which linked to changes in activity of integrins, was not found in TEM8. In the present of Mg²⁺ , TEM8 remains in a conformation analogous to an integrin open (high-affinity) conformation. © 2017 Wiley Periodicals, Inc.

Rehabilitating drug-induced long-QT promoters: in-silico design of hERG-neutral cisapride analogues with retained pharmacological activity

Background

The human ether-a-go-go related gene 1 (hERG1), which codes for a potassium ion channel, is a key element in the cardiac delayed rectified potassium current, I_Kr, and plays an important role in the normal repolarization of the heart’s action potential. Many approved drugs have been withdrawn from the market due to their prolongation of the QT interval. Most of these drugs have high potencies for their principal targets and are often irreplaceable, thus “rehabilitation” studies for decreasing their high hERG1 blocking affinities, while keeping them active at the binding sites of their targets, have been proposed to enable these drugs to re-enter the market.

Methods

In this proof-of-principle study, we focus on cisapride, a gastroprokinetic agent withdrawn from the market due to its high hERG1 blocking affinity. Here we tested an a priori strategy to predict a compound’s cardiotoxicity using de novo drug design with molecular docking and Molecular Dynamics (MD) simulations to generate a strategy for the rehabilitation of cisapride.

Results

We focused on two key receptors, a target interaction with the (adenosine) receptor and an off-target interaction with hERG1 channels. An analysis of the fragment interactions of cisapride at human A_2A adenosine receptors and hERG1 central cavities helped us to identify the key chemical groups responsible for the drug activity and hERG1 blockade. A set of cisapride derivatives with reduced cardiotoxicity was then proposed using an in-silico two-tier approach. This set was compared against a large dataset of commercially available cisapride analogs and derivatives.

Conclusions

An interaction decomposition of cisapride and cisapride derivatives allowed for the identification of key active scaffolds and functional groups that may be responsible for the unwanted blockade of hERG1.

Ensemble QSAR: a QSAR method based on conformational ensembles and metric descriptors

Quantitative structure-activity relationship (QSAR) is the most versatile tool in computer-assisted molecular design. One conceptual drawback seen in QSAR approaches is the "one chemical-one structure-one parameter value" dogma where the model development is based on physicochemical description for a single molecular conformation, while ignoring the rest of the conformational space. It is well known that molecules have several low-energy conformations populated at physiological temperature, and each conformer makes a significant impact on associated properties such as biological activity. At the level of molecular interaction, the dynamics around the molecular structure is of prime essence rather than the average structure. As a step toward understanding the role of these discrete microscopic states in biological activity, we have put together a theoretically rigorous and computationally tractable formalism coined as eQSAR. In this approach, the biological activity is modeled as a function of physicochemical description for a selected set of low-energy conformers, rather than that's for a single lowest energy conformation. Eigenvalues derived from the "Physicochemical property integrated distance matrices" (PD-matrices) that encompass both 3D structure and physicochemical properties, have been used as descriptors; is a novel addition. eQSAR is validated on three peptide datasets and explicitly elaborated for bradykinin-potentiating peptides. The conformational ensembles were generated by a simple molecular dynamics and consensus dynamics approaches. The eQSAR models are statistically significant and possess the ability to select the most biologically relevant conformation(s) with the relevant physicochemical attributes that have the greatest meaning for description of the biological activity.

NMR and MD investigations of human galectin-1/oligosaccharide complexes

The specific recognition of carbohydrates by lectins plays a major role in many cellular processes. Galectin-1 belongs to a family of 15 structurally related beta-galactoside binding proteins that are able to control a variety of cellular events, including cell cycle regulation, adhesion, proliferation, and apoptosis. The three-dimensional structure of galectin-1 has been solved by x-ray crystallography in the free form and in complex with various carbohydrate ligands. In this work, we used a combination of two-dimensional NMR titration experiments and molecular-dynamics simulations with explicit solvent to study the mode of interaction between human galectin-1 and five galactose-containing ligands. Isothermal titration calorimetry measurements were performed to determine their affinities for galectin-1. The contribution of the different hexopyranose units in the protein-carbohydrate interaction was given particular consideration. Although the galactose moiety of each oligosaccharide is necessary for binding, it is not sufficient by itself. The nature of both the reducing sugar in the disaccharide and the interglycosidic linkage play essential roles in the binding to human galectin-1.

Structural mining: self-consistent design on flexible protein-peptide docking and transferable binding affinity potential

A flexible protein-peptide docking method has been designed to consider not only ligand flexibility but also the flexibility of the protein. The method is based on a Monte Carlo annealing process. Simulations with a distance root-mean-square (dRMS) virtual energy function revealed that the flexibility of protein side chains was as important as ligand flexibility for successful protein-peptide docking. On the basis of mean field theory, a transferable potential was designed to evaluate distance-dependent protein-ligand interactions and atomic solvation energies. The potential parameters were developed using a self-consistent process based on only 10 known complex structures. The effectiveness of each intermediate potential was judged on the basis of a Z score, approximating the gap between the energy of the native complex and the average energy of a decoy set. The Z score was determined using experimentally determined native structures and decoys generated by docking with the intermediate potentials. Using 6600 generated decoys and the Z score optimization criterion proposed in this work, the developed potential yielded an acceptable correlation of R(2) = 0.77, with binding free energies determined for known MHC I complexes (Class I Major Histocompatibility protein HLA-A(*)0201) which were not present in the training set. Test docking on 25 complexes further revealed a significant correlation between energy and dRMS, important for identifying native-like conformations. The near-native structures always belonged to one of the conformational classes with lower predicted binding energy. The lowest energy docked conformations are generally associated with near-native conformations, less than 3.0 Angstrom dRMS (and in many cases less than 1.0 Angstrom) from the experimentally determined structures.

Structure-based phenotyping predicts HIV-1 protease inhibitor resistance

Mutations in HIV-1 drug targets lead to resistance and consequent therapeutic failure of antiretroviral drugs. Phenotypic resistance assays are time-consuming and costly, and genotypic rules-based interpretations may fail to predict the effects of multiple mutations. We have developed a computational procedure that rapidly evaluates changes in the binding energy of inhibitors to mutant HIV-1 PR variants. Models of WT complexes were produced from crystal structures. Mutant complexes were built by amino acid substitutions in the WT complexes with subsequent energy minimization of the ligand and PR binding site residues. Accuracy of the models was confirmed by comparison with available crystal structures and by prediction of known resistance-related mutations. PR variants from clinical isolates were modeled in complex with six FDA-approved PIs, and changes in the binding energy (DeltaE(bind)) of mutant versus WT complexes were correlated with the ratios of phenotypic 50% inhibitory concentration (IC(50)) values. The calculated DeltaE(bind) of five PIs showed significant correlations (R(2) = 0.7-0.8) with IC(50) ratios from the Virco Antivirogram assay, and the DeltaE(bind) of six PIs showed good correlation (R(2) = 0.76-0.85) with IC(50) ratios from the Virologic PhenoSense assay. DeltaE(bind) cutoffs corresponding to a four-fold increase in IC(50) were used to define the structure-based phenotype as susceptible, resistant, or equivocal. Blind predictions for 78 PR variants gave overall agreement of 92% (kappa = 0.756) and 86% (kappa = 0.666) with PhenoSense and Antivirogram phenotypes, respectively. The structural phenotyping predicted drug resistance of clinical HIV-1 PR variants with an accuracy approaching that of frequently used cell-based phenotypic assays.

Structural basis for distinctions between substrate and inhibitor specificities for feline immunodeficiency virus and human immunodeficiency virus proteases

We used feline immunodeficiency virus (FIV) protease (PR) as a mutational framework to define determinants for the observed substrate and inhibitor specificity distinctions between FIV and human immunodeficiency virus (HIV) PRs. Multiple-substitution mutants were constructed by replacing the residues in and around the active site of FIV PR with the structurally equivalent residues of HIV-1 PR. Mutants included combinations of three critical regions (FIV numbering, with equivalent HIV numbering in superscript): I37(32)V in the active core region; N55(46)M, M56(47)I, and V59(50)I in the flap region; and L97(80)T, I98(81)P, Q99(82)V, P100(83)N, and L101(84)I in the 90s loop region. Significant alterations in specificity were observed, consistent with the involvement of these residues in determining the substrate-inhibitor specificity distinctions between FIV and HIV PRs. Two previously identified residues, I35 and I57 of FIV PR, were intolerant to substitution and yielded inactive PRs. Therefore, we attempted to recover the activity by introducing secondary mutations. The addition of G62(53)F and K63(54)I, located at the top of the flap and outside the active site, compensated for the activity lost in the I57(48)G substitution mutants. An additional two substitutions, D105(88)N and N88(74)T, facilitated recovery of activity in mutants that included the I35(30)D substitution. Determination of K(i) values of potent HIV-1 PR inhibitors against these mutants showed that inhibitor specificity paralleled that of HIV-1 PR. The findings indicate that maintenance of both substrate and inhibitor specificity is a function of interactions between residues both inside and outside the active site. Thus, mutations apparently peripheral to the active site can have a dramatic influence on inhibitor efficacy.

Elucidating the inhibiting mode of AHPBA derivatives against HIV-1 protease and building predictive 3D-QSAR models

The Lamarckian genetic algorithm of AutoDock 3.0 has been used to dock 27 3(S)-amino-2(S)-hydroxyl-4-phenylbutanoic acids (AHPBAs) into the active site of HIV-1 protease (HIVPR). The binding mode was demonstrated in the aspects of the inhibitor's conformation, subsite interaction, and hydrogen bonding. The data of geometrical parameters (tau(1), tau(2), and tau(3) listed in Table 2) and root mean square deviation values as compared with the known inhibitor, kni272,(28) show that both kinds of inhibitors interact with HIVPR in a very similar way. The r(2) value of 0.860 indicates that the calculated binding free energies correlate well with the inhibitory activities. The structural and energetic differences in inhibitory potencies of AHPBAs were reasonably explored. Using the binding conformations of AHPBAs, consistent and highly predictive 3D-QSAR models were developed by performing CoMFA, CoMSIA, and HQSAR analyses. The reasonable r(corss)(2) values were 0.613, 0.530, and 0.717 for CoMFA, CoMSIA, and HQSAR models, respectively. The predictive ability of these models was validated by kni272 and a set of nine compounds that were not included in the training set. Mapping these models back to the topology of the active site of HIVPR leads to a better understanding of vital AHPBA-HIVPR interactions. Structural-based investigations and the final 3D-QSAR results provide clear guidelines and accurate activity predictions for novel HIVPR inhibitors.

Macromolecular electrostatics: continuum models and their growing pains

Theoretical understanding of macromolecular electrostatics has advanced substantially over the past year. Continuum models have given promising results for calculating protein-ligand binding free energy differences, as well as pK(a)s and redox properties, particularly with explicit treatment of multiple conformers. Generalized Born and other techniques have led to the first molecular dynamics simulations of proteins and RNA with continuum solvent. Continuum and microscopic descriptions of dielectric relaxation have been critically compared.

Estimation of binding affinities for selective thrombin inhibitors via Monte Carlo simulations

Monte Carlo simulations have been performed on a series of 20 active-site-directed thrombin inhibitors to determine the interactions and energetics associated with the binding of these compounds. Physicochemical descriptors of potential value in the prediction of binding affinities were averaged during simulations of each inhibitor unbound in water and bound to thrombin. Regression equations based on 3-5 descriptors are able to reproduce the experimental binding affinities, which cover a 7 kcal/mol range, with rms errors of 1.0-1.3 kcal/mol, and yield correlation coefficients, r(2), of 0.7-0.8. On the basis of these results, the quantities most important in determining the binding affinities are: (1) the enhancement of van der Waals interactions in going from solution to the bound state, (2) the intramolecular strain induced in the inhibitor upon binding, (3) the number of hydrogen bonds lost in the binding process, and (4) the number of rotatable bonds in the inhibitor. The descriptors are physically reasonable and, in combination with the insights gained from analysis of the simulation structures, suggest directions for the development of improved thrombin inhibitors.

Alteration of substrate and inhibitor specificity of feline immunodeficiency virus protease

Feline immunodeficiency virus (FIV) protease is structurally very similar to human immunodeficiency virus (HIV) protease but exhibits distinct substrate and inhibitor specificities. We performed mutagenesis of subsite residues of FIV protease in order to define interactions that dictate this specificity. The I37V, N55M, M56I, V59I, and Q99V mutants yielded full activity. The I37V, N55M, V59I, and Q99V mutants showed a significant increase in activity against the HIV-1 reverse transcriptase/integrase and P2/nucleocapsid junction peptides compared with wild-type (wt) FIV protease. The I37V, V59I, and Q99V mutants also showed an increase in activity against two rapidly cleaved peptides selected by cleavage of a phage display library with HIV-1 protease. Mutations at Q54K, I98P, and L101I dramatically reduced activity. Mutants containing a I35D or I57G substitution showed no activity against either FIV or HIV substrates. FIV proteases all failed to cut HIV-1 matrix/capsid, P1/P6, P6/protease, and protease/reverse transcriptase junctions, indicating that none of the substitutions were sufficient to change the specificity completely. The I37V, N55M, M56I, V59I, and Q99V mutants, compared with wt FIV protease, all showed inhibitor specificity more similar to that of HIV-1 protease. The data also suggest that FIV protease prefers a hydrophobic P2/P2' residue like Val over Asn or Glu, which are utilized by HIV-1 protease, and that S2/S2' might play a critical role in distinguishing FIV and HIV-1 protease by specificity. The findings extend our observations regarding the interactions involved in substrate binding and aid in the development of broad-based inhibitors.