QM/MM-PBSA method to estimate free energies for reactions in proteins

We have developed a method to estimate free energies of reactions in proteins, called QM/MM-PBSA. It estimates the internal energy of the reactive site by quantum mechanical (QM) calculations, whereas bonded, electrostatic, and van der Waals interactions with the surrounding protein are calculated at the molecular mechanics (MM) level. The electrostatic part of the solvation energy of the reactant and the product is estimated by solving the Poisson-Boltzmann (PB) equation, and the nonpolar part of the solvation energy is estimated from the change in solvent-accessible surface area (SA). Finally, the change in entropy is estimated from the vibrational frequencies. We test this method for five proton-transfer reactions in the active sites of [Ni,Fe] hydrogenase and copper nitrite reductase. We show that QM/MM-PBSA reproduces the results of a strict QM/MM free-energy perturbation method with a mean absolute deviation (MAD) of 8-10 kJ/mol if snapshots from molecular dynamics simulations are used and 4-14 kJ/mol if a single QM/MM structure is used. This is appreciably better than the original QM/MM results or if the QM energies are supplemented with a point-charge model, a self-consistent reaction field, or a PB model of the protein and the solvent, which give MADs of 22-36 kJ/mol for the same test set.     

MM and QM/MM Modeling of Threonyl-tRNA Synthetase: Model Testing and Simulations

Aminoacyl-tRNA synthetases are centrally important enzymes in protein synthesis. We have investigated threonyl-tRNA synthetase from E. coli, complexed with reactants, using molecular mechanics and combined quantum mechanical/molecular mechanical (QM/MM) techniques. These modeling methods have the potential to provide molecular level understanding of enzyme catalytic processes. Modeling of this enzyme presents a number of challenges. The procedure of system preparation and testing is described in detail. For example, the number of metal ions at the active site, and their positions, were investigated. Molecular dynamics simulations suggest that the system is most stable when it contains only one magnesium ion, and the zinc ion is removed. Two different QM/MM methods were tested in models based on the findings of MM molecular dynamics simulations. AM1/CHARMM calculations resulted in unrealistic structures for the phosphates in this system. This is apparently due to an error of AM1. PM3/CHARMM calculations proved to be more suitable for this enzyme system. These results will provide a useful basis for future modeling investigations of the enzyme mechanism and dynamics.     

Fragmentation-based QM/MM simulations: length dependence of chain dynamics and hydrogen bonding of polyethylene oxide and polyethylene in aqueous solutions

Molecular fragmentation quantum mechanics (QM) calculations have been combined with molecular mechanics (MM) to construct the fragmentation QM/MM method for simulations of dilute solutions of macromolecules. We adopt the electrostatics embedding QM/MM model, where the low-cost generalized energy-based fragmentation calculations are employed for the QM part. Conformation energy calculations, geometry optimizations, and Born-Oppenheimer molecular dynamics simulations of poly(ethylene oxide), PEO(n) (n = 6-20), and polyethylene, PE(n) ( n = 9-30), in aqueous solution have been performed within the framework of both fragmentation and conventional QM/MM methods. The intermolecular hydrogen bonding and chain configurations obtained from the fragmentation QM/MM simulations are consistent with the conventional QM/MM method. The length dependence of chain conformations and dynamics of PEO and PE oligomers in aqueous solutions is also investigated through the fragmentation QM/MM molecular dynamics simulations.     

Reaction path potential for complex systems derived from combined ab initio quantum mechanical and molecular mechanical calculations

