Mixed implicit/explicit solvation models in quantum mechanical calculations of binding enthalpy for protein–ligand complexes

An approach to quantum mechanical investigation of interactions in protein–ligand complexes has been developed that treats the solvation effect in a mixed scheme combining implicit and explicit solvent models. In this approach, the first solvation shell of the solvent around the solute is modeled with a limited number of hydrogen bonded explicit solvent molecules. The influence of the remaining bulk solvent is treated as a surrounding continuum in the conductor‐like screening model (COSMO). The enthalpy term of the binding free energy for the protein–ligand complexes was calculated using the semiempirical PM3 method implemented in the MOPAC package, applied to a trimmed model of the protein–ligand complex constructed with special rules. The dependence of the accuracy of binding enthalpy calculations on size of the trimmed model and number of optimized parameters was evaluated. Testing of the approach was performed for 12 complexes of different ligands with trypsin, thrombin, and ribonuclease with experimentally known binding enthalpies. The root‐mean‐square deviation (RMSD) of the calculated binding enthalpies from experimental data was found as ～1 kcal/mol over a large range. © 2006 Wiley Periodicals, Inc. Int J Quantum Chem, 2006     

Converging ligand‐binding free energies obtained with free‐energy perturbations at the quantum mechanical level

In this article, the convergence of quantum mechanical (QM) free‐energy simulations based on molecular dynamics simulations at the molecular mechanics (MM) level has been investigated. We have estimated relative free energies for the binding of nine cyclic carboxylate ligands to the octa‐acid deep‐cavity host, including the host, the ligand, and all water molecules within 4.5 Å of the ligand in the QM calculations (158–224 atoms). We use single‐step exponential averaging (ssEA) and the non‐Boltzmann Bennett acceptance ratio (NBB) methods to estimate QM/MM free energy with the semi‐empirical PM6‐DH2X method, both based on interaction energies. We show that ssEA with cumulant expansion gives a better convergence and uses half as many QM calculations as NBB, although the two methods give consistent results. With 720,000 QM calculations per transformation, QM/MM free‐energy estimates with a precision of 1 kJ/mol can be obtained for all eight relative energies with ssEA, showing that this approach can be used to calculate converged QM/MM binding free energies for realistic systems and large QM partitions. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Computation of the binding affinities of catechol‐O‐methyltransferase inhibitors: Multisubstate relative free energy calculations

Alchemical free energy simulations are amongst the most accurate techniques for the computation of the free energy changes associated with noncovalent protein–ligand interactions. A procedure is presented to estimate the relative binding free energies of several ligands to the same protein target where multiple, low‐energy configurational substates might coexist, as opposed to one unique structure. The contributions of all individual substates were estimated, explicitly, with the free energy perturbation method, and combined in a rigorous fashion to compute the overall relative binding free energies and dissociation constants. It is shown that, unless the most stable bound forms are known a priori, inaccurate results may be obtained if the contributions of multiple substates are ignored. The method was applied to study the complex formed between human catechol‐O‐methyltransferase and BIA 9‐1067, a newly developed tight‐binding inhibitor that is currently under clinical evaluation for the therapy of Parkinson's disease. Our results reveal an exceptionally high‐binding affinity (K_d in subpicomolar range) and provide insightful clues on the interactions and mechanism of inhibition. The inhibitor is, itself, a slowly reacting substrate of the target enzyme and is released from the complex in the form of O‐methylated product. By comparing the experimental catalytic rate (k_cat) and the estimated dissociation rate (k_off) constants of the enzyme‐inhibitor complex, one can conclude that the observed inhibition potency (K_i) is primarily dependent on the catalytic rate constant of the inhibitor's O‐methylation, rather than the rate constant of dissociation of the complex. © 2012 Wiley Periodicals, Inc.     

Coarse‐grained molecular dynamics simulations of protein–ligand binding

