QM/MM−PBSA Method To Estimate Free Energies for Reactions in Proteins

Fuentes L, Pimenta A, Nagem RAP, Chávez-Olórtegui C, Schneider FS, Molina F, Sanchez EF, Suárez D, Ferreira RS. Assessing the Interactions between Snake Venom Metalloproteinases and Hydroxamate Inhibitors Using Kinetic and ITC Assays, Molecular Dynamics Simulations and MM/PBSA-Based Scoring Functions

Bothrops species are the main cause of snake bites in rural communities of tropical developing countries of Central and South America. Envenomation by Bothrops snakes is characterized by prominent local inflammation, hemorrhage and necrosis as well as systemic hemostatic disturbances. These pathological effects are mainly caused by the major toxins of the viperidae venoms, the snake venom metalloproteinases (SVMPs). Despite the antivenom therapy efficiency to block the main toxic effects on bite victims, this treatment shows limited efficacy to prevent tissue necrosis. Thus, drug-like inhibitors of these toxins have the potential to aid serum therapy of accidents inflicted by viper snakes. Broad-spectrum metalloprotease inhibitors bearing a hydroxamate zinc-binding group are potential candidates to improve snake bites therapy and could also be used to study toxin-ligand interactions. Therefore, in this work, we used both docking calculations and molecular dynamics simulations to assess the interactions between six hydroxamate inhibitors and two P-I SVMPs selected as models: Atroxlysin-I (hemorrhagic) from Bothrops atrox, and Leucurolysin-a (nonhemorrhagic) from Bothrops leucurus. We also employed a large variety of end-point free energy methods in combination with entropic terms to produce scoring functions of the relative affinities of the inhibitors for the toxins. Then we identified the scoring functions that best correlated with experimental data obtained from kinetic activity assays. In addition, to the characterization of these six molecules as inhibitors of the toxins, this study sheds light on the main enzyme-inhibitor interactions, explaining the broad-spectrum behavior of the inhibitors, and identifies the energetic and entropic terms that improve the performance of the scoring functions.

Computational methods for calculation of protein-ligand binding affinities in structure-based drug design

Drug design is an expensive and time-consuming process. Any method that allows reducing the time the costs of the drug development project can have great practical value for the pharmaceutical industry. In structure-based drug design, affinity prediction methods are of great importance. The majority of methods used to predict binding free energy in protein-ligand complexes use molecular mechanics methods. However, many limitations of these methods in describing interactions exist. An attempt to go beyond these limits is the application of quantum-mechanical description for all or only part of the analyzed system. However, the extensive use of quantum mechanical (QM) approaches in drug discovery is still a demanding challenge. This chapter briefly reviews selected methods used to calculate protein-ligand binding affinity applied in virtual screening (VS), rescoring of docked poses, and lead optimization stage, including QM methods based on molecular simulations.

Accurate Binding Free Energy Method from End-State MD Simulations

Herein, we introduce a new strategy to estimate binding free energies using end-state molecular dynamics simulation trajectories. The method is adopted from linear interaction energy (LIE) and ANI-2x neural network potentials (machine learning) for the atomic simulation environment (ASE). It predicts the single-point interaction energies between ligand–protein and ligand–solvent pairs at the accuracy of the wb97x/6-31G* level for the conformational space that is sampled by molecular dynamics (MD) simulations. Our results on 54 protein–ligand complexes show that the method can be accurate and have a correlation of R = 0.87–0.88 to the experimental binding free energies, outperforming current end-state methods with reduced computational cost. The method also allows us to compare BFEs of ligands with different scaffolds. The code is available free of charge (documentation and test files) at https://github.com/otayfuroglu/deepQM.

Comparative assessment of QM-based and MM-based models for prediction of protein-ligand binding affinity trends

Methods which accurately predict protein-ligand binding strengths are critical for drug discovery. In the last two decades, advances in chemical modelling have enabled steadily accelerating progress in the discovery and optimization of structure-based drug design. Most computational methods currently used in this context are based on molecular mechanics force fields that often have deficiencies in describing the quantum mechanical (QM) aspects of molecular binding. In this study, we show the competitiveness of our QM-based Molecules-in-Molecules (MIM) fragmentation method for characterizing binding energy trends for seven different datasets of protein-ligand complexes. By using molecular fragmentation, the MIM method allows for accelerated QM calculations. We demonstrate that for classes of structurally similar ligands bound to a common receptor, MIM provides excellent correlation to experiment, surpassing the more popular Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) and Molecular Mechanics Generalized Born Surface Area (MM/GBSA) methods. The MIM method offers a relatively simple, well-defined protocol by which binding trends can be ascertained at the QM level and is suggested as a promising option for lead optimization in structure-based drug design.

