Comment on “A stationary‐wave model of enzyme catalysis” by Carlo Canepa

Energy barriers for enzyme‐catalyzed reactions calculated with quantum mechanics/molecular mechanics (QM/MM) and empirical valence bond (EVB) methods can be in excellent agreement with activation energies derived from experiments, supporting the applicability of transition state theory for enzymic reactions. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011     

Implementation and performance of the artificial force induced reaction method in the GRRM17 program

This article reports implementation and performance of the artificial force induced reaction (AFIR) method in the upcoming 2017 version of GRRM program (GRRM17). The AFIR method, which is one of automated reaction path search methods, induces geometrical deformations in a system by pushing or pulling fragments defined in the system by an artificial force. In GRRM17, three different algorithms, that is, multicomponent algorithm (MC‐AFIR), single‐component algorithm (SC‐AFIR), and double‐sphere algorithm (DS‐AFIR), are available, where the MC‐AFIR was the only algorithm which has been available in the previous 2014 version. The MC‐AFIR does automated sampling of reaction pathways between two or more reactant molecules. The SC‐AFIR performs automated generation of global or semiglobal reaction path network. The DS‐AFIR finds a single path between given two structures. Exploration of minimum energy structures within the hypersurface in which two different electronic states degenerate, and an interface with the quantum mechanics/molecular mechanics method, are also described. A code termed SAFIRE will also be available, as a visualization software for complicated reaction path networks. © 2017 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Producing DFT/MM enzyme reaction trajectories from SCC‐DFTB/MM driving forces to probe the underlying electronics of a glycosyltransferase reaction

The SCC‐DFTB/MIO/CHARMM free energy surface for a glycosyltransferase, TcTS, is benchmarked against a DFT/MM reaction trajectory using the same CHARMM MM force field ported to the NWChem package. The popular B3LYP functional, against which the MIO parameter set was parameterized is used to optimize TS structures and run DFT reaction dynamics. A novel approach was used to generate reaction forces from a SCC‐DFTB/MIO/CHARMM reaction surface to drive B3LYP/6‐31G/MM and B3LYP/6‐31G(d)/MM reaction trajectories. Although TS structures compare favorably, differences stemming primarily from a minimal basis set approximation prevented a successful 6‐31G(d) FEARCF reaction dynamics trajectory. None the less, the dynamic evolution of the B3LYP/6‐31G/MM‐computed electron density provided an opportunity to perform NBO analysis along the reaction trajectory. Here, we illustrate that a successful ab initio reaction trajectory is computationally accessible when the underlying potential energy function of the semi‐empirical method used to produce driving forces is sufficiently close to the ab initio potential. © 2017 Wiley Periodicals, Inc.     

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.     

Ab initio QM/MM study of excited state electron transfer between pyrene and 4,4′‐bis(dimethylamino)‐diphenylmethane with different solvent systems: Role of hydrogen bonding within solvent molecules

The exciplex is a charge transfer species formed in the process of electron transfer between an electron donor and an electron acceptor and hence is very sensitive to solvent polarity. In order to understand the role of solvent in exciplex formation between pyrene (PY) and 4,4′‐bis(dimethylamino)diphenylmethane (DMDPM), we used two types of solvent approximations: an implicit solvent model and an explicit solvent model. The difference in energies between the excited and the meta‐stable Frank–Condon state (ΔE) of the structures were assumed to correspond to the emission maximum of the exciplex in different solvents. The ΔE values show the trend of stabilization of the exciplex with an increase in solvent polarity. This trend in stabilization is substantially more prominent in the explicit solvent model than that with the implicit solvent model. The ΔE value obtained in methanol reflects equal stabilization compared to that in a more polar solvent, N,N‐dimethylformamide. This extra stabilization of the exciplex may be explained on the basis of the H‐bonding capability of the protic solvent, methanol. © 2004 Wiley Periodicals, Inc. Int J Quantum Chem, 2005     

Implementation of an adaptive umbrella sampling method for the calculation of multidimensional potential of mean force of chemical reactions in solution

We describe the implementation of an adaptive umbrella sampling method, making use of the weighted histogram analysis method, for computing multidimensional potential of mean force for chemical reaction in solution. The approach is illustrated by investigating the effect of aqueous solution on the free energy surface for the proton transfer reaction of [H(3)N-H-NH(3)](+) using a combined quantum mechanical and molecular mechanical AM1/TIP3P potential.     

Theoretical investigation of enantioselectivity of cage‐like supramolecular assembly: The insights into the shape complementarity and host–guest interaction