Combined ab initio quantum mechanical and molecular mechanical calculations have been widely used for modeling chemical reactions in complex systems such as enzymes, with most applications being based on the determination of a minimum energy path connecting the reactant through the transition state to the product in the enzyme environment. However, statistical mechanics sampling and reaction dynamics calculations with a combined ab initio quantum mechanical (QM) and molecular mechanical (MM) potential are still not feasible because of the computational costs associated mainly with the ab initio quantum mechanical calculations for the QM subsystem. To address this issue, a reaction path potential energy surface is developed here for statistical mechanics and dynamics simulation of chemical reactions in enzymes and other complex systems. The reaction path potential follows the ideas from the reaction path Hamiltonian of Miller, Handy and Adams for gas phase chemical reactions but is designed specifically for large systems that are described with combined ab initio quantum mechanical and molecular mechanical methods. The reaction path potential is an analytical energy expression of the combined quantum mechanical and molecular mechanical potential energy along the minimum energy path. An expansion around the minimum energy path is made in both the nuclear and the electronic degrees of freedom for the QM subsystem internal energy, while the energy of the subsystem described with MM remains unchanged from that in the combined quantum mechanical and molecular mechanical expression and the electrostatic interaction between the QM and MM subsystems is described as the interaction of the MM charges with the QM charges. The QM charges are polarizable in response to the changes in both the MM and the QM degrees of freedom through a new response kernel developed in the present work. The input data for constructing the reaction path potential are energies, vibrational frequencies, and electron density response properties of the QM subsystem along the minimum energy path, all of which can be obtained from the combined quantum mechanical and molecular mechanical calculations. Once constructed, it costs much less for its evaluation. Thus, the reaction path potential provides a potential energy surface for rigorous statistical mechanics and reaction dynamics calculations of complex systems. As an example, the method is applied to the statistical mechanical calculations for the potential of mean force of the chemical reaction in triosephosphate isomerase.     

How does the cAMP-dependent protein kinase catalyze the phosphorylation reaction: an ab initio QM/MM study

We have carried out density functional theory QM/MM calculations on the catalytic subunit of cAMP-dependent protein kinase (PKA). The QM/MM calculations indicate that the phosphorylation reaction catalyzed by PKA is mainly dissociative, and Asp166 serves as the catalytic base to accept the proton delivered by the substrate peptide. Among the key interactions in the active site, the Mg(2+) ions, glycine rich loop, and Lys72 are found to stabilize the transition state through electrostatic interactions. On the other hand, Lys168, Asn171, Asp184, and the conserved waters bound to Mg(2+) ions do not directly contribute to lower the energy barrier of the phosphorylation reaction, and possible roles for these residues are proposed. The QM/MM calculations with different QM/MM partition schemes or different initial structures yield consistent results. In addition, we have carried out 12 ns molecular dynamics simulations on both wild type and K168A mutated PKA, respectively, to demonstrate that the catalytic role of Lys168 is to keep ATP and substrate peptide in the near-attack reactive conformation.     

Quantum mechanically derived AMBER-compatible heme parameters for various states of the cytochrome P450 catalytic cycle

Molecular mechanics (MM) methods are computationally affordable tools for screening chemical libraries of novel compounds for sites of P450 metabolism. One challenge for MM methods has been the absence of a consistent and transferable set of parameters for the heme within the P450 active site. Experimental data indicate that mammalian P450 enzymes vary greatly in the size, architecture, and plasticity of their active sites. Thus, obtaining X-ray-based geometries for the development of accurate MM parameters for the major classes of hepatic P450 remains a daunting task. Our previous work with preliminary gas-phase quantum mechanics (QM)-derived atomic partial charges greatly improved the accuracy of docking studies of raloxifene to CYP3A4. We have therefore developed and tested a consistent set of transferable MM parameters based on gas-phase QM calculations of two model systems of the heme-a truncated (T-HM) and a full (F-HM) for four states of the P450 catalytic cycle. Our results indicate that the use of the atomic partial charges from the F-HM further improves the accuracy of docked predictions for raloxifene to CYP3A4. Different patterns for substrate docking are also observed depending on the choice of heme model and state. Newly parameterized heme models are tested in implicit and explicitly solvated MD simulations in the absence and presence of enzyme structures, for CYP3A4, and appear to be stable on the nanosecond simulation timescale. The new force field for the various heme states may aid the community for simulations of P450 enzymes and other heme-containing enzymes.     

The elusive 5'-deoxyadenosyl radical in coenzyme-B12-mediated reactions

Vitamin B(12) and its biologically active counterparts possess the only examples of carbon-cobalt bonds in living systems. The role of such motifs as radical reservoirs has potential application in future catalytic and electronic nanodevices. To fully understand radical generation in coenzyme B(12) (dAdoCbl)-dependent enzymes, however, major obstacles still need to be overcome. In this work, we have used Car-Parrinello molecular dynamics (CPMD) simulations, in a mixed quantum mechanics/molecular mechanics (QM/MM) framework, to investigate the initial stages of the methylmalonyl-CoA-mutase-catalyzed reaction. We demonstrate that the 5'-deoxyadenosyl radical (dAdo(?)) exists as a distinct entity in this reaction, consistent with the results of extensive experimental and some previous theoretical studies. We report free energy calculations and first-principles trajectories that help understand how B(12) enzymes catalyze coenzyme activation and control highly reactive radical intermediates.     