QM/MM Energy Decomposition Using the Interacting Quantum Atoms Approach

The interacting quantum atoms (IQA) method decomposes the quantum mechanical (QM) energy of a molecular system in terms of one- and two-center (atomic) contributions within the context of the quantum theory of atoms in molecules. Here, we demonstrate that IQA, enhanced with molecular mechanics (MM) and Poisson-Boltzmann surface-area (PBSA) solvation methods, is naturally extended to the realm of hybrid QM/MM methodologies, yielding intra- and inter-residue energy terms that characterize all kinds of covalent and noncovalent bonding interactions. To test the robustness of this approach, both metal-water interactions and QM/MM boundary artifacts are characterized in terms of the IQA descriptors derived from QM regions of varying size in Zn(II)- and Mg(II)-water clusters. In addition, we analyze a homologous series of inhibitors in complex with a matrix metalloproteinase (MMP-12) by carrying out QM/MM-PBSA calculations on their crystallographic structures followed by IQA energy decomposition. Overall, these applications not only show the advantages of the IQA QM/MM approach but also address some of the challenges lying ahead for expanding the QM/MM methodology.

Assessing the Performance of Screening MM/PBSA in Protein-Ligand Interactions

Accurate calculation of the binding free energies between a protein and a ligand is the primary objective of structure-based drug design, but it still remains a challenging problem. In this work, we apply the screening molecular mechanics/Poisson Boltzmann surface area (MM/PBSA) method to calculate the binding affinity of protein-ligand interactions. Our results show that the performance of the screening MM/PBSA is better than that of the standard MM/PBSA, especially in a charged-ligand system. In addition, we also investigate the effect of the solute dielectric constant on the results, and find that the optimal solute dielectric constants are different between the neutral-ligand system and the charged-ligand system. Moreover, we also evaluate the effect of the atomic-charge methods on the performance of the screening MM/PBSA. The present study demonstrates that the screening MM/PBSA should be a reliable method for calculating binding energy of biosystems.

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.

An accurate free energy estimator: based on MM/PBSA combined with interaction entropy for protein-ligand binding affinity

The molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) method is constantly used to calculate the binding free energy of protein-ligand complexes, and has been shown to effectively balance computational cost against accuracy. The relative binding affinities obtained by the MM/PBSA approach are acceptable, while it usually overestimates the absolute binding free energy. This paper proposes four free energy estimators based on the MM/PBSA for enthalpy change combined with interaction entropy (IE) for entropy change using different weights for individual energy terms. The ΔG_{PBSA_IE} method is determined to be an optimal estimator based on its performance in terms of the correlation between experimental and theoretical values and error estimations. This approach is optimized using high-quality experimental values from a training set containing 84 protein-ligand systems, and the coefficients for the sum of electrostatic energy and polar solvation free energy, van der Waals (vdW) energy, non-polar solvation energy and entropy change are obtained by multivariate linear fitting to the corresponding experimental values. A comparison between the traditional MM/PBSA method and this method shows that the correlation coefficient is improved from 0.46 to 0.72 and the slope of the regression line increases from 0.10 to 1.00. More importantly, the mean absolute error (MAE) is significantly reduced from 22.52 to 1.59 kcal mol^-1. Furthermore, the numerical stability of this method is validated on a test set with a similar correlation coefficient, slope and MAE to those of the training set. Based on the above advantages, the ΔG_{PBSA_IE} method can be a powerful tool for a reliable and accurate estimation of binding free energy and plays a significant role in a detailed energetic investigation of protein-ligand interaction.

Understanding PIM-1 kinase inhibitor interactions with free energy simulation

