Benzimidazoles Containing Piperazine Skeleton at C-2 Position as Promising Tubulin Modulators with Anthelmintic and Antineoplastic Activity.

Benzimidazole anthelmintic drugs hold promise for repurposing as cancer treatments due to their interference with tubulin polymerization and depolymerization, manifesting anticancer properties. We explored the potential of benzimidazole compounds with a piperazine fragment at C-2 as tubulin-targeting agents. In particular, we assessed their anthelmintic activity against isolated Trichinella spiralis muscle larvae and their effects on glioblastoma (U-87 MG) and breast cancer (MDA-MB-231) cell lines. Compound 7c demonstrated exceptional anthelmintic efficacy, achieving a 92.7% reduction in parasite activity at 100 µg/mL after 48 hours. In vitro cytotoxicity analysis of MDA-MB 231 and U87 MG cell lines showed that derivatives 7b, 7d, and 7c displayed lower IC50 values compared to albendazole (ABZ), the control. These piperazine benzimidazoles effectively reduced cell migration in both cell lines, with compound 7c exhibiting the most significant reduction, making it a promising candidate for further study. The binding mode of the most promising compound 7c, was determined using the induced fit docking-molecular dynamics (IFD-MD) approach. Regular docking and IFD were also employed for comparison. The IFD-MD analysis revealed that 7c binds to tubulin in a unique binding cavity near that of ABZ, but the benzimidazole ring was fitted much deeper into the binding pocket. Finally, the absolute free energy of perturbation technique was applied to evaluate the 7c binding affinity, further confirming the observed binding mode.

R346K Mutation in the Mu Variant of SARS-CoV-2 Alters the Interactions with Monoclonal Antibodies from Class 2: A Free Energy Perturbation Study.

The Mu variant of SARS-CoV-2 has been recently classified as a variant of interest (VOI) by the World Health Organization (WHO) but limited data are available at the moment. In particular, special attention was given to the R346K mutation located in the receptor binding domain (RBD). In the current study, we performed free energy perturbation (FEP) calculations to elucidate its possible impact on a set of neutralizing monoclonal antibodies (mAbs) that have been shown to be strong inhibitors of the most other known COVID-19 variants. Our results show that R346K affects class 2 antibodies but its effect is not so significant (0.66 kcal/mol), i.e., it reduces the binding with antibodies by about 3-fold. An identical value was also calculated in the presence of both class 1 and class 2 antibodies (BD-812/836). Further, a similar reduction in the binding (0.4 kcal/mol) was obtained for the BD-821/771 pair of mAbs. For comparison, the addition of the K417N mutation, present in the newly registered Mu variant in July 2021 in the U.K., affected the class 1 mAbs by strongly reducing the binding by 1.29 kcal/mol or about 10-fold. Thus, the resistance effect of the R346K mutation on the Mu variant is possible but not so significant and is due to the additional decrease of antibody neutralization based on the reduced binding of class 2 antibodies.

N501Y and K417N Mutations in the Spike Protein of SARS-CoV-2 Alter the Interactions with Both hACE2 and Human-Derived Antibody: A Free Energy of Perturbation Retrospective Study.

The N501Y and K417N mutations in the spike protein of SARS-CoV-2 and their combination gave rise to questions, but the data on their mechanism of action at the molecular level were limited. In this study, we present free energy perturbation (FEP) calculations, performed at the end of December 2020, for the interactions of the spike S1 receptor-binding domain (RBD) with both the ACE2 receptor and an antibody derived from COVID-19 patients. Our results showed that the S1 RBD-ACE2 interactions were significantly increased whereas those with the STE90-C11 antibody dramatically decreased. The K417N mutation in a combination with N501Y fully abolished the antibody effect. However, Lys417Asn seems to have a compensatory mechanism of action increasing the S1 RBD-ACE2 free energy of binding. This may explain the increased spread of the virus observed in the U.K. and South Africa and also gives rise to an important question regarding the possible human immune response and the success of the already available vaccines. Notably, when the experimental data became available confirming our calculations, it was demonstrated that protein-protein FEP can be a useful tool for providing urgent data to the scientific community.

