Toward quantitative estimates of binding affinities for protein–ligand systems involving large inhibitor compounds: A steered molecular dynamics simulation route

Understanding binding mechanisms between enzymes and potential inhibitors and quantifying protein – ligand affinities in terms of binding free energy is of primary importance in drug design studies. In this respect, several approaches based on molecular dynamics simulations, often combined with docking techniques, have been exploited to investigate the physicochemical properties of complexes of pharmaceutical interest. Even if the geometric properties of a modeled protein – ligand complex can be well predicted by computational methods, it is still challenging to rank with chemical accuracy a series of ligand analogues in a consistent way. In this article, we face this issue calculating relative binding free energies of a focal adhesion kinase, an important target for the development of anticancer drugs, with pyrrolopyrimidine‐based ligands having different inhibitory power. To this aim, we employ steered molecular dynamics simulations combined with nonequilibrium work theorems for free energy calculations. This technique proves very powerful when a series of ligand analogues is considered, allowing one to tackle estimation of protein – ligand relative binding free energies in a reasonable time. In our cases, the calculated binding affinities are comparable with those recovered from experiments by exploiting the Michaelis – Menten mechanism with a competitive inhibitor.     

Computational fragment-based design of Wee1 kinase inhibitors with tricyclic core scaffolds

Wee1 is cell cycle protein comprising a kinase domain and is a validated cancer target. We have designed molecules with variable tricyclic core scaffolds [6-6-5] system and extended them based on the chemical space available in the active site of Wee1 kinase using de novo drug design. The core scaffolds and linking fragments were extracted from pharmacophore-based virtual screening of ZINC and PubChem databases and Ludi library. These molecules bind the hinge region of kinase active site and form hydrogen bonds as confirmed from molecular docking, molecular dynamics simulations, and MM_PBSA calculations. When compared with reference inhibitors, AZD1775 and PHA-848125, the de novo designed molecules also show good docking scores and stability, retained non-covalent interactions, and high binding free energies contributed from active site residues.

     

Elucidating the Dual Mode of Action of Dipeptidyl Enoates in the Inhibition of Rhodesain Cysteine Proteases

A computational study of the two possible inhibition mechanisms of rhodesain cysteine protease by the dipeptidyl enoate Cbz-Phe-Leu-CH=CH−CO₂C₂H₅ has been carried out by means of molecular dynamics simulations with hybrid QM/MM potentials. The low free energy barriers confirm that the Cys25 residue can attack both C_β and C1 atoms of the inhibitor, confirming a dual mode of action in the inhibition of the rhodesain by enoates. According to the results, the inhibition process through the Cys25 attack on the C_β atom of the inhibitor is an exergonic and irreversible process, while the inhibition process when Cys25 attacks on the C1 atom of the inhibitor is and exergonic but reversible process. The interactions between the inhibitor and rhodesain suggest that P2 is the most important fragment to consider in the design of new efficient inhibitors of rhodesain. These results may be useful for the design of new inhibitors of rhodesain and other related cysteine proteases based on dipeptidyl enoates scaffolds.     

Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations

The Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA) and the Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) methods calculate binding free energies for macromolecules by combining molecular mechanics calculations and continuum solvation models. To systematically evaluate the performance of these methods, we report here an extensive study of 59 ligands interacting with six different proteins. First, we explored the effects of the length of the molecular dynamics (MD) simulation, ranging from 400 to 4800 ps, and the solute dielectric constant (1, 2, or 4) on the binding free energies predicted by MM/PBSA. The following three important conclusions could be observed: (1) MD simulation length has an obvious impact on the predictions, and longer MD simulation is not always necessary to achieve better predictions. (2) The predictions are quite sensitive to the solute dielectric constant, and this parameter should be carefully determined according to the characteristics of the protein/ligand binding interface. (3) Conformational entropy often show large fluctuations in MD trajectories, and a large number of snapshots are necessary to achieve stable predictions. Next, we evaluated the accuracy of the binding free energies calculated by three Generalized Born (GB) models. We found that the GB model developed by Onufriev and Case was the most successful model in ranking the binding affinities of the studied inhibitors. Finally, we evaluated the performance of MM/GBSA and MM/PBSA in predicting binding free energies. Our results showed that MM/PBSA performed better in calculating absolute, but not necessarily relative, binding free energies than MM/GBSA. Considering its computational efficiency, MM/GBSA can serve as a powerful tool in drug design, where correct ranking of inhibitors is often emphasized.     