The proviral integration site of the Moloney leukemia virus (PIM) family includes three homologous members. PIM-1 kinase is an important target in effective therapeutic interventions of lymphomas, prostate cancer and leukemia. In the current work, we performed free energy calculations to calculate the binding affinities of several inhibitors targeting this protein. The alchemical method with integration and perturbation-based estimators and the end-point methods were compared. The computational results indicated that the alchemical method can accurately predict the binding affinities, while the end-point methods give relatively unreliable predictions. Decomposing the free energy difference into enthalpic and entropic components with MBAR reweighting enabled us to investigate the detailed thermodynamic parameters with which the entropy-enthalpy compensation in this protein-ligand binding case is identified. We then studied the conformational ensemble, and the important protein-ligand interactions were identified. The current work sheds light on the understanding of the PIM-1-kinase-inhibitor interactions at the atomic level and will be useful in the further development of potential drugs.

Affinity Calculations of Cyclodextrin Host-Guest Complexes: Assessment of Strengths and Weaknesses of End-Point Free Energy Methods

The end-point methods like MM/PBSA or MM/GBSA estimate the free energy of a biomolecule by combining its molecular mechanics energy with solvation free energy and entropy terms. On the one hand, their performance largely depends on the particular system of interest, and despite numerous attempts to improve their reliability that have resulted in many variants, there is still no clear alternative to improve their accuracy. On the other hand, the relatively small cyclodextrin host-guest complexes, for which high-quality binding calorimetric data are usually available, are becoming reference models for testing the accuracy of free energy methods. In this work, we further assess the performance of various MM/PBSA-like approaches as applied to cyclodextrin complexes. To this end, we select a set of complexes between β-cyclodextrin and 57 small organic molecules that has been previously studied with the binding energy distribution analysis method in combination with an implicit solvent model ( Wickstrom, L.; He, P.; Gallicchio, E.; Levy, R. M. J. Chem. Theory Comput. 2013 , 9 , 3136 - 3150 ). For each complex, a conventional 1.0 μs molecular dynamics simulation in explicit solvent is performed. Then we employ semiempirical quantum chemical calculations, several variants of the MM-PB(GB)SA methods, entropy estimations, etc., to assess the reliability of the end-point affinity calculations. The best end-point protocol in this study, which combines DFTB3 energies with entropy corrections, yields estimations of the binding free energies that still have substantial errors (RMSE = 2.2 kcal/mol), but it exhibits a good prediction capacity in terms of ligand ranking ( R² = 0.66) that is close to or even better than that of rigorous free energy methodologies. Our results can be helpful to discriminate between the intrinsic limitations of the end-point methods and other sources of error, such as the underlying energy and continuum solvation methods.

From bench to bedside, via desktop. Recent advances in the application of cutting-edge in silico tools in the research of drugs targeting bromodomain modules

The discipline of drug discovery has greatly benefited by computational tools and in silico algorithms aiming at rationalization of many related processes, from the stage of early hit identification to the preclinical phases of drug candidate validation. The various methodologies referred to as molecular modeling tools span a broad spectrum of applications, from straightforward approaches such as virtual screening of compound libraries to more advanced techniques involving the precise estimation of free energy upon binding of the candidate drug to its macromolecular target. To this end, we report an overview of specific studies where implementation of such sophisticated modeling algorithms has shown to be indispensable for addressing challenging systems and biological questions otherwise difficult to answer. We focus our attention on the emerging field of bromodomain inhibitors. Bromodomains are small modules involved in epigenetic signaling and currently comprise high-priority targets for developing both drug candidates and chemical probes for basic biomedical research. We attempt a critical presentation of selected cases utilizing cutting-edge in silico methodologies, with our main emphasis being on absolute or relative free energy simulations, on implementation of quantum-mechanics level calculations and on characterization of solvent thermodynamics. We discuss the advantages and strengths as well as the drawbacks and weaknesses of computational tools utilized in those works and we attempt to comment on specific issues related to their integration into the regular medicinal chemistry practice. Our conclusion is that while such methods indeed represent highly promising resources for further advancing the discipline, their application is not always trivial.

Solvation energies of the proton in methanol revisited and temperature effects

Various functionals assessing solvation free energies and enthalpies of the proton in methanol.

Determination of the Bridging Ligand in the Active Site of Tyrosinase

