Binding affinities by alchemical perturbation using QM/MM with a large QM system and polarizable MM model

The most general way to improve the accuracy of binding‐affinity calculations for protein–ligand systems is to use quantum‐mechanical (QM) methods together with rigorous alchemical‐perturbation (AP) methods. We explore this approach by calculating the relative binding free energy of two synthetic disaccharides binding to galectin‐3 at a reasonably high QM level (dispersion‐corrected density functional theory with a triple‐zeta basis set) and with a sufficiently large QM system to include all short‐range interactions with the ligand (744–748 atoms). The rest of the protein is treated as a collection of atomic multipoles (up to quadrupoles) and polarizabilities. Several methods for evaluating the binding free energy from the 3600 QM calculations are investigated in terms of stability and accuracy. In particular, methods using QM calculations only at the endpoints of the transformation are compared with the recently proposed non‐Boltzmann Bennett acceptance ratio (NBB) method that uses QM calculations at several stages of the transformation. Unfortunately, none of the rigorous approaches give sufficient statistical precision. However, a novel approximate method, involving the direct use of QM energies in the Bennett acceptance ratio method, gives similar results as NBB but with better precision, ～3 kJ/mol. The statistical error can be further reduced by performing a greater number of QM calculations. © 2015 Wiley Periodicals, Inc.     

Convergence in the QM‐only and QM/MM modeling of enzymatic reactions: A case study for acetylene hydratase

We report systematic quantum mechanics‐only (QM‐only) and QM/molecular mechanics (MM) calculations on an enzyme‐catalyzed reaction to assess the convergence behavior of QM‐only and QM/MM energies with respect to the size of the chosen QM region. The QM and MM parts are described by density functional theory (typically B3LYP/def2‐SVP) and the CHARMM force field, respectively. Extending our previous work on acetylene hydratase with QM regions up to 157 atoms (Liao and Thiel, J. Chem. Theory Comput. 2012, 8, 3793), we performed QM/MM geometry optimizations with a QM region M4 composed of 408 atoms, as well as further QM/MM single‐point calculations with even larger QM regions up to 657 atoms. A charge deletion analysis was conducted for the previously used QM/MM model ( M3a , with a QM region of 157 atoms) to identify all MM residues with strong electrostatic contributions to the reaction energetics (typically more than 2 kcal/mol), which were then included in M4 . QM/MM calculations with this large QM region M4 lead to the same overall mechanism as the previous QM/MM calculations with M3a , but there are some variations in the relative energies of the stationary points, with a mean absolute deviation (MAD) of 2.7 kcal/mol. The energies of the two relevant transition states are close to each other at all levels applied (typically within 2 kcal/mol), with the first (second) one being rate‐limiting in the QM/MM calculations with M3a ( M4 ). QM‐only gas‐phase calculations give a very similar energy profile for QM region M4 (MAD of 1.7 kcal/mol), contrary to the situation for M3a where we had previously found significant discrepancies between the QM‐only and QM/MM results (MAD of 7.9 kcal/mol). Extension of the QM region beyond M4 up to M7 (657 atoms) leads to only rather small variations in the relative energies from single‐point QM‐only and QM/MM calculations (MAD typically about 1–2 kcal/mol). In the case of acetylene hydratase, a model with 408 QM atoms thus seems sufficient to achieve convergence in the computed relative energies to within 1–2 kcal/mol.Copyright © 2013 Wiley Periodicals, Inc.     

Predicting hydration free energies with a hybrid QM/MM approach: an evaluation of implicit and explicit solvation models in SAMPL4

The correct representation of solute-water interactions is essential for the accurate simulation of most biological phenomena. Several highly accurate quantum methods are available to deal with solvation by using both implicit and explicit solvents. So far, however, most evaluations of those methods were based on a single conformation, which neglects solute entropy. Here, we present the first test of a novel approach to determine hydration free energies that uses molecular mechanics (MM) to sample phase space and quantum mechanics (QM) to evaluate the potential energies. Free energies are determined by using re-weighting with the Non-Boltzmann Bennett (NBB) method. In this context, the method is referred to as QM-NBB. Based on snapshots from MM sampling and accounting for their correct Boltzmann weight, it is possible to obtain hydration free energies that incorporate the effect of solute entropy. We evaluate the performance of several QM implicit solvent models, as well as explicit solvent QM/MM for the blind subset of the SAMPL4 hydration free energy challenge. While classical free energy simulations with molecular dynamics give root mean square deviations (RMSD) of 2.8 and 2.3 kcal/mol, the hybrid approach yields an improved RMSD of 1.6 kcal/mol. By selecting an appropriate functional and basis set, the RMSD can be reduced to 1 kcal/mol for calculations based on a single conformation. Results for a selected set of challenging molecules imply that this RMSD can be further reduced by using NBB to reweight MM trajectories with the SMD implicit solvent model.     