Enantioselectivity in the aza‐Cope rearrangement of a guest molecule encapsulated in a cage‐like supramolecular assembly [Ga₄L₆]^12? [L = 1,5‐bis(2',3'‐dihydroxybenzamido)naphthalene] is investigated using density functional theory and ab initio molecular orbital calculations. Reaction pathways leading to R‐ and S‐enantiomers encapsulated in the [Ga₄L₆]^12? are explored. The reaction barriers and the stabilities of the prochiral structures differed in the [Ga₄L₆]^12?, resulting that the product with an R structure is favorably produced in the Δ‐structure [Ga₄L₆]^12?. The large energy difference in the prochiral structures in the [Ga₄L₆]^12? was attributed to the deformation of the bulky substituent. The host–guest interaction energy raises the reaction barrier for the product with an S structure. The previous study suggested that the different stability of the prochiral substrates in the assembly was the origin of the enantioselectivity, and the suggestion is supported by our computational finding. In addition, our results show that the difference in the reaction barriers also importantly contributes to the enantioselectivity. © 2015 Wiley Periodicals, Inc.     

Increasing the time step with mass scaling in Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations

Born‐Oppenheimer ab initio QM/MM molecular dynamics simulation with umbrella sampling is a state‐of‐the‐art approach to calculate free energy profiles of chemical reactions in complex systems. To further improve its computational efficiency, a mass‐scaling method with the increased time step in MD simulations has been explored and tested. It is found that by increasing the hydrogen mass to 10 amu, a time step of 3 fs can be employed in ab initio QM/MM MD simulations. In all our three test cases, including two solution reactions and one enzyme reaction, the resulted reaction free energy profiles with 3 fs time step and mass scaling are found to be in excellent agreement with the corresponding simulation results using 1 fs time step and the normal mass. These results indicate that for Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling, the mass‐scaling method can significantly reduce its computational cost while has little effect on the calculated free energy profiles. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Application of QM simulations and multivariate analysis in the study of alkene reactivity in the zeolite H‐ZSM5

Reported herein are the results of an investigation into the effect of the extended framework of the zeolite ZSM‐5 on the reaction energetics and structures of (a) the physisorbed complex formed between the zeolite and six alkenes, (b) the corresponding chemisorbed alkoxide intermediate and (c) the transition states (TS) connecting the two. For this, quantum mechanical (QM) simulations of ZSM‐5 in the presence and absence of the zeolite framework have been employed. A 46T density functional theory (DFT) cluster model and a 3T:46T DFT:UFF ONIOM model are used to represent the former scenario and a simple 3T DFT cluster model for the latter. The structural implications of neglecting the zeolite framework have been rigorously compared using the multivariate statistical method principal components analysis (PCA). This method allows one to assess the correlated nature of the changes in structure along the reaction coordinate, for multiple different alkenes, in a facile, reliable way. Copyright © 2008 John Wiley & Sons, Ltd.     

Mechanism of chemical reactions in the active site of aspartate N-acetyltransferase NAT8L revealed by molecular modeling

The results of a computational study of the synthesis of a key brain metabolite, N-acetyl-l-aspartate, catalyzed by aspartate N-acetyltransferase, encoded by the NAT8L gene, are reported. The reaction Gibbs energy profiles were computed using molecular dynamics simulations with interaction potentials estimated on-the-fly by the quantum mechanics/molecular mechanics QM(PBE0/6-31G**)/MM(CHARMM) approach. The revealed reaction mechanism includes four elementary steps with corresponding activation energies not exceeding 14 kcal mol^?1     

Minimum MD simulation length required to achieve reliable results in free energy perturbation calculations: Case study of relative binding free energies of fructose‐1,6‐bisphosphatase inhibitors

In an attempt to establish the criteria for the length of simulation to achieve the desired convergence of free energy calculations, two studies were carried out on chosen complexes of FBPase‐AMP mimics. Calculations were performed for varied length of simulations and for different starting configurations using both conventional‐ and QM/MM‐FEP methods. The results demonstrate that for small perturbations, 1248 ps simulation time could be regarded a reasonable yardstick to achieve convergence of the results. As the simulation time is extended, the errors associated with free energy calculations also gradually tapers off. Moreover, when starting the simulation from different initial configurations of the systems, the results are not changed significantly, when performed for 1248 ps. This study carried on FBPase‐AMP mimics corroborates well with our previous successful demonstration of requirement of simulation time for solvation studies, both by conventional and ab initio FEP. The establishment of aforementioned criteria of simulation length serves a useful benchmark in drug design efforts using FEP methodologies, to draw a meaningful and unequivocal conclusion. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011     

Reaction mechanism of NO with hydrolysates of NAMI-A: an MD simulation by combining the QM/MM(ABEEM) with the MD-FEP method