Coarse‐grained molecular dynamics (CGMD) simulations with the MARTINI force field were performed to reproduce the protein–ligand binding processes. We chose two protein–ligand systems, the levansucrase–sugar (glucose or sucrose), and LinB–1,2‐dichloroethane systems, as target systems that differ in terms of the size and shape of the ligand‐binding pocket and the physicochemical properties of the pocket and the ligand. Spatial distributions of the Coarse‐grained (CG) ligand molecules revealed potential ligand‐binding sites on the protein surfaces other than the real ligand‐binding sites. The ligands bound most strongly to the real ligand‐binding sites. The binding and unbinding rate constants obtained from the CGMD simulation of the levansucrase–sucrose system were approximately 10 times greater than the experimental values; this is mainly due to faster diffusion of the CG ligand in the CG water model. We could obtain dissociation constants close to the experimental values for both systems. Analysis of the ligand fluxes demonstrated that the CG ligand molecules entered the ligand‐binding pockets through specific pathways. The ligands tended to move through grooves on the protein surface. Thus, the CGMD simulations produced reasonable results for the two different systems overall and are useful for studying the protein–ligand binding processes. © 2014 Wiley Periodicals, Inc.     

Effect of interprotein polarization on protein–protein binding energy

Molecular dynamics simulation in explicit water for the binding of the benchmark barnase‐barstar complex was carried out to investigate the effect polarization of interprotein hydrogen bonds on its binding free energy. Our study is based on the AMBER force field but with polarized atomic charges derived from fragment quantum mechanical calculation for the protein complex. The quantum‐derived atomic charges include the effect of polarization of interprotein hydrogen bonds, which was absent in the standard force fields that were used in previous theoretical calculations of barnase‐barstar binding energy. This study shows that this polarization effect impacts both the static (electronic) and dynamic interprotein electrostatic interactions and significantly lowers the free energy of the barnase‐barstar complex. © 2012 Wiley Periodicals, Inc.     

Ab initio quantum mechanics-based free energy perturbation method for calculating relative solvation free energies

A free energy perturbation (FEP) method was developed that uses ab initio quantum mechanics (QM) for treating the solute molecules and molecular mechanics (MM) for treating the surroundings. Like our earlier results using AM1 semi empirical QMs, the ab initio QM/MM-based FEP method was shown to accurately calculate relative solvation free energies for a diverse set of small molecules that differ significantly in structure, aromaticity, hydrogen bonding potential, and electron density. Accuracy was similar to or better than conventional FEP methods. The QM/MM-based methods eliminate the need for time-consuming development of MM force field parameters, which are frequently required for drug-like molecules containing structural motifs not adequately described by MM. Future automation of the method and parallelization of the code for Linux 128/256/512 clusters is expected to enhance the speed and increase its use for drug design and lead optimization.     

Semiempirical calculations of binding enthalpy for protein–ligand complexes

A systematic semiempirical quantum mechanical study of the interactions between proteins and ligands has been performed to determine the ability of this approach for the accurate estimation of the enthalpic contribution to the binding free energy of the protein–ligand systems. This approach has been applied for eight test protein–ligand complexes with experimentally known binding enthalpies. The calculations were performed using the semiempirical PM3 approach incorporated in the MOPAC 97, ZAVA originally elaborated in Algodign, and MOPAC 2002 with MOZYME facility packages. Special attention was paid to take into account structural water molecules, which were located in the protein–ligand binding site. It was shown that the results of binding enthalpy calculations fit experimental data within ～2 kcal/mol in the presented approach. © 2003 Wiley Periodicals, Inc. Int J Quantum Chem, 2004     

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.     

Accurate prediction of relative binding affinities of a series of HIV-1 protease inhibitors using semi-empirical quantum mechanical charge

A major challenge in computer-aided drug design is the accurate estimation of ligand binding affinity. Here, a new approach that combines the adaptive steered molecular dynamics (ASMD) and partial atomic charges calculated by semi-empirical quantum mechanics (SQMPC), namely ASMD-SQMPC, is suggested to predict the ligand binding affinities, with 24 HIV-1 protease inhibitors as testing examples. In the ASMD-SQMPC, the relative binding free energy (ΔG) is reflected by the average maximum potential of mean force (<PMF>_max) between bound and unbound states. The correlation coefficient (R²) between the <PMF>_max and experimentally determined ΔG is 0.86, showing a significant improvement compared with the conventional ASMD (R² = 0.52). Therefore, this study provides an efficient approach to predict the relative ΔG and reveals the significance of precise partial atomic charges in the theoretical simulations.     

Two‐dimensional replica‐exchange method for predicting protein–ligand binding structures