Fragment‐based quantum mechanical calculation of protein–protein binding affinities

The electrostatically embedded generalized molecular fractionation with conjugate caps (EE‐GMFCC) method has been successfully utilized for efficient linear‐scaling quantum mechanical (QM) calculation of protein energies. In this work, we applied the EE‐GMFCC method for calculation of binding affinity of Endonuclease colicin–immunity protein complex. The binding free energy changes between the wild‐type and mutants of the complex calculated by EE‐GMFCC are in good agreement with experimental results. The correlation coefficient (R) between the predicted binding energy changes and experimental values is 0.906 at the B3LYP/6‐31G*‐D level, based on the snapshot whose binding affinity is closest to the average result from the molecular mechanics/Poisson–Boltzmann surface area (MM/PBSA) calculation. The inclusion of the QM effects is important for accurate prediction of protein–protein binding affinities. Moreover, the self‐consistent calculation of PB solvation energy is required for accurate calculations of protein–protein binding free energies. This study demonstrates that the EE‐GMFCC method is capable of providing reliable prediction of relative binding affinities for protein–protein complexes. © 2018 Wiley Periodicals, Inc.     

Binding free energies in the SAMPL6 octa-acid host–guest challenge calculated with MM and QM methods

We have estimated free energies for the binding of eight carboxylate ligands to two variants of the octa-acid deep-cavity host in the SAMPL6 blind-test challenge (with or without endo methyl groups on the four upper-rim benzoate groups, OAM and OAH, respectively). We employed free-energy perturbation (FEP) for relative binding energies at the molecular mechanics (MM) and the combined quantum mechanical (QM) and MM (QM/MM) levels, the latter obtained with the reference-potential approach with QM/MM sampling for the MM → QM/MM FEP. The semiempirical QM method PM6-DH+ was employed for the ligand in the latter calculations. Moreover, binding free energies were also estimated from QM/MM optimised structures, combined with COSMO-RS estimates of the solvation energy and thermostatistical corrections from MM frequencies. They were performed at the PM6-DH+ level of theory with the full host and guest molecule in the QM system (and also four water molecules in the geometry optimisations) for 10–20 snapshots from molecular dynamics simulations of the complex. Finally, the structure with the lowest free energy was recalculated using the dispersion-corrected density-functional theory method TPSS-D3, for both the structure and the energy. The two FEP approaches gave similar results (PM6-DH+/MM slightly better for OAM), which were among the five submissions with the best performance in the challenge and gave the best results without any fit to data from the SAMPL5 challenge, with mean absolute deviations (MAD) of 2.4–5.2 kJ/mol and a correlation coefficient (R²) of 0.77–0.93. This is the first time QM/MM approaches give binding free energies that are competitive to those obtained with MM for the octa-acid host. The QM/MM-optimised structures gave somewhat worse performance (MAD?=?3–8 kJ/mol and R²?=?0.1–0.9), but the results were improved compared to previous studies of this system with similar methods.     

Performance assessment of semiempirical molecular orbital methods in describing halogen bonding: quantum mechanical and quantum mechanical/molecular mechanical-molecular dynamics study

