Prediction of Accurate Binding Modes Using Combination of Classical and Accelerated Molecular Dynamics and Free-Energy Perturbation Calculations: An Application to Toxicity Studies.

Estimating the correct binding modes of ligands in protein-ligand complexes is crucial not only in the drug discovery process but also for elucidating potential toxicity mechanisms. In the current paper, we propose a computational modeling workflow using the combination of docking, classical molecular dynamics (cMD), accelerated molecular dynamics (aMD) and free-energy perturbation (FEP+ protocol) for identification of possible ligand binding modes. It was applied for investigation of selected perfluorocarboxyl acids (PFCAs) in the PPARγ nuclear receptor. Although both regular and induced fit docking failed to reproduce the experimentally determined binding mode of the ligands when docked into a non-native X-ray structure, cMD and aMD simulations successfully identified the most probable binding conformations. Moreover, multiple binding modes were identified for all of these compounds and the shorter-chain PFCAs continuously moved between a few energetically favorable binding conformations. On the basis of MD predictions of binding conformations, we applied the default and also redesigned FEP+ sampling protocols, which accurately reproduced experimental differences in the binding energies. Thus, the preliminary MD simulations can also provide helpful information about correct setup of the FEP+ calculations. These results show that the PFCA binding modes were accurately predicted and that the FEP+ protocol can be used to estimate free energies of binding of flexible ligands that are not typical druglike compounds. Our in silico workflow revealed the specific ligand-residue interactions within the ligand binding domain and the main characteristics of the PFCAs, and it was concluded that these compounds are week PPARγ partial agonists. This work also suggests a common pipeline for identification of ligand binding modes, ligand-protein dynamics description, and relative free-energy calculations.

Physicologically Based Toxicokinetic Models of Tebuconazole and Application in Human Risk Assessment.

A series of physiologically based toxicokinetic (PBTK) models for tebuconazole were developed in four species, rat, rabbit, rhesus monkey, and human. The developed models were analyzed with respect to the application of the models in higher tier human risk assessment, and the prospect of using such models in risk assessment of cumulative and aggregate exposure is discussed. Relatively simple and biologically sound models were developed using available experimental data as parameters for describing the physiology of the species, as well as the absorption, distribution, metabolism, and elimination (ADME) of tebuconazole. The developed models were validated on in vivo half-life data for rabbit with good results, and on plasma and tissue concentration-time course data of tebuconazole after i.v. administration in rabbit. In most cases, the predicted concentration levels were seen to be within a factor of 2 compared to the experimental data, which is the threshold set for the use of PBTK simulation results in risk assessment. An exception to this was seen for one of the target organs, namely, the liver, for which tebuconazole concentration was significantly underestimated, a trend also seen in model simulations for the liver after other nonoral exposure scenarios. Possible reasons for this are discussed in the article. Realistic dietary and dermal exposure scenarios were derived based on available exposure estimates, and the human version of the PBTK model was used to simulate the internal levels of tebuconazole and metabolites in the human body for these scenarios. By a variant of the models where the R(-)- and S(+)-enantiomers were treated as two components in a binary mixture, it was illustrated that the inhibition between the two tebuconazole enantiomers did not affect the simulation results for these realistic exposure scenarios. The developed models have potential as an important tool in risk assessment.

Discovery of a novel selective PPARγ ligand with partial agonist binding properties by integrated in silico/in vitro work flow.

Full agonists to the peroxisome proliferator-activated receptor (PPAR)γ, such as Rosiglitazone, have been associated with a series of undesired side effects, such as weight gain, fluid retention, cardiac hypertrophy, and hepatotoxicity. Nevertheless, PPARγ is involved in the expression of genes that control glucose and lipid metabolism and is an important target for drugs against type 2 diabetes, dyslipidemia, atherosclerosis, and cardiovascular disease. In an effort to identify novel PPARγ ligands with an improved pharmacological profile, emphasis has shifted to selective ligands with partial agonist binding properties. Toward this end we applied an integrated in silico/in vitro workflow, based on pharmacophore- and structure-based virtual screening of the ZINC library, coupled with competitive binding and transactivation assays, and adipocyte differentiation and gene expression studies. Hit compound 9 was identified as the most potent ligand (IC50 = 0.3 µM) and a relatively poor inducer of adipocyte differentiation. The binding mode of compound 9 was confirmed by molecular dynamics simulation, and the calculated free energy of binding was -8.4 kcal/mol. A novel functional group, the carbonitrile group, was identified to be a key substituent in the ligand-protein interactions. Further studies on the transcriptional regulation properties of compound 9 revealed a gene regulatory profile that was to a large extent unique, however functionally closer to that of a partial agonist.

