A water-swap reaction coordinate for the calculation of absolute protein-ligand binding free energies

The accurate prediction of absolute protein-ligand binding free energies is one of the grand challenge problems of computational science. Binding free energy measures the strength of binding between a ligand and a protein, and an algorithm that would allow its accurate prediction would be a powerful tool for rational drug design. Here we present the development of a new method that allows for the absolute binding free energy of a protein-ligand complex to be calculated from first principles, using a single simulation. Our method involves the use of a novel reaction coordinate that swaps a ligand bound to a protein with an equivalent volume of bulk water. This water-swap reaction coordinate is built using an identity constraint, which identifies a cluster of water molecules from bulk water that occupies the same volume as the ligand in the protein active site. A dual topology algorithm is then used to swap the ligand from the active site with the identified water cluster from bulk water. The free energy is then calculated using replica exchange thermodynamic integration. This returns the free energy change of simultaneously transferring the ligand to bulk water, as an equivalent volume of bulk water is transferred back to the protein active site. This, directly, is the absolute binding free energy. It should be noted that while this reaction coordinate models the binding process directly, an accurate force field and sufficient sampling are still required to allow for the binding free energy to be predicted correctly. In this paper we present the details and development of this method, and demonstrate how the potential of mean force along the water-swap coordinate can be improved by calibrating the soft-core Coulomb and Lennard-Jones parameters used for the dual topology calculation. The optimal parameters were applied to calculations of protein-ligand binding free energies of a neuraminidase inhibitor (oseltamivir), with these results compared to experiment. These results demonstrate that the water-swap coordinate provides a viable and potentially powerful new route for the prediction of protein-ligand binding free energies.     

Exploring hierarchical refinement techniques for induced fit docking with protein and ligand flexibility

We present a series of molecular‐mechanics‐based protein refinement methods, including two novel ones, applied as part of an induced fit docking procedure. The methods used include minimization; protein and ligand sidechain prediction; a hierarchical ligand placement procedure similar to a‐priori protein loop predictions; and a minimized Monte Carlo approach using normal mode analysis as a move step. The results clearly indicate the importance of a proper opening of the active site backbone, which might not be accomplished when the ligand degrees of freedom are prioritized. The most accurate method consisted of the minimized Monte Carlo procedure designed to open the active site followed by a hierarchical optimization of the sidechain packing around a mobile flexible ligand. The methods have been used on a series of 88 protein‐ligand complexes including both cross‐docking and apo‐docking members resulting in complex conformations determined to within 2.0 Å heavy‐atom RMSD in 75% of cases where the protein backbone rearrangement upon binding is less than 1.0 Å α‐carbon RMSD. We also demonstrate that physics‐based all‐atom potentials can be more accurate than docking‐style potentials when complexes are sufficiently refined. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Active site pressurization: a new tool for structure-guided drug design and other studies of protein flexibility

We present a new molecular dynamics methodology to assist in structure-based drug design and other studies that seek to predict protein deformability. Termed Active Site Pressurization (ASP), the new methodology simply injects a resin into the ligand binding-site of a protein during the course of a molecular dynamics simulation such that novel, energetically reasonable protein conformations are generated in an unbiased way that may be better representations of the ligand binding conformation than are currently available. Here we apply two different versions of the ASP methodology to three proteins, cytochrome P450cam, PcrA helicase, and glycogen synthase kinase 3beta (GSK3beta), and show that the method is capable of inducing significant conformational changes when compared to the X-ray crystal structures. Application of the ASP methodology therefore provides a view of binding site flexibility that is a rich source of data for inclusion in a variety of further investigations, including high-throughput virtual screening, lead hopping, revealing alternative modes of deformation, and revealing hidden exit and entrance tunnels.     

Modeling loop backbone flexibility in receptor‐ligand docking simulations

The relevance of receptor conformational change during ligand binding is well documented for many pharmaceutically relevant receptors, but is still not fully accounted for in in silico docking methods. While there has been significant progress in treatment of receptor side chain flexibility sampling of backbone flexibility remains challenging because the conformational space expands dramatically and the scoring function must balance protein–protein and protein–ligand contributions. Here, we investigate an efficient multistage backbone reconstruction algorithm for large loop regions in the receptor and demonstrate that treatment of backbone receptor flexibility significantly improves binding mode prediction starting from apo structures and in cross docking simulations. For three different kinase receptors in which large flexible loops reconstruct upon ligand binding, we demonstrate that treatment of backbone flexibility results in accurate models of the complexes in simulations starting from the apo structure. At the example of the DFG‐motif in the p38 kinase, we also show how loop reconstruction can be used to model allosteric binding. Our approach thus paves the way to treat the complex process of receptor reconstruction upon ligand binding in docking simulations and may help to design new ligands with high specificity by exploitation of allosteric mechanisms. © 2012 Wiley Periodicals, Inc.     