Scaffold mining of kinase hinge binders in crystal structure database

Protein kinases are the second most prominent group of drug targets, after G-protein-coupled receptors. Despite their distinct inhibition mechanisms, the majority of kinase inhibitors engage the conserved hydrogen bond interactions with the backbone of hinge residues. We mined Pfizer internal crystal structure database (CSDb) comprising of several thousand of public as well as internal X-ray binary complexes to compile an inclusive list of hinge binding scaffolds. The minimum ring scaffolds with directly attached hetero-atoms and functional groups were extracted from the full compounds by applying a rule-based filtering procedure employing a comprehensive annotation of ATP-binding site of the human kinase complements. The results indicated large number of kinase inhibitors of diverse chemical structures are derived from a relatively small number of common scaffolds, which serve as the critical recognition elements for protein kinase interaction. Out of the nearly 4,000 kinase-inhibitor complexes in the CSDb we identified approximately 600 unique scaffolds. Hinge scaffolds are overwhelmingly flat with very little sp3 characteristics, and are less lipophilic than their corresponding parent compounds. Examples of the most common as well as the uncommon hinge scaffolds are presented. Although the most common scaffolds are found in complex with multiple kinase targets, a large number of them are uniquely bound to a specific kinase, suggesting certain scaffolds could be more promiscuous than the others. The compiled collection of hinge scaffolds along with their three-dimensional binding coordinates could serve as basis set for hinge hopping, a practice frequently employed to generate novel invention as well as to optimize existing leads in medicinal chemistry.     

Novel and Potential Small Molecule Scaffolds as DYRK1A Inhibitors by Integrated Molecular Docking-Based Virtual Screening and Dynamics Simulation Study

The dual-specificity tyrosine phosphorylation-regulated kinase 1A (DYRK1A) is a novel, promising and emerging biological target for therapeutic intervention in neurodegenerative diseases, especially in Alzheimer’s disease (AD). The molMall database, comprising rare, diverse and unique compounds, was explored for molecular docking-based virtual screening against the DYRK1A protein, in order to find out potential inhibitors. Ligands exhibiting hydrogen bond interactions with key amino acid residues such as Ile165, Lys188 (catalytic), Glu239 (gk+1), Leu241 (gk+3), Ser242, Asn244, and Asp307, of the target protein, were considered potential ligands. Hydrogen bond interactions with Leu241 (gk+3) were considered key determinants for the selection. High scoring structures were also docked by Glide XP docking in the active sites of twelve DYRK1A related protein kinases, viz. DYRK1B, DYRK2, CDK5/p25, CK1, CLK1, CLK3, GSK3β, MAPK2, MAPK10, PIM1, PKA, and PKCα, in order to find selective DYRK1A inhibitors. MM/GBSA binding free energies of selected ligand–protein complexes were also calculated in order to remove false positive hits. Physicochemical and pharmacokinetic properties of the selected six hit ligands were also computed and related with the proposed limits for orally active CNS drugs. The computational toxicity webserver ProTox-II was used to predict the toxicity profile of selected six hits (molmall IDs 9539, 11352, 15938, 19037, 21830 and 21878). The selected six docked ligand–protein systems were exposed to 100 ns molecular dynamics (MD) simulations to validate their mechanism of interactions and stability in the ATP pocket of human DYRK1A kinase. All six ligands were found to be stable in the ATP binding pocket of DYRK1A kinase.     

Grand canonical Monte Carlo simulation of ligand-protein binding

A new application of the grand canonical thermodynamics ensemble to compute ligand-protein binding is described. The described method is sufficiently rapid that it is practical to compute ligand-protein binding free energies for a large number of poses over the entire protein surface, thus identifying multiple putative ligand binding sites. In addition, the method computes binding free energies for a large number of poses. The method is demonstrated by the simulation of two protein-ligand systems, thermolysin and T4 lysozyme, for which there is extensive thermodynamic and crystallographic data for the binding of small, rigid ligands. These low-molecular-weight ligands correspond to the molecular fragments used in computational fragment-based drug design. The simulations correctly identified the experimental binding poses and rank ordered the affinities of ligands in each of these systems.     