QSAR model for human pregnane X receptor (PXR) binding: screening of environmental chemicals and correlations with genotoxicity, endocrine disruption and teratogenicity.

The pregnane X receptor (PXR) has a key role in regulating the metabolism and transport of structurally diverse endogenous and exogenous compounds. Activation of PXR has the potential to initiate adverse effects, causing drug-drug interactions, and perturbing normal physiological functions. Therefore, identification of PXR ligands would be valuable information for pharmaceutical and toxicological research. In the present study, we developed a quantitative structure-activity relationship (QSAR) model for the identification of PXR ligands using data based on a human PXR binding assay. A total of 631 molecules, representing a variety of chemical structures, constituted the training set of the model. Cross-validation of the model showed a sensitivity of 82%, a specificity of 85%, and a concordance of 84%. The developed model provided knowledge about molecular descriptors that may influence the binding of molecules to PXR. The model was used to screen a large inventory of environmental chemicals, of which 47% was found to be within domain of the model. Approximately 35% of the chemicals within domain were predicted to be PXR ligands. The predicted PXR ligands were found to be overrepresented among chemicals predicted to cause adverse effects, such as genotoxicity, teratogenicity, estrogen receptor activation and androgen receptor antagonism compared to chemicals not causing these effects. The developed model may be useful as a tool for predicting potential PXR ligands and for providing mechanistic information of toxic effects of chemicals.

Identification of cytochrome P450 2D6 and 2C9 substrates and inhibitors by QSAR analysis.

This paper presents four new QSAR models for CYP2C9 and CYP2D6 substrate recognition and inhibitor identification based on human clinical data. The models were used to screen a large data set of environmental chemicals for CYP activity, and to analyze the frequency of CYP activity among these compounds. A large fraction of these chemicals were found to be CYP active, and thus potentially capable of affecting human physiology. 20% of the compounds within applicability domain of the models were predicted to be CYP2C9 substrates, and 17% to be inhibitors. The corresponding numbers for CYP2D6 were 9% and 21%. Where the majority of CYP2C9 active compounds were predicted to be both a substrate and an inhibitor at the same time, the CYP2D6 active compounds were primarily predicted to be only inhibitors. It was demonstrated that the models could identify compound classes with a high occurrence of specific CYP activity. An overrepresentation was seen for poly-aromatic hydrocarbons (group of procarcinogens) among CYP2C9 active and mutagenic compounds compared to CYP2C9 inactive and mutagenic compounds. The mutagenicity was predicted with a QSAR model based on Ames in vitro test data.

Prediction of drug efficacy for cancer treatment based on comparative analysis of chemosensitivity and gene expression data.

The NCI60 database is the largest available collection of compounds with measured anti-cancer activity. The strengths and limitations for using the NCI60 database as a source of new anti-cancer agents are explored and discussed in relation to previous studies. We selected a sub-set of 2333 compounds with reliable experimental half maximum growth inhibitions (GI(50)) values for 30 cell lines from the NCI60 data set and evaluated their growth inhibitory effect (chemosensitivity) with respect to tissue of origin. This was done by identifying natural clusters in the chemosensitivity data set and in a data set of expression profiles of 1901 genes for the corresponding tumor cell lines. Five clusters were identified based on the gene expression data using self-organizing maps (SOM), comprising leukemia, melanoma, ovarian and prostate, basal breast, and luminal breast cancer cells, respectively. The strong difference in gene expression between basal and luminal breast cancer cells was reflected clearly in the chemosensitivity data. Although most compounds in the data set were of low potency, high efficacy compounds that showed specificity with respect to tissue of origin could be found. Furthermore, eight potential topoisomerase II inhibitors were identified using a structural similarity search. Finally, a set of genes with expression profiles that were significantly correlated with anti-cancer drug activity was identified. Our study demonstrates that the combined data sets, which provide comprehensive information on drug activity and gene expression profiles of tumor cell lines studied, are useful for identifying potential new active compounds.

The phosphorylation specificity of B-RAF WT, B-RAF D594V, B-RAF V600E and B-RAF K601E kinases: an in silico study.