A comparative QM/MM simulation study of the reaction mechanisms of human and Plasmodium falciparum HG(X)PRTases

QM/MM hybrid potential free-energy simulations are performed to compare the reaction mechanisms of human hypoxanthine guanine phosphoribosyl transferase (HGPRTase) and the corresponding enyzme from Plasmodium falciparum (Pf), hypoxanthine guanine xanthine phosphoribosyl transferase (HGXPRTase). These enzymes share 44% of sequence identity but display very different affinities for xanthine. The calculations show that in both enzymes phosphoribosyl transfer proceeds via a dissociative mechanism from an anionic form of the substrate. Nevertheless, there are significant differences in the geometries of critical structures along the reaction paths which it may be possible to exploit for the design of specific inhibitors against the Pf enzyme.     

Development of a quantum mechanics-based free-energy perturbation method: use in the calculation of relative solvation free energies

Free-energy perturbation (FEP) is considered the most accurate computational method for calculating relative solvation and binding free-energy differences. Despite some success in applying FEP methods to both drug design and lead optimization, FEP calculations are rarely used in the pharmaceutical industry. One factor limiting the use of FEP is its low throughput, which is attributed in part to the dependence of conventional methods on the user's ability to develop accurate molecular mechanics (MM) force field parameters for individual drug candidates and the time required to complete the process. In an attempt to find an FEP method that could eventually be automated, we developed a method that uses quantum mechanics (QM) for treating the solute, MM for treating the solute surroundings, and the FEP method for computing free-energy differences. The thread technique was used in all transformations and proved to be essential for the successful completion of the calculations. Relative solvation free energies for 10 structurally diverse molecular pairs were calculated, and the results were in close agreement with both the calculated results generated by conventional FEP methods and the experimentally derived values. While considerably more CPU demanding than conventional FEP methods, this method (QM/MM-based FEP) alleviates the need for development of molecule-specific MM force field parameters and therefore may enable future automation of FEP-based calculations. Moreover, calculation accuracy should be improved over conventional methods, especially for calculations reliant on MM parameters derived in the absence of experimental data.     

An interface between the universal force field and the effective fragment potential method

In order to properly describe reactions in heterogeneous catalyst systems, the reactants, solvent, and bulk effects of the surface must be taken into account. Embedded-cluster QM (quantum mechanics)/MM (molecular mechanics) methods can treat reactions on surfaces (the gas-surface interface), and the effective fragment potential method (EFP) can accurately treat the solvent effects on reactions (the gas-liquid interface). In order to create a QM/MM/EFP hybrid method for treatment of heterogeneous catalytic systems in the presence of a solvent (the liquid-surface interface), an EFP-MM interaction potential has been developed. Example calculations on small clusters of silica and water have been carried out.     

QM/MM analysis suggests that Alkaline Phosphatase (AP) and nucleotide pyrophosphatase/phosphodiesterase slightly tighten the transition state for phosphate diester hydrolysis relative to solution: implication for catalytic promiscuity in the AP superfamily

