Comprehensive and Automated Linear Interaction Energy Based Binding-Affinity Prediction for Multifarious Cytochrome P450 Aromatase Inhibitors.

Cytochrome P450 aromatase (CYP19A1) plays a key role in the development of estrogen dependent breast cancer, and aromatase inhibitors have been at the front line of treatment for the past three decades. The development of potent, selective and safer inhibitors is ongoing with in silico screening methods playing a more prominent role in the search for promising lead compounds in bioactivity-relevant chemical space. Here we present a set of comprehensive binding affinity prediction models for CYP19A1 using our automated Linear Interaction Energy (LIE) based workflow on a set of 132 putative and structurally diverse aromatase inhibitors obtained from a typical industrial screening study. We extended the workflow with machine learning methods to automatically cluster training and test compounds in order to maximize the number of explained compounds in one or more predictive LIE models. The method uses protein-ligand interaction profiles obtained from Molecular Dynamics (MD) trajectories to help model search and define the applicability domain of the resolved models. Our method was successful in accounting for 86% of the data set in 3 robust models that show high correlation between calculated and observed values for ligand-binding free energies (RMSE < 2.5 kJ mol-1), with good cross-validation statistics.

Prediction of ligand binding affinity and orientation of xenoestrogens to the estrogen receptor by molecular dynamics simulations and the linear interaction energy method.