1H-benzimidazole-2-yl hydrazones as tubulin-targeting agents: Synthesis, structural characterization, anthelmintic activity and antiproliferative activity against MCF-7 breast carcinoma cells and molecular docking studies.

In the present study, fifteen benzimidazolyl-2-hydrazones 7a-7o of fluoro-, hydroxy- and methoxy-substituted benzaldehydes and 1,3-benzodioxole-5-carbaldehyde were synthesized and their structure was identified by IR, NMR, and elemental analysis. The compounds 7j 2-(3-hydroxybenzylidene)-1-(5(6)-methyl-1H-benzimidazol-2-yl)hydrazone and 7i 2-(3-hydroxybenzylidene)-1-(1H-benzimidazol-2-yl)hydrazone have exerted the strongest anthelmintic activity (100% after 24 h incubation period at 37 °C) against isolated muscle larvae of Trichinella spiralis in an in vitro experiment. The in vitro cytotoxicity assay towards MCF-7 breast cancer cells and mouse embryo fibroblasts 3T3 showed that the studied benzimidazolyl-2-hydrazones exhibit low to moderate cytotoxic effects. The ability of the studied benzimidazolyl-2-hydrazones to modulate microtubule polymerization was confirmed and suggested that their anthelmintic action is mediated through inhibition of the tubulin polymerization likewise the other known benzimidazole anthelmitics. It was also shown that the four most promising benzimidazolyl-2-hydrazones do not affect significantly the AChE activity even at high tested concentration, thus indicating that they do not have the potential for neurotoxic effects. The binding mode of compounds 7j and 7n in the colchicine-binding site of tubulin were clarified by molecular docking simulations. Taken together, these results demonstrate that for the synthesized benzimidazole derivatives the anthelmintic activity against T. spiralis and the inhibition of tubulin polymerization are closely related.

Discovery of new AKT1 inhibitors by combination of in silico structure based virtual screening approaches and biological evaluations.

Communicated by Ramaswamy H. Sarma.

An Improved Free Energy Perturbation FEP+ Sampling Protocol for Flexible Ligand-Binding Domains.

Recent improvements to the free energy perturbation (FEP) calculations, especially FEP+ , established their utility for pharmaceutical lead optimization. Herein, we propose a modified version of the FEP/REST (i.e., replica exchange with solute tempering) sampling protocol, based on detail studies on several targets by probing a large number of perturbations with different sampling schemes. Improved FEP+ binding affinity predictions for regular flexible-loop motions and considerable structural changes can be obtained by extending the prior to REST (pre-REST) sampling time from 0.24 ns/λ to 5 ns/λ and 2 × 10 ns/λ, respectively. With this new protocol, much more precise ∆∆G values of the individual perturbations, including the sign of the transformations and decreased error were obtained. We extended the REST simulations from 5 ns to 8 ns to achieve reasonable free energy convergence. Implementing REST to the entire ligand as opposed to solely the perturbed region, and also some important flexible protein residues (pREST region) in the ligand binding domain (LBD) has considerably improved the FEP+ results in most of the studied cases. Preliminary molecular dynamics (MD) runs were useful for establishing the correct binding mode of the compounds and thus precise alignment for FEP+ . Our improved protocol may further increase the FEP+ accuracy.

Discovery of GlyT2 Inhibitors Using Structure-Based Pharmacophore Screening and Selectivity Studies by FEP+ Calculations.

In recent years, mammalian Glycine transporter 2 (GlyT2) has emerged as a promising target for the development of compounds against chronic pain states. In our current work, we discovered a new set of promising hits that inhibit the glycine transporter at nano- and micromolar activity and have excellent selectivity over GlyT1 (as shown by in vitro studies) using a newly designed virtual screening (VS) protocol that combines a structure-based pharmacophore and docking screens with a success rate of 75%. Furthermore, the free energy perturbation calculations and molecular dynamics (MD) studies revealed the GlyT2 amino acid residues critical for the binding and selectivity of both Glycine and our Hit1 compound. The FEP+ results well-matched with the available literature mutational data proving the quality of the generated GlyT2 structure. On the basis of these results, we propose that our hit compounds may lead to new chronic pain agents to address unmet and challenging clinical needs.