Fully automated molecular mechanics based induced fit protein-ligand docking method

We describe a method for docking a ligand into a protein receptor while allowing flexibility of the protein binding site. The method employs a multistep procedure that begins with the generation of protein and ligand conformations. An initial placement of the ligand is then performed by computing binding site hotspots. This initial placement is followed by a protein side-chain refinement stage that models protein flexibility. The final step of the process is an energy minimization of the ligand pose in the presence of the rigid receptor. Thus the algorithm models flexibility of the protein at two stages, before and after ligand placement. We validated this method by performing docking and cross docking studies of eight protein systems for which crystal structures were available for at least two bound ligands. The resulting rmsd values of the 21 docked protein-ligand complexes showed values of 2 A or less for all but one of the systems examined. The method has two critical benefits for high throughput virtual screening studies. First, no user intervention is required in the docking once the initial binding site selection has been made in the protein. Second, the initial protein conformation generation needs to be performed only once for a given binding region. Also, the method may be customized in various ways depending on the particular scenario in which dockings are being performed. Each of the individual steps of the method is fully independent making it straightforward to explore different variants of the high level workflow to further improve accuracy and performance.     

A structure-guided approach for protein pocket modeling and affinity prediction

Binding affinity prediction is frequently addressed using computational models constructed solely with molecular structure and activity data. We present a hybrid structure-guided strategy that combines molecular similarity, docking, and multiple-instance learning such that information from protein structures can be used to inform models of structure–activity relationships. The Surflex-QMOD approach has been shown to produce accurate predictions of binding affinity by constructing an interpretable physical model of a binding site with no experimental binding site structural information. We introduce a method to integrate protein structure information into the model induction process in order to construct more robust physical models. The structure-guided models accurately predict binding affinities over a broad range of compounds while producing more accurate representations of the protein pockets and ligand binding modes. Structure-guidance for the QMOD method yielded significant performance improvements, both for affinity and pose prediction, especially in cases where predictions were made on ligands very different from those used for model induction.     

Assessing protein loop flexibility by hierarchical Monte Carlo sampling

Loop flexibility is often crucial to protein biological function in solution. We report a new Monte Carlo method for generating conformational ensembles for protein loops and cyclic peptides. The approach incorporates the triaxial loop closure method which addresses the inverse kinematic problem for generating backbone move sets that do not break the loop. Sidechains are sampled together with the backbone in a hierarchical way, making it possible to make large moves that cross energy barriers. As an initial application, we apply the method to the flexible loop in triosephosphate isomerase that caps the active site, and demonstrate that the resulting loop ensembles agree well with key observations from previous structural studies. We also demonstrate, with 3 other test cases, the ability to distinguish relatively flexible and rigid loops within the same protein.     

喹啉加氧酶中氧气扩散的通道分析

通过生物信息学分析、量化计算优化、CAVER和MDpocket预测、随机加速分子动力学及伞状抽样动力学模拟等方法,对喹啉加氧酶(HOD)中的氧气扩散途径进行了计算预测.结果表明,氧气在HOD中的反应位点包埋在蛋白内部,而HOD中有数条可能的通道供氧气进出,其中长度最短的通道具有最高的优先度,不仅在随机加速动力学模拟中具有最高的氧气逸出概率,而且伞状抽样方法计算得到的自由能也最低.此通道的内端位于底物Re面的氧气结合位点,较好地解释了HOD的相关实验数据.     

A pose prediction approach based on ligand 3D shape similarity

Molecular docking predicts the best pose of a ligand in the target protein binding site by sampling and scoring numerous conformations and orientations of the ligand. Failures in pose prediction are often due to either insufficient sampling or scoring function errors. To improve the accuracy of pose prediction by tackling the sampling problem, we have developed a method of pose prediction using shape similarity. It first places a ligand conformation of the highest 3D shape similarity with known crystal structure ligands into protein binding site and then refines the pose by repacking the side-chains and performing energy minimization with a Monte Carlo algorithm. We have assessed our method utilizing CSARdock 2012 and 2014 benchmark exercise datasets consisting of co-crystal structures from eight proteins. Our results revealed that ligand 3D shape similarity could substitute conformational and orientational sampling if at least one suitable co-crystal structure is available. Our method identified poses within 2 Å RMSD as the top-ranking pose for 85.7 % of the test cases. The median RMSD for our pose prediction method was found to be 0.81 Å and was better than methods performing extensive conformational and orientational sampling within target protein binding sites. Furthermore, our method was better than similar methods utilizing ligand 3D shape similarity for pose prediction.     