The performance of semiempirical molecular-orbital methods--MNDO, MNDO-d, AM1, RM1, PM3 and PM6--in describing halogen bonding was evaluated, and the results were compared with molecular mechanical (MM) and quantum mechanical (QM) data. Three types of performance were assessed: (1) geometrical optimizations and binding energy calculations for 27 halogen-containing molecules complexed with various Lewis bases (Two of the tested methods, AM1 and RM1, gave results that agree with the QM data.); (2) charge distribution calculations for halobenzene molecules, determined by calculating the solvation free energies of the molecules relative to benzene in explicit and implicit generalized Born (GB) solvents (None of the methods gave results that agree with the experimental data.); and (3) appropriateness of the semiempirical methods in the hybrid quantum-mechanical/molecular-mechanical (QM/MM) scheme, investigated by studying the molecular inhibition of CK2 protein by eight halobenzimidazole and -benzotriazole derivatives using hybrid QM/MM molecular-dynamics (MD) simulations with the inhibitor described at the QM level by the AM1 method and the rest of the system described at the MM level. The pure MM approach with inclusion of an extra point of positive charge on the halogen atom approach gave better results than the hybrid QM/MM approach involving the AM1 method. Also, in comparison with the pure MM-GBSA (generalized Born surface area) binding energies and experimental data, the calculated QM/MM-GBSA binding energies of the inhibitors were improved by replacing the G(GB,QM/MM) solvation term with the corresponding G(GB,MM) term.     

The QM‐MM interface for CHARMM‐deMon

We present a new QM/MM interface for fast and efficient simulations of organic and biological molecules. The CHARMM/deMon interface has been developed and tested to perform minimization and atomistic simulations for multi‐particle systems. The current features of this QM/MM interface include readability for molecular dynamics, tested compatibility with Free Energy Perturbation simulations (FEP) using the dual topology/single coordinate method. The current coupling scheme uses link atoms, but further extensions of the code to incorporate other available schemes are planned. We report the performance of different levels of theory for the treatment of the QM region, while the MM region was represented by a classical force‐field (CHARMM27) or a polarizable force‐field based on a simple Drude model. The current QM/MM implementation can be coupled to the dual‐thermostat method and the VV2 integrator to run molecular dynamics simulations. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Use of a QM/MM-based FEP method to evaluate the anomalous hydration behavior of simple alkyl amines and amides: application to the design of FBPase inhibitors for the treatment of type-2 diabetes

Standard molecular mechanics (MM) force fields predict a nearly linear decrease in hydration free energy with each successive addition of a methyl group to ammonia or acetamide, whereas a nonadditive relationship is observed experimentally. In contrast, the non-additive hydration behavior is reproduced directly using a quantum mechanics (QM)/MM-based free-energy perturbation (FEP) method wherein the solute partial atomic charges are updated at every window. Decomposing the free energies into electrostatic and van der Waals contributions and comparing the results with the corresponding free energies obtained using a conventional FEP method and a QM/MM method wherein the charges are not updated suggests that inaccuracies in the electrostatic free energies are the primary reason for the inability of the conventional FEP method to predict the experimental findings. The QM/MM-based FEP method was subsequently used to evaluate inhibitors of the diabetes drug target fructose-1,6-bisphosphatase adenosine 5'-monophosphate and 6-methylamino purine riboside 5'-monophosphate. The predicted relative binding free energy was consistent with the experimental findings, whereas the relative binding free energy predicted using the conventional FEP method differed from the experimental finding by an amount consistent with the overestimated relative solvation free energies calculated for alkylamines. Accordingly, the QM/MM-based FEP method offers potential advantages over conventional FEP methods, including greater accuracy and reduced user input. Moreover, since drug candidates often contain either functionality that is inadequately treated by MM (e.g., simple alkylamines and alkylamides) or new molecular scaffolds that require time-consuming development of MM parameters, these advantages could enable future automation of FEP calculations as well as greatly increase the use and impact of FEP calculations in drug discovery.     

Methodological aspects of QM/MM calculations: A case study on matrix metalloproteinase‐2