Elucidation of the orientation of selected drugs with 2-hydroxylpropyl-ß-cyclodextrin using 2D-NMR spectroscopy and molecular modeling.

This project aims to study the nature of interaction and orientation of selected drugs such as dexamethorphan HBr (DXM), diphenhydramine HCl (DPH), and lidocaine HCl (LDC) inclusion complexes with hydroxyl-propyl ß-cyclodextrin (HP-ß-CD) using 1HNMR spectroscopy, 2D-NMR ROESY and molecular-modeling techniques. Freeze-drying technique was used to formulate the inclusion complexes between DXM, DPH and LDC with HP-ß-CD (1:1 M ratio) in solid state. Inclusion complex formation was initially characterized by Fourier transform-infrared spectroscopy (FT-IR), differential scanning calorimetry (DSC), X-ray diffraction (XRD) and scanning electron microscopy (SEM) techniques. Further characterization of inclusion complexes to determine the interaction of DXM, DPH and LDC with HP-ß-CD was performed using the 1HNMR spectroscopy, 2D-NMR ROESY and molecular modeling techniques. Inclusion complexes of DXM, DPH and LDC with HP-ß-CD were successfully prepared using the freeze-drying technique. Preliminary studies with FT-IR, DSC, XRD and SEM indicated the formation of inclusion complexes of DXM, DPH and LDC with HP-ß-CD at 1:1 M ratio. 1HNMR study showed a change in proton chemical shift upon complexation. 2D-NMR ROESY (two-dimensional) spectroscopy gave an insight into the spatial arrangement between the host and guest atoms. 2D-ROESY experiments further predicted the direction of orientation of guest molecules, indicating the probability that amino moieties of DXM, DPH and LDC are inside the hydrophobic HP-ß-CD cavity. Cross-peaks of inclusion complexes demonstrated intermolecular nuclear Overhauser effects (NOE) between the amino protons in DXM, DPH and LDC and H-atoms of HP-ß-CD. Molecular modeling studies further confirmed the NMR data, providing a structural basis of the individual complex formations. Microsecond time-level molecular dynamics and metadynamics simulations indicate much stronger binding of DXM to HP-ß-CD and more dynamic behavior for DPH and LDC. In particular, LDC can exhibit multiple binding modes, and even spent some time (∼1-2%) out of the carrier, proving the dynamic nature of the complex. To conclude, 2D-NMR and molecular dynamic simulations elucidate the formation of inclusion complexes and intermolecular interactions of DXM, DPH and LDC with HP-ß-CD.

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.

PPARγ non-covalent antagonists exhibit mutable binding modes with a similar free energy of binding: a case study.

The structural and dynamical properties of PPARγ receptor in a complex with either partial or full agonists have been intensively studied but little is known about the receptor antagonistic conformation. A composition of microsecond accelerated molecular dynamics (aMD) simulation show that like partial agonists a non-covalent PPARγ full antagonist can bind in different modes of similar population size and free energies of binding. Four different and periodically exchanging ligand conformations are detected and described. The studied antagonist interacts with different receptor substructures and affects both the co-activator and the Cdk5 phosphorylation sites and, presumably, the natural complex with the DNA. However, no significant changes in the conformational states of the activation helix 12, and in particular an antagonist orientation, have been recorded. Finally, our results show also that the aMD approach can be successfully used in recovering the possible binding modes, considering fully the receptor flexibility, and is not dependent on the starting conformation.

PPARγ helix 12 exhibits an antagonist conformation.

