Binding free energy calculation with QM/MM hybrid methods for Abl-Kinase inhibitor

Electronic and structural study of T315I mutated form in DFG-out conformation of BCR-ABL inhibitors

In this work, the four main drugs for the treatment of chronic myeloid leukemia were analyzed, being imatinib, dasatinib, nilotinib and ponatinib followed by four derivative molecules of nilotinib and ponatinib. For these derivative molecules, the fluorine atoms were replaced by hydrogen and chlorine atoms in order to shade light to the structural effects on this set of inhibitors. Electronic studies were performed at density functional theory level with the B3LYP functional and 6-311+G(d,p) basis set. The frontier molecular orbitals, gap HOMO-LUMO, and NBO were analyzed and compared to docking studies for mutant T315I tyrosine kinase protein structure code 3IK3, in the DFG-out conformation. Structural similarities were pointed out, such as the presence of groups common to all inhibitors and modifications raised up on new generations of imatinib-based inhibitors. One of them is the trifluoromethyl group present in nilotinib and later included in ponatinib, in addition to the 1-methylpiperazin-1-ium group that is present in imatinib and ponatinib. The frontier molecular orbitals of imatinib and ponatinib are contributing to the same amino acid residues, and the ineffectiveness of imatinib against the T315I mutation was discussed.Communicated by Ramaswamy H. Sarma.

Elucidation of Ligand/Protein Interactions between BCR-ABL Tyrosine Kinase and Some Commercial Anticancer Drugs Via DFT Methods

In this study, the interactions of 7 commercially available BCR-ABL tyrosine kinase enzyme inhibitors with amino acids in the active site of the relevant enzyme were investigated quantum mechanically. Here a per-residue study was carried out. Interaction energies were calculated by using the coordinates of the critical residues in the binding site of the enzyme and the drug molecules docked in this region. DFT methods were used during the QM processes. All interaction energies were calculated via M06-2X functional and 6-31G (d,p) basis set in vacuum. Based on the results obtained, it was tried to be determined which of the important residues in the binding cavity of the enzyme could better interact with the examined ligands. It is thought that this study may contribute to the development of tyrosine kinase enzyme inhibitors.

Protein–ligand free energies of binding from full-protein DFT calculations: convergence and choice of exchange–correlation functional

Quantum mechanical binding free energies based on thousands of full-protein DFT calculations are tractable, reproducible and converge well.

Simulating protein-ligand binding with neural network potentials

Drug molecules adopt a range of conformations both in solution and in their protein-bound state. The strain and reduced flexibility of bound drugs can partially counter the intermolecular interactions that drive protein-ligand binding. To make accurate computational predictions of drug binding affinities, computational chemists have attempted to develop efficient empirical models of these interactions, although these methods are not always reliable. Machine learning has allowed the development of highly-accurate neural-network potentials (NNPs), which are capable of predicting the stability of molecular conformations with accuracy comparable to state-of-the-art quantum chemical calculations but at a billionth of the computational cost. Here, we demonstrate that these methods can be used to represent the intramolecular forces of protein-bound drugs within molecular dynamics simulations. These simulations are shown to be capable of predicting the protein-ligand binding pose and conformational component of the absolute Gibbs energy of binding for a set of drug molecules. Notably, the conformational energy for anti-cancer drug erlotinib binding to its target was found to be considerably overestimated by a molecular mechanical model, while the NNP predicts a more moderate value. Although the ANI-1ccX NNP was not trained to describe ionic molecules, reasonable binding poses are predicted for charged ligands, but this method is not suitable for modeling charged ligands in solution.

Designing new generation of potent inhibitors against membrane-type matrix metalloproteinase-2: a computational effort against multiple myeloma

Matrix metalloproteinases (MMPs) play important roles in cancer progression and, despite their inhibitors have failed in the clinical trials, they have always been considered as suitable targets for the treatment of tumor. We have recently shown that membrane type (MT) 2-MMPs, is selectively expressed in multiple myeloma (MM) cells and mediates the metastatic characteristics of these cells. In this study, we designed efficient inhibitors against MT2-MMP using state-of-art molecular modeling methods. First, the 3D structure of MT2-MMP was predicted. Then, the proposed potent inhibitors against two regions of the catalytic domain of MT2-MMP (active site and MT-LOOP) were identified through molecular docking, QM-MM and molecular dynamics simulations from a set of compounds in Analyticon library, IBS library, Maybridge screening fragment library and drugbank library. Moreover, ADME estimation showed that pharmacokinetic properties of inhibitors are in the acceptable range for humans. Finally, our data suggested that compounds 'structures.722' (dobutamine) and 'M2' are suitable candidates to inhibit MT2-MMP for further examination in the laboratory.Communicated by Ramaswamy H. Sarma.