We have developed a two‐dimensional replica‐exchange method for the prediction of protein–ligand binding structures. The first dimension is the umbrella sampling along the reaction coordinate, which is the distance between a protein binding pocket and a ligand. The second dimension is the solute tempering, in which the interaction between a ligand and a protein and water is weakened. The second dimension is introduced to make a ligand follow the umbrella potential more easily and enhance the binding events, which should improve the sampling efficiency. As test cases, we applied our method to two protein‐ligand complex systems (MDM2 and HSP 90‐alpha). Starting from the configuration in which the protein and the ligand are far away from each other in each system, our method predicted the ligand binding structures in excellent agreement with the experimental data from Protein Data Bank much faster with the improved sampling efficiency than the replica‐exchange umbrella sampling method that we have previously developed. © 2013 Wiley Periodicals, Inc.     

A method to calculate the one‐electron reduction potentials for nitroaromatic compounds based on gas‐phase quantum mechanics

Nitroaromatic compounds (NACs) are widespread environmental contaminants, and the one‐electron reduction potential (E) is an important parameter used in modeling their environmental fate. We have identified a method that is both accurate and efficient to predict E values for NACs, using gas‐phase quantum mechanics (QM) calculations combined with empirical correlations. First, the adiabatic electron affinity (EA) at 0 K is calculated using the B98/MG3S method, and the predictions are scaled by a factor of 0.802 to account for systematic errors in the density functional calculations. Second, the E values are predicted from a linear correlation between E and EA. Using this method, E values were predicted with a mean absolute deviation from measured values of 0.021 V for the 14 NACs used to obtain the correlation and 0.029 V for six additional NACs. This represents a substantial improvement in accuracy over predictions by other QM methods, which are affected by large errors in solvation or aqueous‐phase calculations for some compounds. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010.     

Three pillars for achieving quantum mechanical molecular dynamics simulations of huge systems: Divide‐and‐conquer,density‐functional tight‐binding,and massively parallel computation

The linear‐scaling divide‐and‐conquer (DC) quantum chemical methodology is applied to the density‐functional tight‐binding (DFTB) theory to develop a massively parallel program that achieves on‐the‐fly molecular reaction dynamics simulations of huge systems from scratch. The functions to perform large scale geometry optimization and molecular dynamics with DC‐DFTB potential energy surface are implemented to the program called DC‐DFTB‐K. A novel interpolation‐based algorithm is developed for parallelizing the determination of the Fermi level in the DC method. The performance of the DC‐DFTB‐K program is assessed using a laboratory computer and the K computer. Numerical tests show the high efficiency of the DC‐DFTB‐K program, a single‐point energy gradient calculation of a one‐million‐atom system is completed within 60 s using 7290 nodes of the K computer. © 2016 Wiley Periodicals, Inc.     

Site‐directed analysis on protein hydrophobicity

Hydrophobicity of a protein is considered to be one of the major intrinsic factors dictating the protein aggregation propensity. Understanding how protein hydrophobicity is determined is, therefore, of central importance in preventing protein aggregation diseases and in the biotechnological production of human therapeutics. Traditionally, protein hydrophobicity is estimated based on hydrophobicity scales determined for individual free amino acids, assuming that those scales are unaltered when amino acids are embedded in a protein. Here, we investigate how the hydrophobicity of constituent amino acid residues depends on the protein context. To this end, we analyze the hydration free energy—free energy change on hydration quantifying the hydrophobicity—of the wild‐type and 21 mutants of amyloid‐beta protein associated with Alzheimer's disease by performing molecular dynamics simulations and integral‐equation calculations. From detailed analysis of mutation effects on the protein hydrophobicity, we elucidate how the protein global factor such as the total charge as well as underlying protein conformations influence the hydrophobicity of amino acid residues. Our results provide a unique insight into the protein hydrophobicity for rationalizing and predicting the protein aggregation propensity on mutation, and open a new avenue to design aggregation‐resistant proteins as biotherapeutics. © 2014 Wiley Periodicals, Inc.     

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.     

Sequence‐based prediction of protein–peptide binding sites using support vector machine

Protein–peptide interactions are essential for all cellular processes including DNA repair, replication, gene‐expression, and metabolism. As most protein – peptide interactions are uncharacterized, it is cost effective to investigate them computationally as the first step. All existing approaches for predicting protein – peptide binding sites, however, are based on protein structures despite the fact that the structures for most proteins are not yet solved. This article proposes the first machine‐learning method called SPRINT to make Sequence‐based prediction of Protein – peptide Residue‐level Interactions. SPRINT yields a robust and consistent performance for 10‐fold cross validations and independent test. The most important feature is evolution‐generated sequence profiles. For the test set (1056 binding and non‐binding residues), it yields a Matthews’ Correlation Coefficient of 0.326 with a sensitivity of 64% and a specificity of 68%. This sequence‐based technique shows comparable or more accurate than structure‐based methods for peptide‐binding site prediction. SPRINT is available as an online server at: http://sparks-lab.org/ . © 2016 Wiley Periodicals, Inc.     