Although the peroxisome proliferator-activated receptor γ (PPARγ) is one of the most studied nuclear receptors (NR), it is still unknown whether its activation helix (helix 12, H12) could exhibit an antagonist conformation as previously demonstrated for most of the NRs. The high H12 flexibility in the apo PPARγ form and the lack of appropriate antagonist ligands complicate the structural and dynamics description by most of the experimental techniques. Based on intensive (≈12 µs) accelerated molecular dynamics (aMD) simulations together with metadynamics and conventional MD runs, we reveal that H12 could exist in an antagonist conformation. This H12 state and the well-known agonist configuration have virtually identical free energy. Notably, significant deviations in the H12 conformations are detected in a homodimer. In chain A the activation helix is stabilized only in a full agonist conformation whereas in chain B, due to agonist to antagonist state exchanges, H12 is oriented toward helix 4. In summary, the results provide an explanation of the observed asymmetry in most of the PPARγ homodimer crystal structures. They also suggest selection guidance for protein moieties and structure candidates that would best serve as potential ligand binding sites to achieve a stable antagonist form of the receptor.

Structural and Dynamical Insight into PPARγ Antagonism: In Silico Study of the Ligand-Receptor Interactions of Non-Covalent Antagonists.

The structural and dynamical properties of the peroxisome proliferator-activated receptor γ (PPARγ) nuclear receptor have been broadly studied in its agonist state but little is known about the key features required for the receptor antagonistic activity. Here we report a series of molecular dynamics (MD) simulations in combination with free energy estimation of the recently discovered class of non-covalent PPARγ antagonists. Their binding modes and dynamical behavior are described in details. Two key interactions have been detected within the cavity between helices H3, H11 and the activation helix H12, as well as with H12. The strength of the ligand-amino acid residues interactions has been analyzed in relation to the specificity of the ligand dynamical and antagonistic features. According to our results, the PPARγ activation helix does not undergo dramatic conformational changes, as seen in other nuclear receptors, but rather perturbations that occur through a significant ligand-induced reshaping of the ligand-receptor and the receptor-coactivator binding pockets. The H12 residue Tyr473 and the charge clamp residue Glu471 play a central role for the receptor transformations. Our results also demonstrate that MD can be a helpful tool for the compound phenotype characterization (full agonists, partial agonists or antagonists) when insufficient experimental data are available.

Activation helix orientation of the estrogen receptor is mediated by receptor dimerization: evidence from molecular dynamics simulations.

In recent years, the nuclear receptors (NR) dynamics have been studied extensively by various approaches. However, the transition path of helix 12 (H12) to an agonist or an antagonist conformation and the exchange pathway between these states is not clear yet. A number of accelerated molecular dynamics (aMD) runs were performed on both an ERα monomer and a homodimer with a total length of 2.2 µs. We have been able to sample reasonably well the H12 conformational landscape to reproduce precisely both the agonist and the antagonist conformations, starting from an unfolded position, and to describe the transition path between them, even in the presence of an agonist ligand. These conformations were the most prevalent, suggesting that the extended H12 state is not likely to exist and that the natural ERα H12 position might exist in both the agonist and antagonist states. Remarkably, the H12 transition occurs and is regulated only in a dimer form and the proper agonist or antagonist H12 conformation can be achieved solely in one of the dimer subunits. These results clearly demonstrate that clusters of the two well-known H12 states exist by themselves in the protein free energy landscape, i.e. they are not constituted directly by the ligands, and dimerization favors the switch between them. Conversely, in a monomer, no transitions have been observed. Thus, the dimer formation helps the constitution of populations of discrete H12 conformational states and reshapes the conformational landscape. Further analyses have shown that these observations can be explained by specific interface and long range protein-protein interactions, resulting in conformational fluctuations in helices 5 and 11. Based on these results, a new ERα activation/deactivation mechanism and a sequence of binding events during receptor activity modulation have been suggested according to which ligands control the H12 conformation via alterations of the inter-dimer interactions. These findings agree with the HDX and fluorescence experiments and provide an explanation on a structural basis of these data, demonstrating that the dynamics of H12 are not altered greatly upon ligand binding and large fluctuations at the end of H11 are present.