Tyrosinase is a type-3 copper enzyme that is widely distributed in plants, fungi, insects, and mammals. Developing high potent inhibitors against tyrosinase is of great interest in diverse fields including tobacco curing, food processing, bio-insecticides development, cosmetic development, and human healthcare-related research. In the crystal structure of Agaricus bisporus mushroom tyrosinase, there is an oxygen atom bridging the two copper ions in the active site. It is unclear whether the identity of this bridging oxygen is a water molecule or a hydroxide anion. In the present study, we theoretically determine the identity of this critical bridging oxygen by performing first-principles hybrid quantum mechanics/molecular mechanics/Poisson-Boltzmann-surface area (QM/MM-PBSA) calculations along with a thermodynamic cycle that aim to improve the accuracy. Our results show that the binding with water molecule is energy favored and the QM/MM-optimized structure is very close to the crystal structure, whereas the binding with hydroxide anions causes the increase of energy and significant structural changes of the active site, indicating that the identity of the bridging oxygen must be a water molecule rather than a hydroxide anion. The different binding behavior between water and hydroxide anions may explain why molecules with a carboxyl group or too many negative charges have lower inhibitory activity. In light of this, the design of high potent active inhibitors against tyrosinase should satisfy both the affinity to the copper ions and the charge neutrality of the entire molecule.

Smart, Innovative and Intelligent Technologies Used in Drug Designing

Smart and intelligent computational methods are essential nowadays for designing, manufacturing and optimizing new drugs. New and innovative computational tools and algorithms are consistently developed and applied for the development of novel therapeutic compounds in many research projects. Rapid developments in the architecture of computers have also provided complex calculations to be performed in a smart, intelligent and timely manner for desired quality outputs. Research groups worldwide are developing drug discovery platforms and innovative tools following smart manufacturing ideas using highly advanced biophysical, statistical and mathematical methods for accelerated discovery and analysis of smaller molecules. This chapter discusses novel innovative applications in drug discovery involving use of structure-based drug design which utilizes geometrical knowledge of the three-dimensional protein structures. It discusses statistical and physics based methods such as quantum mechanics and classical molecular dynamics which can also play a major role in improving the performance and in prediction of computational drug discovery. Lastly, the authors provide insights on recent developments in cloud computing with significant increase in smart and intelligent computational power thus allowing larger data sets to be analyzed simultaneously on multi processor cloud systems. Future directions for the research are outlined.

Pathways for Arene Oxidation in Non-Heme Diiron Enzymes: Lessons from Computational Studies on Benzoyl Coenzyme A Epoxidase