An organometallic inhibitor for glycogen synthase kinase 3

Replacing natural products with kinetically inert metal complexes may lead to a new class of therapeutics in which a metal center plays the role of an innocent bystander, organizing the orientation of the organic ligands in the receptor space. As an example of this approach, a ruthenium complex is described that copies the binding mode of indolocarbazole protein kinase inhibitors and serves as a reversible, low-nanomolar inhibitor for glycogen synthase kinase 3 (GSK-3).     

Molecular dynamics simulation and free energy calculation studies of the binding mechanism of allosteric inhibitors with p38α MAP kinase

p38 MAP kinase is a promising target for anti-inflammatory treatment. The classical kinase inhibitors imatinib and sorafenib as well as BI-1 and BIRB-796 were reported to bind in the DFG-out form of human p38α, known as type II or allosteric kinase inhibitors. Although DFG-out conformation has attracted great interest in the design of type II kinase inhibitors, the structural requirements for binding and mechanism of stabilization of DFG-out conformation remain unclear. As allosteric inhibition is important to the selectivity of kinase inhibitor, herein the binding modes of imatinib, sorafenib, BI-1 and BIRB-796 to p38α were investigated by molecular dynamics simulation. Binding free energies were calculated by molecular mechanics/Poisson-Boltzmann surface area method. The predicted binding affinities can give a good explanation of the activity difference of the studied inhibitors. Furthermore, binding free energies decomposition analysis and further structural analysis indicate that the dominating effect of van der Waals interaction drives the binding process, and key residues, such as Lys53, Gly71, Leu75, Ile84, Thr106, Met109, Leu167, Asp168, and Phe169, play important roles by forming hydrogen bond, salt bridge, and hydrophobic interactions with the DFG-out conformation of p38α. Finally, we also conducted a detailed analysis of BI-1, imatinib, and sorafenib binding to p38α in comparison with BIRB-796 exploited for gaining potency as well as selectivity of p38 inhibitors. These results are expected to be useful for future rational design of novel type II p38 inhibitors.     

Alpha7 nicotinic acetylcholine receptor agonists: prediction of their binding affinity through a molecular mechanics Poisson-Boltzmann surface area approach

A group of agonists for the alpha7 neuronal nicotinic acetylcholine receptors (nAChRs) was investigated, and their free energies of binding DeltaG(bind) were calculated by applying the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) approach. This method, based on molecular dynamics simulations of fully solvated protein-ligand complexes, allowed us to estimate the contribution of both polar and nonpolar terms as well as the entropy to the overall free energy of binding. The calculated results were in a good agreement with the experimentally determined DeltaG(bind) values, thereby pointing to the MM-PBSA protocol as a valuable computational tool for the rational design of specific agents targeting the neuronal alpha7 nAChR subtypes.     

Postprocessing of docked protein-ligand complexes using implicit solvation models

Molecular docking plays an important role in drug discovery as a tool for the structure-based design of small organic ligands for macromolecules. Possible applications of docking are identification of the bioactive conformation of a protein-ligand complex and the ranking of different ligands with respect to their strength of binding to a particular target. We have investigated the effect of implicit water on the postprocessing of binding poses generated by molecular docking using MM-PB/GB-SA (molecular mechanics Poisson-Boltzmann and generalized Born surface area) methodology. The investigation was divided into three parts: geometry optimization, pose selection, and estimation of the relative binding energies of docked protein-ligand complexes. Appropriate geometry optimization afforded more accurate binding poses for 20% of the complexes investigated. The time required for this step was greatly reduced by minimizing the energy of the binding site using GB solvation models rather than minimizing the entire complex using the PB model. By optimizing the geometries of docking poses using the GB(HCT+SA) model then calculating their free energies of binding using the PB implicit solvent model, binding poses similar to those observed in crystal structures were obtained. Rescoring of these poses according to their calculated binding energies resulted in improved correlations with experimental binding data. These correlations could be further improved by applying the postprocessing to several of the most highly ranked poses rather than focusing exclusively on the top-scored pose. The postprocessing protocol was successfully applied to the analysis of a set of Factor Xa inhibitors and a set of glycopeptide ligands for the class II major histocompatibility complex (MHC) A(q) protein. These results indicate that the protocol for the postprocessing of docked protein-ligand complexes developed in this paper may be generally useful for structure-based design in drug discovery.     