Combination of genetic screening and molecular dynamics as a useful tool for identification of disease-related mutations: ZASP PDZ domain G54S mutation case.

Cypher/ZASP (LDB3 gene) is known to interact with a network of proteins. It binds to α-actinin and the calcium voltage channels (LTCC) via its PDZ domain. Here we report the identification of a highly conserved ZASP G54S mutation classified as a variant of unknown significance in a sample of an adult with hypertrophic cardiomyopathy (HCM). The initial bioinformatics calculations strongly evaluated G54S as damaging. Furthermore, we employed accelerated and classical molecular dynamics and free energy calculations to study the structural impact of this mutation on the ZASP apo form and to address the question of whether it can be linked to HCM. Seventeen independent MD runs and simulations of 2.5 µs total were performed and showed that G54S perturbs the α2 helix position via destabilization of the adjacent loop linked to the ß5 sheet. This also leads to the formation of a strong H-bond between peptide target residues Leu17 and Gln66, thus restricting both the α-actinin2 and LTCC C-terminal peptides to access their natural binding site and reducing in this way their binding capacity. On the basis of these observations and the adult's clinical data, we propose that ZASP(G54S) and presumably other ZASP PDZ domain mutations can cause HCM. To the best of our knowledge, this is the first reported ZASP PDZ domain mutation that might be linked to HCM. The integrated workflow used in this study can be applied for the identification and description of other mutations that might be related to particular diseases.

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.

Structural insight into the UNC-45-myosin complex.

The UNC-45 chaperone protein interacts with and affects the folding, stability, and the ATPase activity of myosins. It plays a critical role in the cardiomyopathy development and in the breast cancer tumor growth. Here we propose the first structural model of the UNC-45-myosin complex using various in silico methods. Initially, the human UNC-45B binding epitope was identified and the protein was docked to the cardiac myosin (MYH7) motor domain. The final UNC45B-MYH7 structure was obtained by performing of total 630 ns molecular dynamics simulations. The results indicate a complex formation, which is mainly stabilized by electrostatic interactions. Remarkably, the contact surface area is similar to that of the myosin-actin complex. A significant interspecies difference in the myosin binding epitope is observed. Our results reveal the structural basis of MYH7 exons 15-16 hypertrophic cardiomyopathy mutations and provide directions for drug targeting.

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.

A combination of 3D-QSAR, docking, local-binding energy (LBE) and GRID study of the species differences in the carcinogenicity of benzene derivatives chemicals.

A combination of 3D-QSAR, docking, local-binding energy (LBE) and GRID methods was applied as a tool to study and predict the mechanism of action of 100 carcinogenic benzene derivatives. Two 3D-QSAR models were obtained: (i) model of mouse carcinogenicity on the basis of 100 chemicals (model 1) and (ii) model of the differences in mouse and rat carcinogenicity on the basis of 73 compounds (model 2). 3D-QSAR regression maps indicated the important differences in species carcinogenicity, and the molecular positions associated with them. In order to evaluate the role of P450 metabolic process in carcinogenicity, the following approaches were used. The 3D models of CYP2E1 for mouse and rat were built up. A docking study was applied and the important ligand-protein residues interactions and oxidation positions of the molecules were identified. A new approach for quantitative assessment of metabolism pathways was developed, which enabled us to describe the species differences in CYP2E1 metabolism, and how it can be related to differences in the carcinogenic potential for a subset of compounds. The binding energies of the important substituents (local-binding energy-LBE) were calculated, in order to quantitatively demonstrate the contribution of the substituents in metabolic processes. Furthermore, a computational procedure was used for determining energetically favourable binding sites (GRID examination) of the enzymes. The GRID procedure allowed the identification of some important differences, related to species metabolism in CYP2E1. Comparing GRID, 3D-QSAR maps and LBE results, a similarity was identified, indicating a relationship between P450 metabolic processes and the differences in the carcinogenicity.