We address methodological issues in quantum mechanics/molecular mechanics (QM/MM) calculations on a zinc‐dependent enzyme. We focus on the first stage of peptide bond cleavage by matrix metalloproteinase‐2 (MMP‐2), that is, the nucleophilic attack of the zinc‐coordinating water molecule on the carbonyl carbon atom of the scissile fragment of the substrate. This step is accompanied by significant charge redistribution around the zinc cation, bond cleavage, and bond formation. We vary the size and initial geometry of the model system as well as the computational protocol to demonstrate the influence of these choices on the results obtained. We present QM/MM potential energy profiles for a set of snapshots randomly selected from QM/MM‐based molecular dynamics simulations and analyze the differences in the computed profiles in structural terms. Since the substrate in MMP‐2 is located on the protein surface, we investigate the influence of the thickness of the water layer around the enzyme on the QM/MM energy profile. Thin water layers (0–2 Å) give unrealistic results because of structural reorganizations in the active‐site region at the protein surface. A 12 Å water layer appears to be sufficient to capture the effect of the solvent; the corresponding QM/MM energy profile is very close to that obtained from QM/MM/SMBP calculations using the solvent macromolecular boundary potential (SMBP). We apply the optimized computational protocol to explain the origin of the different catalytic activity of the Glu116Asp mutant: the energy barrier for the first step is higher, which is rationalized on structural grounds. © 2016 Wiley Periodicals, Inc.     

Mechanism of proteolysis in matrix metalloproteinase‐2 revealed by QM/MM modeling

The mechanism of enzymatic peptide hydrolysis in matrix metalloproteinase‐2 (MMP‐2) was studied at atomic resolution through quantum mechanics/molecular mechanics (QM/MM) simulations. An all‐atom three‐dimensional molecular model was constructed on the basis of a crystal structure from the Protein Data Bank (ID: 1QIB), and the oligopeptide Ace‐Gln‐Gly～Ile‐Ala‐Gly‐Nme was considered as the substrate. Two QM/MM software packages and several computational protocols were employed to calculate QM/MM energy profiles for a four‐step mechanism involving an initial nucleophilic attack followed by hydrogen bond rearrangement, proton transfer, and C? N bond cleavage. These QM/MM calculations consistently yield rather low overall barriers for the chemical steps, in the range of 5–10 kcal/mol, for diverse QM treatments (PBE0, B3LYP, and BB1K density functionals as well as local coupled cluster treatments) and two MM force fields (CHARMM and AMBER). It, thus, seems likely that product release is the rate‐limiting step in MMP‐2 catalysis. This is supported by an exploration of various release channels through QM/MM reaction path calculations and steered molecular dynamics simulations. © 2015 Wiley Periodicals, Inc.     

The CHARMM–TURBOMOLE interface for efficient and accurate QM/MM molecular dynamics,free energies,and excited state properties

The quantum mechanical (QM)/molecular mechanical (MM) interface between Chemistry at HARvard Molecular Mechanics (CHARMM) and TURBOMOLE is described. CHARMM provides an extensive set of simulation algorithms, like molecular dynamics (MD) and free energy perturbation, and support for mature nonpolarizable and Drude polarizable force fields. TURBOMOLE provides fast QM calculations using density functional theory or wave function methods and excited state properties. CHARMM–TURBOMOLE is well‐suited for extended QM/MM MD simulations using first principles methods with large (triple‐ζ) basis sets. We demonstrate these capabilities with a QM/MM simulation of Mg²⁺(aq), where the MM outer sphere water molecules are represented using the SWM4‐NDP Drude polarizable force field and the ion and inner coordination sphere are represented using QM PBE, PBE0, and MP2 methods. The relative solvation free energies of Mg²⁺ and Zn²⁺ were calculated using thermodynamic integration. We also demonstrate the features for excited state properties. We calculate the time‐averaged solution absorption spectrum of indole, the emission spectrum of the indole excited state, and the electronic circular dichroism spectrum of an oxacepham. © 2014 Wiley Periodicals, Inc.     

Contributions of pauli repulsions to the energetics and physical properties computed in QM/MM methods

Conventional combined quantum mechanical/molecular mechanical (QM/MM) methods lack explicit treatment of Pauli repulsions between the quantum‐mechanical and molecular‐mechanical subsystems. Instead, classical Lennard‐Jones (LJ) potentials between QM and MM nuclei are used to model electronic Pauli repulsion and long‐range London dispersion, despite the fact that the latter two are inherently of quantum nature. Use of the simple LJ potential in QM/MM methods can reproduce minimal geometries and energies of many molecular clusters reasonably well, as compared to full QM calculations. However, we show here that the LJ potential cannot correctly describe subtle details of the electron density of the QM subsystem because of the neglect of Pauli repulsions between the QM and MM subsystems. The inaccurate electron density subsequently affects the calculation of electronic and magnetic properties of the QM subsystem. To explicitly consider Pauli interactions with QM/MM methods, we propose a method to use empirical effective potentials on the MM atoms. The test case of the binding energy and magnetic properties of a water dimer shows promising results for the general application of effective potentials to mimic Pauli repulsions in QM/MM calculations. © 2013 Wiley Periodicals, Inc.     