Molecular docking and ADMET analysis of hydroxamic acids as HDAC2 inhibitors

Histone deacetylase (HDAC2) belongs to the hydrolase family and a promising target for cancers. We reported 96 hydroxamic compounds optimized using hydrogen-donors, hydrophobic and electron withdrawing groups followed by molecular docking studies. The optimized compounds show good LibDock score and H-bond interaction in the active site of HDAC2. We selected 20 compounds as the best HDAC2 inhibitors based on the LibDock score, binding energy and hydrogen bonding. ADMET predictions on these compounds show good absorption, BBB penetration and no liver toxicity. We subsequently report four compounds selected as best HDAC2 inhibitors based on the LibDock, binding energy, H-bonding and ADMET properties.

Structural insights into the anti-cancer activity of quercetin on G-tetrad, mixed G-tetrad, and G-quadruplex DNA using quantum chemical and molecular dynamics simulations

Human telomerase referred as 'terminal transferase' is a nucleoprotein enzyme which inhibits the disintegration of telomere length and act as a drug target for the anticancer therapy. The tandem repeating structure of telomere sequence forms the guanine-rich quadruplex structures that stabilize stacked tetrads. In our present work, we have investigated the interaction of quercetin with DNA tetrads using DFT. Geometrical analysis revealed that the influence of quercetin drug induces the structural changes into the DNA tetrads. Among DNA tetrads, the quercetin stacked with GCGC tetrad has the highest interaction energy of -88.08 kcal/mol. The binding mode and the structural stability are verified by the absorption spectroscopy method. The longer wavelength was found at 380 nm and it exhibits bathochromic shift. The findings help us to understand the binding nature of quercetin drug with DNA tetrads and it also inhibits the telomerase activity. Further, the quercetin drug interacted with G-quadruplex DNA by using molecular dynamics (MD) simulation studies for 100 ns simulation at different temperatures and different pH levels (T = 298 K, 320 K and pH = 7.4, 5.4). The structural stability of the quercetin with G-quadruplex structure is confirmed by RMSD. For the acidic condition (pH = 5.4), the binding affinity is higher toward G-quadruplex DNA, this result resembles that the quercetin drug is well interacted with G-quadruplex DNA at acidic condition (pH = 7.4) than the neutral condition. The obtained results show that quercetin drug stabilizes the G-quadruplex DNA, which regulates telomerase enzyme and it potentially acts as a novel anti-cancer agent.Communicated by Ramaswamy H. Sarma.

Combined UMC- DFT prediction of electron-hole coupling in unit cells of pentacene crystals

Pentacene is an organic semiconductor that draws special attention from the scientific community due to the high mobility of its charge carriers. As electron-hole interactions are important aspects in the regard of such property, a computationally inexpensive method to predict the coupling between these quasi-particles is highly desired. In this work, we propose a hybrid methodology of combining Uncoupled Monte Carlo Simulations (UMC) and Density functional Theory (DFT) methodologies to obtain a good compromise between computational feasibility and accuracy. As a first step in considering a Pentacene crystal, we describe its unit cell: the Pentacene Dimer. Because many conformations can be encountered for the dimer and considering the complexity of the system, we make use of UMC in order to find the most probable structures and relative orientations for the Pentacene-Pentacene complex. Following, we carry out electronic structure calculations in the scope of DFT with the goal of describing the electron-hole coupling on the most probable configurations obtained by UMC. The comparison of our results with previously reported data on the literature suggests that the methodology is well suited for describing transfer integrals of organic semiconductors. The observed accuracy together with the smaller computational cost required by our approach allows us to conclude that such methodology might be an important tool towards the description of systems with higher complexity.

Insight on Mutation-Induced Resistance from Molecular Dynamics Simulations of the Native and Mutated CSF-1R and KIT

The receptors tyrosine kinases (RTKs) for the colony stimulating factor-1, CSF-1R, and for the stem cell factor, SCFR or KIT, are important mediators of signal transduction. The abnormal function of these receptors, promoted by gain-of-function mutations, leads to their constitutive activation, associated with cancer or other proliferative diseases. A secondary effect of the mutations is the alteration of receptors' sensitivity to tyrosine kinase inhibitors, compromising effectiveness of these molecules in clinical treatment. In particular, the mutation V560G in KIT increases its sensitivity to Imatinib, while the D816V in KIT, and D802V in CSF-1R, triggers resistance to the drug. We analyzed the Imatinib binding affinity to the native and mutated KIT (mutations V560G, S628N and D816V) and CSF-1R (mutation D802V) by using molecular dynamics simulations and energy calculations of Imatinib•target complexes. Further, we evaluated the sensitivity of the studied KIT receptors to Imatinib by measuring the inhibition of KIT phosphorylation. Our study showed that (i) the binding free energy of Imatinib to the targets is highly correlated with their experimentally measured sensitivity; (ii) the electrostatic interactions are a decisive factor affecting the binding energy; (iii) the most deleterious impact to the Imatinib sensitivity is promoted by D802V (CSF-1R) and D816V (KIT) mutations; (iv) the role of the juxtamembrane region, JMR, in the imatinib binding is accessory. These findings contribute to a better description of the mutation-induced effects alternating the targets sensitivity to Imatinib.