Monitoring Conformational Changes in the NDM‐1 Metallo‐β‐lactamase by 19F NMR Spectroscopy

The New Delhi metallo‐β‐lactamase (NDM‐1) is involved in the emerging antibiotic resistance problem. Development of metallo‐β‐lactamases (MBLs) inhibitors has proven challenging, due to their conformational flexibility. Here we report site‐selective labeling of NDM‐1 with 1,1,1‐trifluoro‐3‐bromo acetone (BFA), and its use to study binding events and conformational changes upon ligand–metal binding using ¹⁹F NMR spectroscopy. The results demonstrate different modes of binding of known NDM‐1 inhibitors, including L ‐ and D ‐captopril by monitoring the changing chemical environment of the active‐site loop of NDM‐1. The method described will be applicable to other MBLs and more generally to monitoring ligand‐induced conformational changes.     

Identification of New Lead Molecules Against UBE2NL Enzyme for Cancer Therapy

Cancer is characterized by abnormal growth of cells. Targeting ubiquitin proteins in the discovery of new anticancer therapeutics is an attractive strategy. The present study uses the structure-based drug discovery methods to identify new lead structures, which are selective to the putative ubiquitin-conjugating enzyme E2N-like (UBE2NL). The 3D structure of the UBE2NL was evaluated using homology modeling techniques. The model was validated using standard in silico methods. The hydrophobic pocket of UBE2NL that aids in binding with its natural receptor ubiquitin-conjugating enzyme E2 variant (UBE2V) was identified through protein-protein docking study. The binding site region of the UBE2NL was identified using active site prediction tools. The binding site of UBE2NL which is responsible for cancer cell progression is considered for docking study. Virtual screening study with the small molecular structural database was carried out against the active site of UBE2NL. The ligand molecules that have shown affinity towards UBE2NL were considered for ADME prediction studies. The ligand molecules that obey the Lipinski’s rule of five and Jorgensen’s rule of three pharmacokinetic properties like human oral absorption etc. are prioritized. The resultant ligand molecules can be considered for the development of potent UBE2NL enzyme inhibitors for cancer therapy.     

HIV蛋白酶抑制剂——利迪链菌素的分子对接研究

用分子对接的方法, 对利迪链菌素的抗HIV蛋白酶活性进行了研究. 为了更准确地反映利迪链菌素分子与酶蛋白结合的情况, 充分考虑受体活性部位的柔性, 采用了FlexX(初步对接)和Flexidock(精确对接)分两步将配体与受体进行对接. 在初步对接中, 设计了不同的受体活性部位来考察是否有结合水分子参与抑制剂与酶的结合. 对一种作用方式已知的非肽类HIV蛋白酶抑制剂Aha006进行的对接研究显示, 分子模拟的结果与实际情况吻合得较好, 证明了本文所采用的方法的可靠性. 利迪链菌素与蛋白酶活性部位的对接结果显示, 配体分子与受体之间的结合没有结合水分子的参与, 两者通过5对氢键作用结合成为稳定的复合物. 利迪链菌素占据结合腔, 覆盖了蛋白酶的活性三联体Asp25-Thr26-Gly27, 从而起到抑制其生物活性的作用.     

Improved mapping of protein binding sites

Computational mapping methods place molecular probes – small molecules or functional groups – on a protein surface in order to identify the most favorable binding positions by calculating an interaction potential. Mapping is an important step in a number of flexible docking and drug design algorithms. We have developed improved algorithms for mapping protein surfaces using small organic molecules as molecular probes. The calculations reproduce the binding of eight organic solvents to lysozyme as observed by NMR, as well as the binding of four solvents to thermolysin, in good agreement with x-ray data. Application to protein tyrosine phosphatase 1B shows that the information provided by the mapping can be very useful for drug design. We also studied why the organic solvents bind in the active site of proteins, in spite of the availability of alternative pockets that can very tightly accommodate some of the probes. A possible explanation is that the binding in the relatively large active site retains a number of rotational states, and hence leads to smaller entropy loss than the binding elsewhere else. Indeed, the mapping reveals that the clusters of the ligand molecules in the protein's active site contain different rotational-translational conformers, which represent different local minima of the free energy surface. In order to study the transitions between different conformers, reaction path and molecular dynamics calculations were performed. Results show that most of the rotational states are separated by low free energy barriers at the experimental temperature, and hence the entropy of binding in the active site is expected to be high.     