LICHEM: A QM/MM program for simulations with multipolar and polarizable force fields

We introduce an initial implementation of the LICHEM software package. LICHEM can interface with Gaussian, PSI4, NWChem, TINKER, and TINKER–HP to enable QM/MM calculations using multipolar/polarizable force fields. LICHEM extracts forces and energies from unmodified QM and MM software packages to perform geometry optimizations, single‐point energy calculations, or Monte Carlo simulations. When the QM and MM regions are connected by covalent bonds, the pseudo‐bond approach is employed to smoothly transition between the QM region and the polarizable force field. A series of water clusters and small peptides have been employed to test our initial implementation. The results obtained from these test systems show the capabilities of the new software and highlight the importance of including explicit polarization. © 2016 Wiley Periodicals, Inc.     

Computational methods for the study of enzymic reaction mechanisms III: a perturbation plus QM/MM approach for calculating relative free energies of protonation

We describe a coupling parameter, that is, perturbation, approach to effectively create and annihilate atoms in the quantum mechanical Hamiltonian within the closed shell restricted Hartree-Fock formalism. This perturbed quantum mechanical atom (PQA) method is combined with molecular mechanics (MM) methods (PQA/MM) within a molecular dynamics simulation, to model the protein environment (MM region) effects that also make a contribution to the overall free energy change. Using the semiempirical PM3 method to model the QM region, the application of this PQA/MM method is illustrated by calculation of the relative protonation free energy of the conserved OD2 (Asp27) and the N5 (dihydrofolate) proton acceptor sites in the active site of Escherichia coli dihydrofolate reductase (DHFR) with the bound nicotinamide adenine dinucleotide phosphate (NADPH) cofactor. For a number of choices for the QM region, the relative protonation free energy was calculated as the sum of contributions from the QM region and the interaction between the QM and MM regions via the thermodynamic integration (TI) method. The results demonstrate the importance of including the whole substrate molecule in the QM region, and the overall protein (MM) environment in determining the relative stabilities of protonation sites in the enzyme active site. The PQA/MM free energies obtained by TI were also compared with those estimated by a less computationally demanding nonperturbative method based on the linear response approximation (LRA). For some choices of QM region, the total free energies calculated using the LRA method were in very close agreement with the PQA/MM values. However, the QM and QM/MM component free energies were found to differ significantly between the two methods.     

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 ²) of 0.35–0.86 after a conformational search. The QM/MM-PBSA method gave a better correlation (r ² = 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 ² = 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 ² = 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.     

Development of massive multilevel molecular dynamics simulation program,platypus (PLATform for dYnamic protein unified simulation), for the elucidation of protein functions

A massively parallel program for quantum mechanical‐molecular mechanical (QM/MM) molecular dynamics simulation, called Platypus (PLATform for dYnamic Protein Unified Simulation), was developed to elucidate protein functions. The speedup and the parallelization ratio of Platypus in the QM and QM/MM calculations were assessed for a bacteriochlorophyll dimer in the photosynthetic reaction center (DIMER) on the K computer, a massively parallel computer achieving 10 PetaFLOPs with 705,024 cores. Platypus exhibited the increase in speedup up to 20,000 core processors at the HF/cc‐pVDZ and B3LYP/cc‐pVDZ, and up to 10,000 core processors by the CASCI(16,16)/6‐31G** calculations. We also performed excited QM/MM‐MD simulations on the chromophore of Sirius (SIRIUS) in water. Sirius is a pH‐insensitive and photo‐stable ultramarine fluorescent protein. Platypus accelerated on‐the‐fly excited‐state QM/MM‐MD simulations for SIRIUS in water, using over 4000 core processors. In addition, it also succeeded in 50‐ps (200,000‐step) on‐the‐fly excited‐state QM/MM‐MD simulations for the SIRIUS in water. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.