Binding energies of tyrosine kinase inhibitors: Error assessment of computational methods for imatinib and nilotinib binding

The binding energies of imatinib and nilotinib to tyrosine kinase have been determined by quantum mechanical (QM) computations, and compared with literature binding energy studies using molecular mechanics (MM). The potential errors in the computational methods include these critical factors. Errors in X-ray structures such as structural distortions and steric clashes give unrealistically high van der Waals energies, and erroneous binding energies.MM optimization gives a very different configuration to the QM optimization for nilotinib, whereas the imatinib ion gives similar configurations. Solvation energies are a major component of the overall binding energy. The QM based solvent model (PCM/SMD) gives different values from those used in the implicit PBSA solvent MM models. A major error in inhibitor—kinase binding lies in the non-polar solvation terms. Solvent transfer free energies and the required empirical solvent accessible surface area factors for nilotinib and imatinib ion to give the transfer free energies have been reverse calculated. These values differ from those used in the MM PBSA studies.An intertwined desolvation—conformational binding selectivity process is a balance of thermodynamic desolvation and intramolecular conformational kinetic control.The configurational entropies (TΔS) are minor error sources.

Computational analysis of the binding specificity of Gleevec to Abl, c-Kit, Lck, and c-Src tyrosine kinases

Gleevec, a well-known cancer therapeutic agent, is an effective inhibitor of several tyrosine kinases, including Abl and c-Kit, but displays less potency to inhibit closely homologous tyrosine kinases, such as Lck and c-Src. Because many structural features of the binding site are highly conserved in these homologous kinases, the molecular determinants responsible for the binding specificity of Gleevec remain poorly understood. To address this issue, free energy perturbation molecular dynamics (FEP/MD) simulations with explicit solvent was used to compute the binding affinity of Gleevec to Abl, c-Kit, Lck, and c-Src. The results of the FEP/MD calculations are in good agreement with experiments, enabling a detailed and quantitative dissection of the absolute binding free energy in terms of various thermodynamic contributions affecting the binding specificity of Gleevec to the kinases. Dominant binding free energy contributions arises from the van der Waals dispersive interaction, compensating about two-thirds of the unfavorable free energy penalty associated with the loss of translational, rotational, and conformational freedom of the ligand upon binding. In contrast, the contributions from electrostatic and repulsive interactions nearly cancel out due to solvent effects. Furthermore, the calculations show the importance of the conformation of the kinase activation loop. Among the kinases examined, Abl provides the most favorable binding environment for Gleevec via optimal protein-ligand interactions and a small free energy cost for loss of the translational, rotational, and conformational freedom upon ligand binding. The FEP/MD calculations additionally reveal that Lck and c-Src provide similar nonbinding interactions with the bound-Gleevec, but the former pays less entropic penalty for the ligand losing its translational, rotational, and conformational motions to bind, examining the empirically observed differential binding affinities of Gleevec between the two Src-family kinases.

Influence of PNA containing 8-aza-7-deazaadenine on structure stability and binding affinity of PNA·DNA duplex: insights from thermodynamics, counter ion, hydration and molecular dynamics analysis