Nitrosylation reaction mechanisms of the hydrolysates of NAMI-A and hydrolysis reactions of ruthenium nitrosyl complexes were investigated in the triplet state and the singlet state. Activation free energies were calculated by combining the QM/MM(ABEEM) method with free energy perturbation theory, and the explicit solvent environment was simulated by an ABEEMσπ polarizable force field. Our results demonstrate that nitrosylation reactions of the hydrolysates of NAMI-A occur in both the triplet and the singlet states. The Ru-N-O angle of the triplet ruthenium nitrosyl complexes is in the range of 132.0°–138.2°. However, all the ruthenium nitrosyl complexes at the singlet state show an almost linear Ru-N-O angle. The nitrosylation reaction happens prior to the hydrolysis reaction for the first-step hydrolysates. The activation free energies of the nitrosylation reactions show that the H₂O-NO exchange reaction of [RuCl₄(Im)(H₂O)] in the singlet spin sate is the most likely one. Comparing with the activation free energies of the hydrolysis reactions of the ruthenium nitrosyl complexes, the results indicate that the rate of the DMSO–H₂O exchange reaction of [RuCl₃(NO)(Im)(DMSO)] is faster than that of [RuCl₃(H₂O)(Im)(DMSO)] in both the triplet spin state and the singlet spin state. © 2018 Wiley Periodicals, Inc.     

Structure,strain, and reorganization energy of blue copper models in the protein

The copper coordination geometry in the blue copper proteins plastocyanin, nitrite reductase, cucumber basic protein, and azurin has been studied by combined density functional (B3LYP) and molecular mechanical methods. Compared to quantum chemical vacuum calculations, a significant improvement of the geometry is seen (toward the experimental structures) not only for the dihedral angles of the ligands but also for the bond lengths and angles around the copper ion. The flexible Cu–S_Met bond is well reproduced in the oxidized structures, whereas it is too long in some of the reduced complexes (too short in vacuum). The change in the geometry compared to the vacuum state costs 33–66 kJ/mol. If the covalent bonds between the ligands and the protein are broken, this energy decreases by ∼25 kJ/mol, which is an estimate of the covalent strain. This is similar to what is found for other proteins, so the blue copper proteins are not more strained than other metalloproteins. The inner‐sphere self‐exchange reorganization energy of all four proteins are ∼30 kJ/mol. This is 30–50 kJ/mol lower than in vacuum. The decrease is caused by dielectric and electrostatic effects in the protein, especially the hydrogen bond(s) to the cysteine copper ligands and not by covalent strain. © 2001 John Wiley & Sons, Inc. Int J Quant Chem 81: 335–347, 2001     

Intra‐ and intermolecular interactions in crystals of polar molecules. A study by the mixed quantum mechanical/molecular mechanical SCMP‐NDDO method

Stabilization energies of crystals of polar molecules were calculated with the recently developed NDDO‐SCMP method that determines the wave function of a subunit embedded in the symmetrical environment constituted by the copies of the subunit. The total stabilization energies were decomposed into four components. The deformation energy is the difference between the energy of the molecule in the geometries adopted in the crystal on the one hand, and in vacuo, on the other hand. Further energy components are derived from the molecular geometry found in the crystal phase. The electrostatic component is the interaction energy of the molecule with the crystal field, corresponding to the charge distribution obtained in vacuo. The polarization component is the energy lowering resulted in the self‐consistent optimization of the wave function in the crystal field. The rest of the stabilization energy is attributed to the dispersion–repulsion component, and is calculated from an empirical potential function. The major novelty of this decomposition scheme is the introduction of the deformation energy. It requires the optimization of the structural parameters, including the molecular geometry, the intermolecular coordinates, and the cell parameters of the crystal. The optimization is performed using the recently implemented forces in the SCMP‐NDDO method, and this new feature is discussed in detail. The calculation of the deformation energy is particularly important to obtain stabilization energies for crystals in which the molecular geometry differs considerably from that corresponding to the energy minimum of the isolated molecule. As an example, crystals of diastereoisomeric salts are investigated. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1679–1690, 2001     

A quantum chemical study on the mechanism of the consecutive addition of HCN to propionitrile

MINDO/3 MO method has been used to study the mechanism of the consecutive addition of HCN to propionitrile. The results obtained for the first five steps show that the reaction is exothermic, and step 1 is the rate determining step.     

Branch migration of Holliday junction in RuvA tetramer complex studied by umbrella sampling simulation using a path‐search algorithm

Branch migration of the Holliday junction takes place at the center of the RuvA tetramer. To elucidate how branch migration occurs, umbrella sampling simulations were performed for complexes of the RuvA tetramer and Holliday junction DNA. Although conventional umbrella sampling simulations set sampling points a priori, the umbrella sampling simulation in this study set the sampling points one by one in order to search for a realistic path of the branch migration during the simulations. Starting from the X‐ray structure of the complex, in which the hydrogen bonds between two base‐pairs were unformed, the hydrogen bonds between the next base‐pairs of the shrinking stems were observed to start to disconnect. At the intermediate stage, three or four of the eight unpaired bases interacted closely with the acidic pins from RuvA. During the final stage, these bases moved away from the pins and formed the hydrogen bonds of the new base‐pairs of the growing stems. The free‐energy profile along this reaction path showed that the intermediate stage was a meta‐stable state between two free‐energy barriers of about 10 to 15 kcal/mol. These results imply that the pins play an important role in stabilizing the interactions between the pins and the unpaired base‐pairs. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010