Update of the ATTRACT force field for the prediction of protein–protein binding affinity

Determining the protein–protein interactions is still a major challenge for molecular biology. Docking protocols has come of age in predicting the structure of macromolecular complexes. However, they still lack accuracy to estimate the binding affinities, the thermodynamic quantity that drives the formation of a complex. Here, an updated version of the protein–protein ATTRACT force field aiming at predicting experimental binding affinities is reported. It has been designed on a dataset of 218 protein–protein complexes. The correlation between the experimental and predicted affinities reaches 0.6, outperforming most of the available protocols. Focusing on a subset of rigid and flexible complexes, the performance raises to 0.76 and 0.69, respectively. © 2017 Wiley Periodicals, Inc.     

Electrostatic component of binding energy: Interpreting predictions from poisson–boltzmann equation and modeling protocols

Macromolecular interactions are essential for understanding numerous biological processes and are typically characterized by the binding free energy. Important component of the binding free energy is the electrostatics, which is frequently modeled via the solutions of the Poisson–Boltzmann Equations (PBE). However, numerous works have shown that the electrostatic component (ΔΔG_elec) of binding free energy is very sensitive to the parameters used and modeling protocol. This prompted some researchers to question the robustness of PBE in predicting ΔΔG_elec. We argue that the sensitivity of the absolute ΔΔG_elec calculated with PBE using different input parameters and definitions does not indicate PBE deficiency, rather this is what should be expected. We show how the apparent sensitivity should be interpreted in terms of the underlying changes in several numerous and physical parameters. We demonstrate that PBE approach is robust within each considered force field (CHARMM‐27, AMBER‐94, and OPLS‐AA) once the corresponding structures are energy minimized. This observation holds despite of using two different molecular surface definitions, pointing again that PBE delivers consistent results within particular force field. The fact that PBE delivered ΔΔG_elec values may differ if calculated with different modeling protocols is not a deficiency of PBE, but natural results of the differences of the force field parameters and potential functions for energy minimization. In addition, while the absolute ΔΔG_elec values calculated with different force field differ, their ordering remains practically the same allowing for consistent ranking despite of the force field used. © 2016 Wiley Periodicals, Inc.     

Carbohydrate‐binding protein identification by coupling structural similarity searching with binding affinity prediction

Carbohydrate‐binding proteins (CBPs) are potential biomarkers and drug targets. However, the interactions between carbohydrates and proteins are challenging to study experimentally and computationally because of their low binding affinity, high flexibility, and the lack of a linear sequence in carbohydrates as exists in RNA, DNA, and proteins. Here, we describe a structure‐based function‐prediction technique called SPOT‐Struc that identifies carbohydrate‐recognizing proteins and their binding amino acid residues by structural alignment program SPalign and binding affinity scoring according to a knowledge‐based statistical potential based on the distance‐scaled finite‐ideal gas reference state (DFIRE). The leave‐one‐out cross‐validation of the method on 113 carbohydrate‐binding domains and 3442 noncarbohydrate binding proteins yields a Matthews correlation coefficient of 0.56 for SPalign alone and 0.63 for SPOT‐Struc (SPalign + binding affinity scoring) for CBP prediction. SPOT‐Struc is a technique with high positive predictive value (79% correct predictions in all positive CBP predictions) with a reasonable sensitivity (52% positive predictions in all CBPs). The sensitivity of the method was changed slightly when applied to 31 APO (unbound) structures found in the protein databank (14/31 for APO versus 15/31 for HOLO). The result of SPOT‐Struc will not change significantly if highly homologous templates were used. SPOT‐Struc predicted 19 out of 2076 structural genome targets as CBPs. In particular, one uncharacterized protein in Bacillus subtilis (1oq1A) was matched to galectin‐9 from Mus musculus. Thus, SPOT‐Struc is useful for uncovering novel carbohydrate‐binding proteins. SPOT‐Struc is available at http://sparks‐lab.org . © 2014 Wiley Periodicals, Inc.