Theoretical study of the truncated hemoglobin HbN: exploring the molecular basis of the NO detoxification mechanism

Mycobacterium tuberculosis is the causative agent of human tuberculosis. The nitric oxide reaction with oxy-truncated hemoglobin N (trHbN) has been proposed to be responsible for the resistance mechanism by which this microorganism can evade the toxic effects of NO. In this work, we explore the molecular basis of the NO detoxification mechanism using a combination of classical and hybrid quantum-classical (QM-MM) simulation techniques. We have investigated the structural flexibility of the protein, the ligand affinity properties, and the nitric oxide reaction with coordinated O2. The analysis of the classical MD trajectory allowed us to identify Phe62 as the gate of the main channel for ligand diffusion to the active site. Moreover, the opening of the channel stems from the interplay between collective backbone motions and local rearrangements in the side chains of the residues that form the bottleneck of the tunnel. Even though the protein environment is not found to make a significant contribution to the heme moiety catalyzed reaction, the binding site influences the physiological function of the enzyme at three different levels. First, by isolating the intermediates formed in the reaction, it prevents nondesired reactions from proceeding. Second, it modulates the ligand (O2, NO) affinity of the protein, which can be ascribed to both distal and proximal effects. Finally, the stabilization of the Tyr33-Gln58 pair upon O2 binding might alter the essential dynamics of the protein, leading in turn to a mechanism for ligand-induced regulation.     

Prediction of binding constants of protein ligands: A fast method for the prioritization of hits obtained from de novo design or 3D database search programs

A dataset of 82 protein–ligand complexes of known 3D structure and binding constant K_i was analysed to elucidate the important factors that determine the strength of protein–ligand interactions. The following parameters were investigated: the number and geometry of hydrogen bonds and ionic interactions between the protein and the ligand, the size of the lipophilic contact surface, the flexibility of the ligand, the electrostatic potential in the binding site, water molecules in the binding site, cavities along the protein–ligand interface and specific interactions between aromatic rings. Based on these parameters, a new empirical scoring function is presented that estimates the free energy of binding for a protein–ligand complex of known 3D structure. The function distinguishes between buried and solvent accessible hydrogen bonds. It tolerates deviations in the hydrogen bond geometry of up to 0.25 Å in the length and up to 30 °Cs in the hydrogen bond angle without penalizing the score. The new energy function reproduces the binding constants (ranging from 3.7 × 10^-2 M to 1 × 10^-14 M, corresponding to binding energies between -8 and -80 kJ/mol) of the dataset with a standard deviation of 7.3 kJ/mol corresponding to 1.3 orders of magnitude in binding affinity. The function can be evaluated very fast and is therefore also suitable for the application in a 3D database search or de novo ligand design program such as LUDI. The physical significance of the individual contributions is discussed.     

Influence of the membrane lipophilic environment on the structure and on the substrate access/egress routes of the human aromatase enzyme. A computational study

Human aromatase (HA), an enzyme located on the membrane of the endoplasmatic reticulum, is of crucial biological importance in the biosynthesis of estrogens. High levels of estrogens are related with important pathologies, conferring to HA a key role as a pharmacological target. In this study we provide, for the first time, an atomistic model of HA embedded on a membrane model to understand the influence of the membrane lipophilic environment on the structural and dynamical properties of HA and on the access/egress pathways of the substrate (androstenedione, ASD) and of the oxygen molecule (involved in the enzymatic process) into/from the HA active site. To this end we used several computational techniques such as force field-based molecular dynamics (MD) simulations, Random Expulsion MD, Steered MD, and Implicit Ligand Sampling. Our results show that the membrane anchoring does not markedly affect the structural properties and the flexibility of the protein, but they clearly point out that the membrane has a marked effect on the access/egress routes of the reactants, stabilizing the formation of different channels for both ASD and O(2) with respect to those observed in pure water solution. Due to the importance of HA in medicine and since access/egress channels may influence its substrate selectivity, a detailed understanding of the role of the membrane in shaping these channels may be of valuable help in drug design.