This paper describes the synthesis of a novel 8-aza-7-deazapurin-2,6-diamine (DPP)-containing peptide nucleic acid (PNA) monomer and Boc protecting group-based oligomerization of PNA, replacing adenine (A) with DPP monomers in the PNA strand. The PNA oligomers were synthesized against the biologically relevant SV40 promoter region (2494-AATTTTTTTTATTTA-2508) of pEGFP-N3 plasmid. The DPP-PNA·DNA duplex showed enhanced stability as compared to normal duplex (A-PNA·DNA). The electronic distribution of DPP monomer suggested that DPP had better electron donor properties over 2,6-diamino purine. UV melting and thermodynamic analysis revealed that the PNA oligomer containing a diaminopyrazolo(3,4-d)pyrimidine moiety (DPP) stabilized the PNA·DNA hybrids compared to A-PNA·DNA. DPP-PNA·DNA duplex showed higher water activity (Δnw = 38.5) in comparison to A-PNA·DNA duplex (Δnw = 14.5). The 50 ns molecular dynamics simulations of PNA·DNA duplex containing DPP or unmodified nucleobase-A showed average H-bond distances in the DPP-dT base pair of 2.90 Å (OH-N bond) and 2.91 Å (NH-N bond), which were comparably shorter than in the A-dT base pair, in which the average distances were 3.18 Å (OH-N bond) and 2.97 Å (NH-N bond), and there was one additional H-bond in the DPP-dT base pair of around 2.98 Å (O2H-N2 bond), supporting the higher stability of DPP-PNA·DNA. The analysis of molecular dynamics simulation data showed that the system binding free energy increased at a rate of approximately -4.5 kcal mol(-1) per DPP base of the PNA·DNA duplex. In summary, increased thermal stability, stronger hydrogen bonding and more stable conformation in the DPP-PNA·DNA duplex make it a better candidate as antisense/antigene therapeutic agents.

Explaining why Gleevec is a specific and potent inhibitor of Abl kinase

Tyrosine kinases present attractive drug targets for specific types of cancers. Gleevec, a well-known therapeutic agent against chronic myelogenous leukemia, is an effective inhibitor of Abl tyrosine kinase. However, Gleevec fails to inhibit closely homologous tyrosine kinases, such as c-Src. Because many structural features of the binding site are conserved, the molecular determinants responsible for binding specificity are not immediately apparent. Some have attributed the difference in binding specificity of Gleevec to subtle variations in ligand-protein interactions (binding affinity control), whereas others have proposed that it is the conformation of the DFG motif, in which ligand binding is only accessible to Abl and not to c-Src (conformational selection control). To address this issue, the absolute binding free energy was computed using all-atom molecular dynamics simulations with explicit solvent. The results of the free energy simulations are in good agreement with experiments, thereby enabling a meaningful decomposition of the binding free energy to elucidate the factors controlling Gleevec's binding specificity. The latter is shown to be controlled by a conformational selection mechanism and also by differences in key van der Waals interactions responsible for the stabilization of Gleevec in the binding pocket of Abl.

Calculation of Host-Guest Binding Affinities Using a Quantum-Mechanical Energy Model

The prediction of protein-ligand binding affinities is of central interest in computer-aided drug discovery, but it is still difficult to achieve a high degree of accuracy. Recent studies suggesting that available force fields may be a key source of error motivate the present study, which reports the first mining minima (M2) binding affinity calculations based on a quantum mechanical energy model, rather than an empirical force field. We apply a semi-empirical quantum-mechanical energy function, PM6-DH+, coupled with the COSMO solvation model, to 29 host-guest systems with a wide range of measured binding affinities. After correction for a systematic error, which appears to derive from the treatment of polar solvation, the computed absolute binding affinities agree well with experimental measurements, with a mean error 1.6 kcal/mol and a correlation coefficient of 0.91. These calculations also delineate the contributions of various energy components, including solute energy, configurational entropy, and solvation free energy, to the binding free energies of these host-guest complexes. Comparison with our previous calculations, which used empirical force fields, point to significant differences in both the energetic and entropic components of the binding free energy. The present study demonstrates successful combination of a quantum mechanical Hamiltonian with the M2 affinity method.

Stability and free energy calculation of LNA modified quadruplex: a molecular dynamics study

Telomeric ends of chromosomes, which comprise noncoding repeat sequences of guanine-rich DNA, which are the fundamental in protecting the cell from recombination and degradation. Telomeric DNA sequences can form four stranded quadruplex structures, which are involved in the structure of telomere ends. The formation and stabilization of telomeric quadruplexes has been shown to inhibit the activity of telomerase, thus establishing telomeric DNA quadrulex as an attractive target for cancer therapeutic intervention. Molecular dynamic simulation offers the prospects of detailed description of the dynamical structure with ion and water at molecular level. In this work we have taken a oligomeric part of human telomeric DNA, d(TAGGGT) to form different monomeric quadruplex structures d(TAGGGT)₄. Here we report the relative stabilities of these structures under K⁺ ion conditions and binding interaction between the strands, as determined by molecular dynamic simulations followed by energy calculation. We have taken locked nucleic acid (LNA) in this study. The free energy molecular mechanics Poission Boltzman surface area calculations are performed for the determination of most stable complex structure between all modified structures. We calculated binding free energy for the combination of different strands as the ligand and receptor for all structures. The energetic study shows that, a mixed hybrid type quadruplex conformation in which two parallel strands are bind with other two antiparallel strands, are more stable than other conformations. The possible mechanism for the inhibition of the cancerous growth has been discussed. Such studies may be helpful for the rational drug designing.