基于虚拟筛选和分子动力学模拟发现新型ULK1抑制剂

Unc-51样自噬激活激酶1（unc-51-like autophagy activating kinase 1,ULK1）作为自噬启动的重要调控因子,是肿瘤治疗的关键靶点之一。首先,以已知ULK1抑制剂为基础构建药效团模型,通过药效团模型筛选、分子对接以及分子力学广义波恩表面积（Molecular Mechanics/Generalized Born Surface Area,MM/GBSA）结合自由能计算等方法,对含有52万多个类药性小分子的数据库进行虚拟筛选,得到具有较高理论亲和力的化合物。随后,50ns的分子动力学模拟验证了蛋白质-配体复合物结合的稳定性,最后10ns的平均结合自由能的计算研究进一步验证了配体的结合能力。结果表明,6个化合物（F5258-0159、F3407-0428、F0529-1100、F0696-3531、F3222-5280、F6525-5596）具有骨架新颖、分子对接分数和结合自由能数值优异及与ULK1的结合状态稳定等特点,可以作为新型潜在的ULK1抑制剂用于肿瘤治疗的研究,也为新型ULK1抑制剂的设计和研发提供新的研究思路。     

Direct calculation of the binding free energies of FKBP ligands

Direct calculations of the absolute free energies of binding for eight ligands to FKBP protein were performed using the Fujitsu BioServer massively parallel computer. Using the latest version of the general assisted model building with energy refinement (AMBER) force field for ligand model parameters and the Bennett acceptance ratio for computing free-energy differences, we obtained an excellent linear fit between the calculated and experimental binding free energies. The rms error from a linear fit is 0.4 kcal/mol for eight ligand complexes. In comparison with a previous study of the binding energies of these same eight ligand complexes, these results suggest that the use of improved model parameters can lead to more predictive binding estimates, and that these estimates can be obtained with significantly less computer time than previously thought. These findings make such direct methods more attractive for use in rational drug design.     

A computational analysis of the binding model of MDM2 with inhibitors

Abstract  

It is a new and promising strategy for anticancer drug design to block the MDM2-p53 interaction using a non-peptide small-molecule inhibitor. We carry out molecular dynamics simulations to study the binding of a set of six non-peptide small-molecule inhibitors with the MDM2. The relative binding free energies calculated using molecular mechanics Poisson–Boltzmann surface area method produce a good correlation with experimentally determined results. The study shows that the van der Waals energies are the largest component of the binding free energy for each complex, which indicates that the affinities of these inhibitors for MDM2 are dominated by shape complementarity. The A-ligands and the B-ligands are the same except for the conformation of 2,2-dimethylbutane group. The quantum mechanics and the binding free energies calculation also show the B-ligands are the more possible conformation of ligands. Detailed binding free energies between inhibitors and individual protein residues are calculated to provide insights into the inhibitor-protein binding model through interpretation of the structural and energetic results from the simulations. The study shows that G1, G2 and G3 group mimic the Phe19, Trp23 and Leu26 residues in p53 and their interactions with MDM2, but the binding model of G4 group differs from the original design strategy to mimic Leu22 residue in p53.     

Computational study of ligand binding in lipid transfer proteins: Structures,interfaces, and free energies of protein‐lipid complexes

Plant nonspecific lipid transfer proteins (nsLTPs) bind a wide variety of lipids, which allows them to perform disparate functions. Recent reports on their multifunctionality in plant growth processes have posed new questions on the versatile binding abilities of these proteins. The lack of binding specificity has been customarily explained in qualitative terms on the basis of a supposed structural flexibility and nonspecificity of hydrophobic protein‐ligand interactions. We present here a computational study of protein‐ligand complexes formed between five nsLTPs and seven lipids bound in two different ways in every receptor protein. After optimizing geometries in molecular dynamics calculations, we computed Poisson‐Boltzmann electrostatic potentials, solvation energies, properties of the protein‐ligand interfaces, and estimates of binding free energies of the resulting complexes. Our results provide the first quantitative information on the ligand abilities of nsLTPs, shed new light into protein‐lipid interactions, and reveal new features which supplement commonly held assumptions on their lack of binding specificity. © 2012 Wiley Periodicals, Inc.     

Standard free energy of releasing a localized water molecule from the binding pockets of proteins: double-decoupling method