Oxygenation of aromatic rings using O₂ is catalyzed by several non-heme carboxylate-bridged diiron enzymes. In order to provide a general mechanistic description for these reactions, computational studies were carried out at the ONIOM(B3LYP/BP86/Amber) level on the non-heme diiron enzyme benzoyl coenzyme A epoxidase, BoxB. The calculations revealed four possible pathways for attacking the aromatic ring: (a) electrophilic (2e^-) attack by a bis(μ-oxo)-diiron(IV) species (Q pathway); (b) electrophilic (2e^-) attack via the σ* orbital of a μ-η²:η²-peroxo-diiron(III) intermediate (Pσ* pathway); (c) radical (1e^-) attack via the π*-orbital of a superoxo-diiron(II,III) species (Pπ* pathway); (d) radical (1e^-) attack of a partially quenched bis(μ-oxo)-diiron(IV) intermediate (Q' pathway). The results allowed earlier work of de Visser on olefin epoxidation by diiron complexes and QM-cluster studies of Liao and Siegbahn on BoxB to be put into a broader perspective. Parallels with epoxidation using organic peracids were also examined. Specifically for the BoxB enzyme, the Q pathway was found to be the most preferred, but the corresponding bis(μ-oxo)-diiron(IV) species is significantly destabilized and not expected to be directly observable. Epoxidation via the Pσ* pathway represents an energetically somewhat higher lying alternative; possible strategies for experimental discrimination are discussed. The selectivity toward epoxidation is shown to stem from a combination of inherent electronic properties of the thioacyl substituent and enzymatic constraints. Possible implications of the results for toluene monooxygenases are considered as well.

Binding-affinity predictions of HSP90 in the D3R Grand Challenge 2015 with docking, MM/GBSA, QM/MM, and free-energy simulations

We have estimated the binding affinity of three sets of ligands of the heat-shock protein 90 in the D3R grand challenge blind test competition. We have employed four different methods, based on five different crystal structures: first, we docked the ligands to the proteins with induced-fit docking with the Glide software and calculated binding affinities with three energy functions. Second, the docked structures were minimised in a continuum solvent and binding affinities were calculated with the MM/GBSA method (molecular mechanics combined with generalised Born and solvent-accessible surface area solvation). Third, the docked structures were re-optimised by combined quantum mechanics and molecular mechanics (QM/MM) calculations. Then, interaction energies were calculated with quantum mechanical calculations employing 970-1160 atoms in a continuum solvent, combined with energy corrections for dispersion, zero-point energy and entropy, ligand distortion, ligand solvation, and an increase of the basis set to quadruple-zeta quality. Fourth, relative binding affinities were estimated by free-energy simulations, using the multi-state Bennett acceptance-ratio approach. Unfortunately, the results were varying and rather poor, with only one calculation giving a correlation to the experimental affinities larger than 0.7, and with no consistent difference in the quality of the predictions from the various methods. For one set of ligands, the results could be strongly improved (after experimental data were revealed) if it was recognised that one of the ligands displaced one or two water molecules. For the other two sets, the problem is probably that the ligands bind in different modes than in the crystal structures employed or that the conformation of the ligand-binding site or the whole protein changes.

Computational Approaches to the Chemical Equilibrium Constant in Protein-ligand Binding

The physiological role played by protein-ligand recognition has motivated the development of several computational approaches to the ligand binding affinity. Some of them, termed rigorous, have a strong theoretical foundation but involve too much computation to be generally useful. Some others alleviate the computational burden by introducing strong approximations and/or empirical calibrations, which also limit their general use. Most importantly, there is no straightforward correlation between the predictive power and the level of approximation introduced. Here, we present a general framework for the quantitative interpretation of protein-ligand binding based on statistical mechanics. Within this framework, we re-derive self-consistently the fundamental equations of some popular approaches to the binding constant and pinpoint the inherent approximations. Our analysis represents a first step towards the development of variants with optimum accuracy/efficiency ratio for each stage of the drug discovery pipeline.

The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities

INTRODUCTION

The molecular mechanics energies combined with the Poisson-Boltzmann or generalized Born and surface area continuum solvation (MM/PBSA and MM/GBSA) methods are popular approaches to estimate the free energy of the binding of small ligands to biological macromolecules. They are typically based on molecular dynamics simulations of the receptor-ligand complex and are therefore intermediate in both accuracy and computational effort between empirical scoring and strict alchemical perturbation methods. They have been applied to a large number of systems with varying success.

AREAS COVERED

The authors review the use of MM/PBSA and MM/GBSA methods to calculate ligand-binding affinities, with an emphasis on calibration, testing and validation, as well as attempts to improve the methods, rather than on specific applications.

EXPERT OPINION

MM/PBSA and MM/GBSA are attractive approaches owing to their modular nature and that they do not require calculations on a training set. They have been used successfully to reproduce and rationalize experimental findings and to improve the results of virtual screening and docking. However, they contain several crude and questionable approximations, for example, the lack of conformational entropy and information about the number and free energy of water molecules in the binding site. Moreover, there are many variants of the method and their performance varies strongly with the tested system. Likewise, most attempts to ameliorate the methods with more accurate approaches, for example, quantum-mechanical calculations, polarizable force fields or improved solvation have deteriorated the results.

Density functional theory calculations on entire proteins for free energies of binding: application to a model polar binding site

In drug optimization calculations, the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method can be used to compute free energies of binding of ligands to proteins. The method involves the evaluation of the energy of configurations in an implicit solvent model. One source of errors is the force field used, which can potentially lead to large errors due to the restrictions in accuracy imposed by its empirical nature. To assess the effect of the force field on the calculation of binding energies, in this article we use large-scale density functional theory (DFT) calculations as an alternative method to evaluate the energies of the configurations in a "QM-PBSA" approach. Our DFT calculations are performed with a near-complete basis set and a minimal parameter implicit solvent model, within the self-consistent calculation, using the ONETEP program on protein-ligand complexes containing more than 2600 atoms. We apply this approach to the T4-lysozyme double mutant L99A/M102Q protein, which is a well-studied model of a polar binding site, using a set of eight small aromatic ligands. We observe that there is very good correlation between the MM and QM binding energies in vacuum but less so in the solvent. The relative binding free energies from DFT are more accurate than the ones from the MM calculations, and give markedly better agreement with experiment for six of the eight ligands. Furthermore, in contrast to MM-PBSA, QM-PBSA is able to correctly predict a nonbinder.

Computational insights into function and inhibition of fatty acid amide hydrolase

The Fatty Acid Amide Hydrolase (FAAH) enzyme is a membrane-bound serine hydrolase responsible for the deactivating hydrolysis of a family of naturally occurring fatty acid amides. FAAH is a critical enzyme of the endocannabinoid system, being mainly responsible for regulating the level of its main cannabinoid substrate anandamide. For this reason, pharmacological inhibition of FAAH, which increases the level of endogenous anandamide, is a promising strategy to cure a variety of diseases including pain, inflammation, and cancer. Much structural, mutagenesis, and kinetic data on FAAH has been generated over the last couple of decades. This has prompted several informative computational investigations to elucidate, at the atomic-level, mechanistic details on catalysis and inhibition of this pharmaceutically relevant enzyme. Here, we review how these computational studies - based on classical molecular dynamics, full quantum mechanics, and hybrid QM/MM methods - have clarified the binding and reactivity of some relevant substrates and inhibitors of FAAH. We also discuss the experimental implications of these computational insights, which have provided a thoughtful elucidation of the complex physical and chemical steps of the enzymatic mechanism of FAAH. Finally, we discuss how computations have been helpful for building structure-activity relationships of potent FAAH inhibitors.

Evidence for substrate binding-induced zwitterion formation in the catalytic Cys-His dyad of the SARS-CoV main protease

The coronavirus main protease (M(pro)) represents an attractive drug target for antiviral therapy of coronavirus (CoV) infections, including severe acute respiratory syndrome (SARS). The SARS-CoV M(pro) and related CoV proteases have several distinct features, such as an uncharged Cys-His catalytic dyad embedded in a chymotrypsin-like protease fold, that clearly separate these enzymes from archetypical cysteine proteases. To further characterize the catalytic system of CoV main proteases and to obtain information about improved inhibitors, we performed comprehensive simulations of the proton-transfer reactions in the SARS-CoV M(pro) active site that lead to the Cys(-)/His(+) zwitterionic state required for efficient proteolytic activity. Our simulations, comprising the free enzyme as well as substrate-enzyme and inhibitor-enzyme complexes, lead us to predict that zwitterion formation is fostered by substrate binding but not inhibitor binding. This indicates that M(pro) employs a substrate-induced catalytic mechanism that further enhances its substrate specificity. Our computational data are in line with available experimental results, such as X-ray geometries, measured pKa values, mutagenesis experiments, and the measured differences between the kinetic parameters of substrates and inhibitors. The data also provide an atomistic picture of the formerly postulated electrostatic trigger involved in SARS-CoV M(pro) activity. Finally, they provide information on how a specific microenvironment may finely tune the activity of M(pro) toward specific viral protein substrates, which is known to be required for efficient viral replication. Our simulations also indicate that the low inhibition potencies of known covalently interacting inhibitors may, at least in part, be attributed to insufficient fostering of the proton-transfer reaction. These findings suggest ways to achieve improved inhibitors.

A New Quantum Calibrated Force Field for Zinc-Protein Complex

A quantum calibrated polarizable-charge transfer force field (QPCT) has been proposed to accurately describe the interaction dynamics of zinc-protein complexes. The parameters of the QPCT force field were calibrated by quantum chemistry calculation and capture the polarization and charge transfer effect. QPCTs are validated by molecular dynamic simulation of the hydration shell of the zinc ion, five proteins containing the most common zinc-binding sites (ZnCys2His2, ZnCys3His1, ZnCys4, Zn2Cys6), as well as protein-ligand binding energy in zinc protein MMP3. The calculated results show excellent agreement with the experimental measurement and with results from QM/MM simulation, demonstrating that QPCT is accurate enough to maintain the correct structural integrity of the zinc binding pocket and provide accurate interaction dynamics of the zinc-residue complex. The current approach can also be extended to the study of interaction dynamics of other metal-containing proteins by recalibrating the corresponding parameters to the specific complexes.

Benchmark Study for the Cysteine-Histidine Proton Transfer Reaction in a Protein Environment: Gas Phase, COSMO, QM/MM Approaches

Proton transfer reactions are of crucial interest for the investigation of proteins. We have investigated the accuracy of commonly used quantum chemical methods for the description of proton transfer reactions in different environments (gas phase, COSMO, QM/MM) using the proton transfer between the catalytic dyad residues cysteine 145 and histidine 41 of SARS coronavirus main protease as a case study. The test includes thermodynamic, kinetic, and structural properties. The study comprises computationally demanding ab initio approaches (HF, CC2, MP2, SCS-CC2, SCS-MP2, CCSD(T)), popular density functional theories (BLYP, B3LYP, M06-2X), and semiempirical methods (MNDO/d, AM1, RM1, PM3, PM6). The approximated coupled cluster approach LCCSD(T) is taken as a reference method. We find that the robustness of the tested methods with respect to the environment correlates well with the level of theory. As an example HF, CC2, MP2, and their SCS variants show similar errors for gas phase, COSMO, or QM/MM computations. In contrast for semiempirical methods, the errors strongly diversify if one goes from gas phase to COSMO or QM/MM. Particular problems are observed for the recent semiempirical methods PM6 and RM1, which show the best performance for gas phase calculations but possess larger errors in conjunction with COSMO. Finally, a combination of SCS-MP2 and B3LYP or M06-2X allows reliable estimates about remaining errors.

Accurate Reaction Energies in Proteins Obtained by Combining QM/MM and Large QM Calculations

We here suggest and test a new method to obtain stable energies in proteins for charge-neutral reactions by running large quantum mechanical (QM) calculations on structures obtained by combined QM and molecular mechanics (QM/MM) geometry optimization on several snapshots from molecular dynamics simulations. As a test case, we use a proton transfer between a metal-bound cysteine residue and a second-sphere histidine residue in the active site of [Ni,Fe] hydrogenase, which has been shown to be very sensitive to the surroundings. We include in the QM calculations all residues within 4.5 Å of the active site, two capped residues on each side of the active-site residues, and all charged groups that are buried inside the protein, which for this enzyme includes three iron-sulfur clusters, in total, 930 atoms. These calculations are performed at the BP86/def2-SV(P) level, but the energies are then extrapolated to the B3LYP/def2-TZVP level with a smaller QM system, and zero-point energy, entropy, and thermal effects are added. We test three approaches to model the remaining atoms of the protein solvent, viz., by standard QM/MM approaches using either mechanical or electrostatic embedding or by using a continuum solvation model for the large QM systems. Quite encouragingly, the three approaches give the same results within 14 kJ/mol, and variations in the size of the QM system do not change the energies by more than 8 kJ/mol, provided that the QM/MM junctions are not moved closer to the QM system. The statistical precision for the average over 10 snapshots is 1-3 kJ/mol.

Reorganization energy for internal electron transfer in multicopper oxidases

We have calculated the reorganization energy for the intramolecular electron transfer between the reduced type 1 copper site and the peroxy intermediate of the trinuclear cluster in the multicopper oxidase CueO. The calculations are performed at the combined quantum mechanics and molecular mechanics (QM/MM) level, based on molecular dynamics simulations with tailored potentials for the two copper sites. We obtain a reorganization energy of 91-133 kJ/mol, depending on the theoretical treatment. The two Cu sites contribute by 12 and 22 kJ/mol to this energy, whereas the solvent contribution is 34 kJ/mol. The rest comes from the protein, involving small contributions from many residues. We have also estimated the energy difference between the two electron-transfer states and show that the reduction of the peroxy intermediate is exergonic by 43-87 kJ/mol, depending on the theoretical method. Both the solvent and the protein contribute to this energy difference, especially charged residues close to the two Cu sites. We compare these estimates with energies obtained from QM/MM optimizations and QM calculations in a vacuum and discuss differences between the results obtained at various levels of theory.

Reduction potentials and acidity constants of Mn superoxide dismutase calculated by QM/MM free-energy methods

We used two theoretical methods to estimate reduction potentials and acidity constants in Mn superoxide dismutase (MnSOD), namely combined quantum mechanical and molecular mechanics (QM/MM) thermodynamic cycle perturbation (QTCP) and the QM/MM-PBSA approach. In the latter, QM/MM energies are combined with continuum solvation energies calculated by solving the Poisson-Boltzmann equation (PB) or by the generalised Born approach (GB) and non-polar solvation energies calculated from the solvent-exposed surface area. We show that using the QTCP method, we can obtain accurate and precise estimates of the proton-coupled reduction potential for MnSOD, 0.30±0.01 V, which compares favourably with experimental estimates of 0.26-0.40 V. However, the calculated potentials depend strongly on the DFT functional used: The B3LYP functional gives 0.6 V more positive potentials than the PBE functional. The QM/MM-PBSA approach leads to somewhat too high reduction potentials for the coupled reaction and the results depend on the solvation model used. For reactions involving a change in the net charge of the metal site, the corresponding results differ by up to 1.3 V or 24 pK(a) units, rendering the QM/MM-PBSA method useless to determine absolute potentials. However, it may still be useful to estimate relative shifts, although the QTCP method is expected to be more accurate.

A QM/MM study of the binding of RAPTA ligands to cathepsin B

We have carried out quantum mechanical (QM) and QM/MM (combined QM and molecular mechanics) calculations, as well as molecular dynamics (MD) simulations to study the binding of a series of six RAPTA (Ru(II)-arene-1,3,5-triaza-7-phosphatricyclo-[3.3.1.1] decane) complexes with different arene substituents to cathepsin B. The recently developed QM/MM-PBSA approach (QM/MM combined with Poisson-Boltzmann solvent-accessible surface area solvation) has been used to estimate binding affinities. The QM calculations reproduce the antitumour activities of the complexes with a correlation coefficient (r (2)) of 0.35-0.86 after a conformational search. The QM/MM-PBSA method gave a better correlation (r (2) = 0.59) when the protein was fixed to the crystal structure, but more reasonable ligand structures and absolute binding energies were obtained if the protein was allowed to relax, indicating that the ligands are strained when the protein is kept fixed. In addition, the best correlation (r (2) = 0.80) was obtained when only the QM energies were used, which suggests that the MM and continuum solvation energies are not accurate enough to predict the binding of a charged metal complex to a charged protein. Taking into account the protein flexibility by means of MD simulations slightly improves the correlation (r (2) = 0.91), but the absolute energies are still too large and the results are sensitive to the details in the calculations, illustrating that it is hard to obtain stable predictions when full flexible protein is included in the calculations.

First Principles-Based Calculations of Free Energy of Binding: Application to Ligand Binding in a Self-Assembling Superstructure

The accurate prediction of ligand binding affinities to a protein remains a desirable goal of computational biochemistry. Many available methods use molecular mechanics (MM) to describe the system, however, MM force fields cannot fully describe the complex interactions involved in binding, specifically electron transfer and polarization. First principles approaches can fully account for these interactions, and with the development of linear-scaling first principles programs, it is now viable to apply first principles calculations to systems containing tens of thousands of atoms. In this paper, a quantum mechanical Poisson-Boltzmann surface area approach is applied to a model of a protein-ligand binding cavity, the "tennis ball" dimer. Results obtained from this approach demonstrate considerable improvement over conventional molecular mechanics Poisson-Boltzmann surface area due to the more accurate description of the interactions in the system. For the first principles calculations in this study, the linear-scaling density functional theory program ONETEP is used, allowing the approach to be applied to receptor-ligand complexes of pharmaceutical interest that typically include thousands of atoms.

Free-energy-based methods for binding profile determination in a congeneric series of CDK2 inhibitors

Free-energy pathway methods show great promise in computing the mode of action and the free energy profile associated with the binding of small molecules with proteins, but are generally very computationally demanding. Here we apply a novel approach based on metadynamics and path collective variables. We show that this combination is able to find an optimal reaction coordinate and the free energy profile of binding with explicit solvent and full flexibility, while minimizing human intervention and computational costs. We apply it to predict the binding affinity of a congeneric series of 5 CDK2 inhibitors. The predicted binding free energy profiles are in accordance with experiment.

Exploring the binding of inhibitors derived from tetrabromobenzimidazole to the CK2 protein using a QM/MM-PB/SA approach

We present an adaptation of the MM-PB/SA method for the estimation of binding free energies in protein-ligand complexes simulated with QM/MM molecular dynamics. The method is applied to understand the binding of a set of tetrabromobenzimidazole inhibitors to the CK2 protein. We find that the QM/MM interaction energy alone cannot always be used as a predictor of the binding affinity, and the inclusion of solvation effects via the PB/SA method is essential in getting reliable results. In agreement with experimental observations, we show that the van der Waals interactions are the driving force for the binding, while the electrostatic interactions orient these inhibitors in the CK2 active site. Additionally a per-residue energy decomposition analysis was applied to determine the individual contributions to the protein-inhibitor interaction. Based on these results, we hypothesize that the inclusion of a sufficiently large polar group on the tetrabromobenzimidazole skeleton could increase the binding affinity. The results show that the QM/MM-PB/SA method can be successfully employed to understand complicated structure-activity relationships and to design new inhibitors.