Several members of the Alkaline Phosphatase (AP) superfamily exhibit a high level of catalytic proffciency and promiscuity in structurally similar active sites. A thorough characterization of the nature of transition state for different substrates in these enzymes is crucial for understanding the molecular mechanisms that govern those remarkable catalytic properties. In this work, we study the hydrolysis of a phosphate diester, MpNPP(-), in solution, two experimentally well-characterized variants of AP (R166S AP, R166S/E322Y AP) and wild type Nucleotide pyrophosphatase/phosphodiesterase (NPP) by QM/MM calculations in which the QM method is an approximate density functional theory previously parametrized for phosphate hydrolysis (SCC-DFTBPR). The general agreements found between these calculations and available experimental data for both solution and enzymes support the use of SCC-DFTBPR/MM for a semiquantitative analysis of the catalytic mechanism and nature of transition state in AP and NPP. Although phosphate diesters are cognate substrates for NPP but promiscuous substrates for AP, the calculations suggest that their hydrolysis reactions catalyzed by AP and NPP feature similar synchronous transition states that are slightly tighter in nature compared to that in solution, due in part to the geometry of the bimetallic zinc motif. Therefore, this study provides the first direct computational support to the hypothesis that enzymes in the AP superfamily catalyze cognate and promiscuous substrates via similar transition states to those in solution. Our calculations do not support the finding of recent QM/MM studies by López-Canut and co-workers, who suggested that the same diester substrate goes through a much looser transition state in NPP/AP than in solution, a result likely biased by the large structural distortion of the bimetallic zinc site in their simulations. Finally, our calculations for different phosphate diester orientations and phosphorothioate diesters highlight that the interpretation of thio-substitution experiments is not always straightforward.     

酶催化过程的全程模拟

酶催化包括底物到活性区的输运、选择催化化学反应及产物释放等复杂过程,由于复杂的蛋白质环境效应,任一化学和非化学过程都有可能是决定酶活性的关键步骤。为了全面认识酶催化活性,我们对几类酶催化过程进行了广泛的组合量子/分子力学(QM/MM)和经典分子力学(MM)动力学模拟(MD)研究,详细地讨论了整个酶催化过程的分子机制、关键残基的作用和蛋白质环境效应,丰富了对酶催化活性的认识。随着多尺度模型和计算模拟方法的进一步完善与发展,有望实现超大复杂生物酶催化过程的全程模拟研究,为酶工程领域的相关研究提供支持。     

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.     

Cope elimination: elucidation of solvent effects from QM/MM simulations

The Cope elimination reactions for threo- and erythro-N,N-dimethyl-3-phenyl-2-butylamine oxide have been investigated using QM/MM calculations in water, THF, and DMSO. The aprotic solvents provide up to million-fold rate accelerations. The effects of solvation on the reactants, transition structures, and rates of reaction are elucidated here using two-dimensional potentials of mean force (PMF) derived from free-energy perturbation calculations in Monte Carlo simulations (MC/FEP). The resultant free energies of activation in solution are in close agreement with experiment. Ab initio calculations at the MP2/6-311+G-(2d,p) level using the PCM continuum solvent model were also carried out; however, only the QM/MM methodology was able to reproduce the large rate increases in proceeding from water to the dipolar aprotic solvents. Solute-solvent interaction energies and radial distribution functions are also analyzed and show that poorer solvation of the reactant in the aprotic solvents is primarily responsible for the observed rate enhancements. It is found that the amine oxide oxygen is the acceptor of three hydrogen bonds from water molecules for the reactant but only one to two weaker ones at the transition state. The overall quantitative success of the computations supports the present QM/MM/MC approach, featuring PDDG/PM3 as the QM method.     

On the effect of a variation of the force field, spatial boundary condition and size of the QM region in QM/MM MD simulations

During the past years, the use of combined quantum-classical, QM/MM, methods for the study of complex biomolecular processes, such as enzymatic reactions and photocycles, has increased considerably. The quality of the results obtained from QM/MM calculations is largely dependent on five aspects to be considered when setting up a molecular model: the QM Hamiltonian, the MM Hamiltonian or force field, the boundary and coupling between the QM and MM regions, the size of the QM region and the boundary condition for the MM region. In this study, we systematically investigate the influence of a variation of the molecular mechanics force field and the size of the QM region in QM/MM MD simulations on properties of the photoactive part of the blue light photoreceptor protein AppA. For comparison, we additionally performed classical MD simulations and studied the effect of a variation of the type of spatial boundary condition. The classical boundary conditions and the force field used in a QM/MM MD simulation are shown to have non-neglegible effects upon the structural and energetic properties of the protein which makes it advisable to minimize computational artifacts in QM/MM MD simulations by application of periodic boundary conditions and a thermodynamically calibrated force field. A comparison of the structural and energetic properties of MD simulations starting from two alternative, different X-ray structures for the blue light utilizing flavin protein in its dark state indicates a slight preference of the two force fields used for the so-called Anderson structure over the Jung structure.     

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.