Phosphorylation of the B-RAF kinase is an essential process in tumour induction and maintenance in several cancers. Herein the phosphorylation specificity of the activation segment of the wild type B-RAF kinase and the B-RAF(D594V), B-RAF(V600E) and B-RAF(K601E) mutants was examined by molecular dynamics (MD) simulations and GRID molecular interaction field analysis. According to our analysis, Thr599 and Ser602 were the only residues in the activation segment in B-RAF(WT) that were well exposed to ATP binding, which is in agreement with the experimental results, and provide a molecular basis of the observed phosphorylation. The phosphorylation specificity was altered significantly for each of the three different mutants studied due to the large conformational changes and subsequent alterations in the electrostatic forces between several residues for each of these mutants. Thus the analysis revealed limited phosphorylation potential of the non-active B-RAF(D594V) mutant and several potential ATP binding sites were identified for the highly active B-RAF(V600E) mutant. The Lys601 residue, which is specific to RAF and not present in the activation segment of other similar kinases, was identified to potentially be of major importance to the observed differences in the phosphorylation specificity of the mutants. Our results indicate that Lys601 might be a specific ATP coordinating residue, contributing to the B-RAF phosphorylation specificity. The overall results can be helpful for the understanding of the B-RAF phosphorylation processes on a molecular level.

An in silico study of the molecular basis of B-RAF activation and conformational stability.

BACKGROUND: B-RAF kinase plays an important role both in tumour induction and maintenance in several cancers and it is an attractive new drug target. However, the structural basis of the B-RAF activation is still not well understood. RESULTS: In this study we suggest a novel molecular basis of B-RAF activation based on molecular dynamics (MD) simulations of B-RAFWT and the B-RAFV600E, B-RAFK601E and B-RAFD594V mutants. A strong hydrogen bond network was identified in B-RAFWT in which the interactions between Lys601 and the well known catalytic residues Lys483, Glu501 and Asp594 play an important role. It was found that several mutations, which directly or indirectly destabilized the interactions between these residues within this network, contributed to the changes in B-RAF activity. CONCLUSION: Our results showed that the above mechanisms lead to the disruption of the electrostatic interactions between the A-loop and the alphaC-helix in the activating mutants, which presumably contribute to the flipping of the activation segment to an active form. Conversely, in the B-RAFD594V mutant that has impaired kinase activity, and in B-RAFWT these interactions were strong and stabilized the kinase inactive form.

Molecular basis of inactive B-RAF(WT) and B-RAF(V600E) ligand inhibition, selectivity and conformational stability: an in silico study.

The B-RAF kinase plays an important role both in tumor induction and maintenance in several cancers. The molecular basis of the inactive B-RAF(WT) and B-RAF(V600E) inhibition and selectivity of a series of inhibitors was examined with a combination of molecular dynamics (MD), free energy MM-PBSA and local-binding energy (LBE) approaches. The conformational stability of the unbounded kinases and in particular the processes of the B-RAF (V600E) mutant activation were analyzed. A unique salt bridge network formed mainly by the catalytic residues was identified in the unbounded B-RAFs. The reorganization of this network and the restriction of the active segment flexibility upon ligand binding inhibit both B-RAF(WT) and B-RAF (V600E), thus appearing as an important factor for ligand selectivity. A significant correlation between the binding energies of the compounds in B-RAF(WT) and their inhibition effects on B-RAF (V600E) was revealed, which can explain the low mutant selectivity observed for numerous inhibitors. Our results suggest that the interactions between the activation segment and the alpha C-helix, as well as between the residues in the salt bridge network, are the major mechanism of the B-RAF (V600E) activation. Overall data revealed the important role of Lys601 for ligand activity, selectivity and protein stabilization, proposing an explanation of the observed strong kinase activation in the K601E mutated form.

Prediction of pH-dependent aqueous solubility of druglike molecules.

In the present work, the Henderson-Hasselbalch (HH) equation has been employed for the development of a tool for the prediction of pH-dependent aqueous solubility of drugs and drug candidates. A new prediction method for the intrinsic solubility was developed, based on artificial neural networks that have been trained on a druglike PHYSPROP subset of 4548 compounds. For the prediction of acid/base dissociation coefficients, the commercial tool Marvin has been used, following validation on a data set of 467 molecules from the PHYSPROP database. The best performing network for intrinsic solubility predictions has a cross-validated root mean square error (RMSE) of 0.70 log S-units, while the Marvin pKa plug-in has an RMSE of 0.71 pH-units. A data set of 27 drugs with experimentally determined pH-solubility curves was assembled from the literature for the validation of the combined pH-dependent model, giving a mean RMSE of 0.79 log S-units. Finally, the combined model has been applied on profiling the solubility space at low pH of five large vendor libraries.