Exposure to environmental estrogens has been proposed as a risk factor for disruption of reproductive development and tumorigenesis of humans and wildlife (McLachlan, J. A.; Korach, K. S.; Newbold, R. R.; Degen, G. H. Diethylstilbestrol and other estrogens in the environment. Fundam. Appl. Toxicol. 1984, 4, 686-691). In recent years, many structurally diverse environmental compounds have been identified as estrogens. A reliable computational method for determining estrogen receptor (ER) binding affinity is of great value for the prediction of estrogenic activity of such compounds and their metabolites. In the presented study, a computational model was developed for prediction of binding affinities of ligands to the ERalpha isoform, using MD simulations in combination with the linear interaction energy (LIE) approach. The linear interaction energy approximation was first described by Aqvist et al. (Aqvist, J.; Medina, C.; Samuelsson, J. E. A new method for predicting binding affinity in computer-aided drug design. Protein Eng. 1994, 7, 385-391) and relies on the assumption that the binding free energy (DeltaG) depends linearly on changes in the van der Waals and electrostatic energy of the system. In the present study, MD simulations of ligands in the ERalpha ligand binding domain (LBD) (Shiau, A. K.; Barstad, D.; Loria, P. M.; Cheng, L.; Kushner, P. J.; Agard, D. A.; Greene, G. L. The structural basis of estrogen receptor/coactivator recognition and the antagonism of this interaction by tamoxifen. Cell 1998, 95, 927-937), as well as ligands free in water, were carried out using the Amber 6.0 force field (http://amber.scripps.edu/). Contrary to previous LIE methods, we took into account every possible orientation of the ligands in the LBD and weighted the contribution of each orientation to the total binding affinity according to a Boltzman distribution. The training set (n = 19) contained estradiol (E2), the synthetic estrogens diethylstilbestrol (DES) and 11beta-chloroethylestradiol (E2-Cl), 16alpha-hydroxy-E2 (estriol, EST), the phytoestrogens genistein (GEN), 8-prenylnaringenin (8PN), and zearalenon (ZEA), four derivatives of benz[a]antracene-3,9-diol, and eight estrogenic monohydroxylated PAH metabolites. We obtained an excellent linear correlation (r(2) = 0.94) between experimental (competitive ER binding assay) and calculated binding energies, with K(d) values ranging from 0.15 mM to 30 pM, a 5 000 000-fold difference in binding affinity. Subsequently, a test set (n = 12) was used to examine the predictive value of our model. This set consisted of the synthetic estrogen 5,11-cis-diethyl-5,6,11,12-tetrahydrochrysene-2,8-diol (THC), daidzein (DAI), equol (EQU) and apigenin (API), chlordecone (KEP), progesterone (PRG), several mono- and dihydroxylated PAH metabolites, and two brominated biphenyls. The predicted binding affinities of these estrogenic compounds were in very good agreement with the experimental values (average deviation of 0.61 +/- 0.4 kcal/mol). In conclusion, our LIE model provides a very good method for prediction of absolute ligand binding affinities, as well as binding orientation of ligands.

Modeling and molecular dynamics of glutamine transaminase K/cysteine conjugate beta-lyase.

The homodimeric, pyridoxal 5'-phosphate (PLP)-dependent enzyme glutamine transaminase K/cysteine conjugate beta-lyase (GTK/beta-lyase) has been implicated in the bioactivation of chemopreventive compounds. This paper describes the first homology model of rat renal GTK/beta-lyase and its active site residues, deduced from molecular dynamics (MD) simulations of the binding mode of 13 structurally diverse cysteine S-conjugates and amino acids after Amber-parametrization of PLP. Comparison with Thermus thermophilus aspartate aminotransferase (tAAT) and Trypanosoma cruzi tyrosine aminotransferase (tTAT), used as templates for modeling GTK/beta-lyase, showed that the PLP-binding site of GTK/beta-lyase is highly conserved. Binding of the ligand alpha-carboxylate-group occurred via the conserved residues Arg(432) and Asn(219), and Asn(50) and Gly(70). Two pockets accommodated the various ligand side chains. A small pocket, located directly above PLP, was of a highly hydrophobic and aromatic character. A larger pocket, formed partly by the substrate access channel, was more hydrophilic and notably involved the salt bridge partners Glu(54) and Arg(99*) (* denotes the other subunit). Ligand-binding residues included Leu(51), Phe(71), Tyr(135), Phe(373) and Phe(312*), and pi-stacking interactions were often observed. Tyr(135) and Asn(50) were prominent in hydrogen bonding with the sulfur-atom of cysteine S-conjugates. The observed binding mode of the ligands corresponded well with their experimentally determined inhibitory potency toward GTK/beta-lyase. The current homology model thus provides a starting point for further validation of the role of active site residues in ligand-binding by means of mutagenesis studies. Ultimately, insight in the binding of ligands to GTK/beta-lyase may result in the rational design of new ligands and selective inhibitors.

Homology modeling of rat and human cytochrome P450 2D (CYP2D) isoforms and computational rationalization of experimental ligand-binding specificities.

The ligand-binding characteristics of rat and human CYP2D isoforms, i.e., rat CYP2D1-4 and human CYP2D6, were investigated by measuring IC(50) values of 11 known CYP2D6 ligands using 7-methoxy-4-(aminomethyl)coumarin (MAMC) as substrate. Like CYP2D6, all rat CYP2D isozymes catalyzed the O-demethylation of MAMC with K(m) and V(max) values ranging between 78 and 145 microM and 0.048 and 1.122 min(-1), respectively. To rationalize observed differences in the experimentally determined IC(50) values, homology models of the CYP2D isoforms were constructed. A homology model of CYP2D6 was generated on the basis of crystallized rabbit CYP2C5 and was validated on its ability to reproduce binding orientations corresponding to metabolic profiles of the substrates and to remain stable during unrestrained molecular dynamics simulations at 300 K. Twenty-two active site residues, sharing up to 59% sequence identity, were identified in the CYP2D binding pockets and included CYP2D6 residues Phe120, Glu216, and Asp301. Electrostatic potential calculations displayed large differences in the negative charge of the CYP2D active sites, which was consistent with observed differences in absolute IC(50) values. MD studies on the binding mode of sparteine, quinidine, and quinine in CYP2D2 and CYP2D6 furthermore concurred well with experimentally determined IC(50) values and metabolic profiles. The current study thus provides new insights into differences in the active site topology of the investigated CYP2D isoforms.