Localized water molecules in the binding pockets of proteins play an important role in noncovalent association of proteins and small drug compounds. At times, the dominant contribution to the binding free energy comes from the release of localized water molecules in the binding pockets of biomolecules. Therefore, to quantify the energetic importance of these water molecules for drug design purposes, we have used the double-decoupling approach to calculate the standard free energy of tying up a water molecule in the binding pockets of two protein complexes. The double-decoupling approach is based on the underlying principle of statistical thermodynamics. We have calculated the standard free energies of tying up the water molecule in the binding pockets of these complexes to be favorable. These water molecules stabilize the protein-drug complexes by interacting with the ligands and binding pockets. Our results offer ideas that could be used in optimizing protein-drug interactions, by designing ligands that are capable of targeting localized water molecules in protein binding sites. The resulting free energy of ligand binding could benefit from the potential free energy gain accompanying the release of these water molecules. Furthermore, we have examined the theoretical background of the double-decoupling method and its connection to the molecular dynamics thermodynamic integration techniques.     

含平面胺配体的反式二价钯配合物与DNA碱基的作用

含大的平面胺配体的二价钯金属配合物在当前的抗肿瘤药物设计中代表着一类具有重要发展前途的先导结构. 由于大的平面胺配体具有较大的空间位阻, 目前主要的问题是这类化合物能否和DNA碱基结合形成单功能和双功能加合物. 我们采用密度泛函理论和等电聚焦连续极化(IEF-PCM)溶剂化模型研究了trans-PdCl2L2(L:2-羟基吡啶)的钯配合物与DNA碱基的作用. 该化合物与DNA形成单功能和双功能加合物反应的活化自由能均低于铂类抗肿瘤药. 所有反应在水溶液中均为放热反应. 结果表明, 这一大的平面胺配体不会阻碍该化合物与DNA碱基形成双功能加合物, 而且该化合物与DNA的单功能和双功能结合的速率会大于铂类化合物.     

New parameterization approaches of the LIE method to improve free energy calculations of PlmII‐Inhibitors complexes

The standard parameterization of the Linear Interaction Energy (LIE) method has been applied with quite good results to reproduce the experimental absolute binding free energies for several protein–ligand systems. However, we found that this parameterization failed to reproduce the experimental binding free energy of Plasmepsin II (PlmII) in complexes with inhibitors belonging to four dissimilar scaffolds. To overcome this fact, we developed three approaches of LIE, which combine systematic approaches to predict the inhibitor‐specific values of α, β, and γ parameters, to gauge their ability to calculate the absolute binding free energies for these PlmII‐Inhibitor complexes. Specifically: (i) we modified the linear relationship between the weighted nonpolar desolvation ratio (WNDR) and the α parameter, by introducing two models of the β parameter determined by the free energy perturbation (FEP) method in the absence of the constant term γ, and (ii) we developed a new parameterization model to investigate the linear correlation between WNDR and the correction term γ. Using these parameterizations, we were able to reproduce the experimental binding free energy from these systems with mean absolute errors lower than 1.5 kcal/mol. © 2010 Wiley Periodicals, Inc. J Comput Chem 2010     

Identification of Novel Covalent XPO1 Inhibitors Based on a Hybrid Virtual Screening Strategy

Nuclear export protein 1 (XPO1), a member of the nuclear export protein-p (Karyopherin-P) superfamily, regulates the transport of “cargo” proteins. To facilitate this important process, which is essential for cellular homeostasis, XPO1 must first recognize and bind the cargo proteins. To inhibit this process, small molecule inhibitors have been designed that inhibit XPO1 activity through covalent binding. However, the scaffolds for these inhibitors are very limited. While virtual screening may be used to expand the diversity of the XPO1 inhibitor skeleton, enormous computational resources would be required to accomplish this using traditional screening methods. In the present study, we report the development of a hybrid virtual screening workflow and its application in XPO1 covalent inhibitor screening. After screening, several promising XPO1 covalent molecules were obtained. Of these, compound 8 performed well in both tumor cell proliferation assays and a nuclear export inhibition assay. In addition, molecular dynamics simulations were performed to provide information on the mode of interaction of compound 8 with XPO1. This research has identified a promising new scaffold for XPO1 inhibitors, and it demonstrates an effective and resource-saving workflow for identifying new covalent inhibitors.