Prediction methods and databases within chemoinformatics: emphasis on drugs and drug candidates.

MOTIVATION: To gather information about available databases and chemoinformatics methods for prediction of properties relevant to the drug discovery and optimization process. RESULTS: We present an overview of the most important databases with 2-dimensional and 3-dimensional structural information about drugs and drug candidates, and of databases with relevant properties. Access to experimental data and numerical methods for selecting and utilizing these data is crucial for developing accurate predictive in silico models. Many interesting predictive methods for classifying the suitability of chemical compounds as potential drugs, as well as for predicting their physico-chemical and ADMET properties have been proposed in recent years. These methods are discussed, and some possible future directions in this rapidly developing field are described.

QSPR models for various physical properties of carbohydrates based on molecular mechanics and quantum chemical calculations.

Quantitative Structure-Property Relationships (QSPR) have been developed for a series of monosaccharides, including the physical properties of partial molar heat capacity, heat of solution, melting point, heat of fusion, glass-transition temperature, and solid state density. The models were based on molecular descriptors obtained from molecular mechanics and quantum chemical calculations, combined with other types of descriptors. Saccharides exhibit a large degree of conformational flexibility, therefore a methodology for selecting the energetically most favorable conformers has been developed, and was used for the development of the QSPR models. In most cases good correlations were obtained for monosaccharides. For five of the properties predictions were made for disaccharides, and the predicted values for the partial molar heat capacities were in excellent agreement with experimental values.

Modeling of the phase equilibria of polystyrene in methylcyclohexane with semi-empirical quantum mechanical methods I.

A method for calculating interaction parameters traditionally used in phase-equilibrium computations in low-molecular systems has been extended for the prediction of solvent activities of aromatic polymer solutions (polystyrene+methylcyclohexane). Using ethylbenzene as a model compound for the repeating unit of the polymer, the intermolecular interaction energies between the solvent molecule and the polymer were simulated. The semiempirical quantum chemical method AM1, and a method for sampling relevant internal orientations for a pair of molecules developed previously were used. Interaction energies are determined for three molecular pairs, the solvent and the model molecule, two solvent molecules and two model molecules, and used to calculated UNIQUAC interaction parameters, a(ij) and a(ji). Using these parameters, the solvent activities of the polystyrene 90,000 amu+methylcyclohexane system, and the total vapor pressures of the methylcyclohexane+ethylbenzene system were calculated. The latter system was compared to experimental data, giving qualitative agreement. Figure Solvent activities for the methylcylcohexane(1)+polystyrene(2) system at 316 K. Parameters aij (blue line) obtained with the AM1 method; parameters aij (pink line) from VLE data for the ethylbenzene+methylcyclohexane system. The abscissa is the polymer weight fraction defined as y2(x1)=(1mx1)M2/[x1M1+(1mx1)M2], where x1 is the solvent mole fraction and Mi are the molecular weights of the components.

Modeling and measurements of solid-liquid and vapor-liquid equilibria of polyols and carbohydrates in aqueous solution.

The solubilities of five saccharides in water have been measured at various temperatures. This includes the monosaccharides xylose and galactose, and the disaccharides maltose monohydrate, cellobiose and trehalose dihydrate. A method that uses interaction energies and interaction parameters calculated with molecular mechanics methods has shown to give good predictions of the phase behavior of a variety of mixtures, including glycols and small saccharides in aqueous solution. The method is completely predictive, as the strength of the molecular interactions is determined with a theoretical method in the absence of any phase equilibrium data. For calculating solubilities, experimental values for the melting points and the heats of fusion of the compounds under study are, however, necessary. The solubilities of the five saccharides listed above, raffinose and meso-erythritol in water were calculated with this method. The calculated solubilities are in reasonably good agreement with experiment, and in the case of meso-erythritol, which is a polyalcohol (polyol), and galactose, the agreement between prediction and experiment is excellent. Also the vapor pressures of water over several polyols and saccharides in aqueous solution have been predicted with this method, giving results in excellent agreement with the experimental values.