Comparison of linear-scaling semiempirical methods and combined quantum mechanical/molecular mechanical methods for enzymic reactions. II. An energy decomposition analysis

Dual QM and MM Approach for Computing Equilibrium Isotope Fractionation Factor of Organic Species in Solution

A dual QM and MM approach for computing equilibrium isotope effects has been described. In the first partition, the potential energy surface is represented by a combined quantum mechanical and molecular mechanical (QM/MM) method, in which a solute molecule is treated quantum mechanically, and the remaining solvent molecules are approximated classically by molecular mechanics. In the second QM/MM partition, differential nuclear quantum effects responsible for the isotope effect are determined by a statistical mechanical double-averaging formalism, in which the nuclear centroid distribution is sampled classically by Newtonian molecular dynamics and the quantum mechanical spread of quantized particles about the centroid positions is treated using the path integral (PI) method. These partitions allow the potential energy surface to be properly represented such that the solute part is free of nuclear quantum effects for nuclear quantum mechanical simulations, and the double-averaging approach has the advantage of sampling efficiency for solvent configuration and for path integral convergence. Importantly, computational precision is achieved through free energy perturbation (FEP) theory to alchemically mutate one isotope into another. The PI-FEP approach is applied to model systems for the 18O enrichment found in cellulose of trees to determine the isotope enrichment factor of carbonyl compounds in water. The present method may be useful as a general tool for studying isotope fractionation in biological and geochemical systems.

Enzymatic Kinetic Isotope Effects from Path-Integral Free Energy Perturbation Theory

Path-integral free energy perturbation (PI-FEP) theory is presented to directly determine the ratio of quantum mechanical partition functions of different isotopologs in a single simulation. Furthermore, a double averaging strategy is used to carry out the practical simulation, separating the quantum mechanical path integral exactly into two separate calculations, one corresponding to a classical molecular dynamics simulation of the centroid coordinates, and another involving free-particle path-integral sampling over the classical, centroid positions. An integrated centroid path-integral free energy perturbation and umbrella sampling (PI-FEP/UM, or simply, PI-FEP) method along with bisection sampling was summarized, which provides an accurate and fast convergent method for computing kinetic isotope effects for chemical reactions in solution and in enzymes. The PI-FEP method is illustrated by a number of applications, to highlight the computational precision and accuracy, the rule of geometrical mean in kinetic isotope effects, enhanced nuclear quantum effects in enzyme catalysis, and protein dynamics on temperature dependence of kinetic isotope effects.

CAN SEMIEMPIRICAL QUANTUM MODELS CALCULATE THE BINDING ENERGY OF HYDROGEN BONDING FOR BIOLOGICAL SYSTEMS?

A modified semiempirical model named RM1BH, which is based on RM1 parameterizations, is proposed to simulate varied biological hydrogen-bonded systems. The RM1BH is formulated by adding Gaussian functions to the core–core repulsion items in original RM1 formula to reproduce the binding energies of hydrogen bonding of experimental and high-level computational results. In the parameterizations of our new model, 35 base-pair dimers, 18 amino acid residue dimers, 14 dimers between a base and an amino acid residue, and 20 other multimers were included. The results performed with RM1BH were compared with experimental values and the benchmark density-functional (B3LYP/6-31G**/BSSE) and Möller–Plesset perturbation (MP2/6-31G**/BSSE) calculations on various biological hydrogen-bonded systems. It was demonstrated that RM1BH model outperforms the PM3 and RM1 models in the calculations of the binding energies of biological hydrogen-bonded systems by very close agreement with the values of both high-level calculations and experiments. These results provide insight into the ideas, methods, and views of semiempirical modifications to investigate the weak interactions of biological systems.

A quantum mechanical approach to ligand binding — Calculation of ligand–protein binding affinities for stromelysin-1 (MMP-3) inhibitors

Linear-scaling quantum mechanical method was applied to calculate binding affinities of six stromelysin-1 (MMP-3) inhibitors with two different zinc binding groups (ZBGs). The entire protein and ligand–protein complexes were calculated using PM5 Hamiltonian, which enables the treatment of metal ion coordination, bond forming/breaking, and proton/charge transfers associated with the ligand binding process by the self-consistent field method. The calculated binding energies reproduce the binding-affinity trend observed experimentally.

Redefinition of rubisco carboxylase reaction reveals origin of water for hydration and new roles for active-site residues

Crystallographic, mutagenesis, kinetic, and computational studies on Rubisco over three decades have revealed much about its catalytic mechanism and the role played by several active-site residues. However, key questions remain unanswered. Specific details of the carboxylase and oxygenase mechanisms, required to underpin the rational re-engineering of Rubisco, are still speculative. Here we address critical gaps in knowledge with a definitive comprehensive computational investigation of the mechanism of carboxylase activity at the Rubisco active site. Density functional theory calculations (B3LYP/6-31G(d,p)) were performed on active-site fragment models of a size up to 77 atoms, not previously possible computationally. All amino acid residues suspected to play roles in the acid-base chemistry in the multistep reaction, and interacting directly with the central Mg (2+) atom and the reactive moiety of substrate and intermediates, were included. The results provide a firm basis for us to propose a novel mechanism for the entire sequence of reactions in the carboxylase catalysis and to define precise roles for the active-site residues, singly and in concert. In this mechanism, the carbamylated LYS201 plays a more limited role than previously proposed but is crucial for initiating the reaction by acting as a base in the enolization. We suggest a wider role for HIS294, with involvement in the carboxylation, hydration, and C2-C3 bond-scission steps, consistent with the suggestion of Harpel et al. (1998) but contrary to the consensus view of Cleland et al. (1998). In contrast to the common assumption that the water molecule for the hydration step comes from within the active site, we propose that the Mg-coordinated water is not dissociated at the start of the gas-addition reaction but rather remains coordinated and is used for the hydration of the C3 carbon atom. New roles are also proposed for LYS175, GLU204, and HIS294. The mechanism suggests roles in the gas-addition step for residues in three spatially distinct regions of the active site, HIS294 and LYS334 in the C-terminal domain of the large subunit (LSU), but also hitherto unsuspected roles for a cluster of three residues (ASN123, GLU60, and TYR20) in the N-terminal domain of the partner LSU of the dimer containing the active site. Our new mechanism is supported by existing experimental data, provides new convincing interpretations of previously puzzling data, and allows new insights into mutational strategies for improving Rubisco activity.

Self-consistent polarization neglect of diatomic differential overlap: application to water clusters

Semiempirical self-consistent field (SCF) methods based on the neglect of diatomic differential overlap (NDDO) formalism have the ability to treat the formation and breaking of chemical bonds but have been found to poorly describe hydrogen bonding and weak electrostatic complexes. In contrast, most empirical potentials are not able to describe bond breaking and formation but have the ability to add missing elements of hydrogen bonding by using classical electrostatic interactions. We present a new method which combines aspects of both NDDO-based SCF techniques and classical descriptions of polarization to describe the diffuse nature of the electronic wavefunction in a self-consistent manner. We develop the "self-consistent polarization neglect of diatomic differential overlap" (SCP-NDDO) theory with the additional description of molecular dispersion developed as a second-order perturbation theory expression. The current study seeks to model water-water interactions as a test case. To this end, we have parametrized the method to accurate ab initio complete basis set limit estimates of small water cluster binding energies of Xantheas and co-workers [J. Chem. Phys. 116, 1493 (2002); 120, 823 (2004)]. Overall agreement with the ab initio binding energies (n=2-6, and 8) is achieved with a rms error of 0.19 kcal/mol. We achieve noticeable improvements in the structure, vibrational frequencies, and energetic predictions of water clusters (n< or =21) relative to standard NDDO-based methods.

Hybrid quantum and classical methods for computing kinetic isotope effects of chemical reactions in solutions and in enzymes

A method for incorporating quantum mechanics into enzyme kinetics modeling is presented. Three aspects are emphasized: 1) combined quantum mechanical and molecular mechanical methods are used to represent the potential energy surface for modeling bond forming and breaking processes, 2) instantaneous normal mode analyses are used to incorporate quantum vibrational free energies to the classical potential of mean force, and 3) multidimensional tunneling methods are used to estimate quantum effects on the reaction coordinate motion. Centroid path integral simulations are described to make quantum corrections to the classical potential of mean force. In this method, the nuclear quantum vibrational and tunneling contributions are not separable. An integrated centroid path integral-free energy perturbation and umbrella sampling (PI-FEP/UM) method along with a bisection sampling procedure was summarized, which provides an accurate, easily convergent method for computing kinetic isotope effects for chemical reactions in solution and in enzymes. In the ensemble-averaged variational transition state theory with multidimensional tunneling (EA-VTST/MT), these three aspects of quantum mechanical effects can be individually treated, providing useful insights into the mechanism of enzymatic reactions. These methods are illustrated by applications to a model process in the gas phase, the decarboxylation reaction of N-methyl picolinate in water, and the proton abstraction and reprotonation process catalyzed by alanine racemase. These examples show that the incorporation of quantum mechanical effects is essential for enzyme kinetics simulations.

Extending the Power of Quantum Chemistry to Large Systems with the Fragment Molecular Orbital Method

Following the brief review of the modern fragment-based methods and other approaches to perform quantum-mechanical calculations of large systems, the theoretical development of the fragment molecular orbital method (FMO) is covered in detail, with the emphasis on the physical properties, which can be computed with FMO. The FMO-based polarizable continuum model (PCM) for treating the solvent effects in large systems and the pair interaction energy decomposition analysis (PIEDA) are described in some detail, and a range of applications of FMO to biological studies is introduced. The factors determining the relative stability of polypeptide conformers (alpha-helix, beta-turn, and extended form) are elucidated using FMO/PCM and PIEDA, and the interactions in the Trp-cage miniprotein construct (PDB: 1L2Y) are analyzed using PIEDA.

Redistributed charge and dipole schemes for combined quantum mechanical and molecular mechanical calculations

Special care is needed in carrying out combined quantum mechanical and molecular mechanical (QM/MM) calculations if the QM/MM boundary passes through a covalent bond. The present paper discusses the importance of correctly handling the MM partial point charges at the QM/MM boundary, and in particular, it contributes in two aspects: (1) Two schemes, namely, the redistributed charge (RC) scheme and the redistributed charge and dipole (RCD) scheme, are introduced to handle link atoms in QM/MM calculations. In both schemes, the point charge at the MM boundary atom that is replaced by the link atom is redistributed to the midpoint of the bonds that connect the MM boundary atom and its neighboring MM atoms. These redistributed charges serve as classical mimics for the auxiliary orbitals associated with the MM host atom in the generalized hybrid orbital (GHO) method. In the RCD scheme, the dipoles of these bonds are preserved by further adjustment of the values of the redistributed charges. The treatments are justified as classical analogues of the QM description given by the GHO method. (2) The new methods are compared quantitatively to similar methods that were suggested by previous work, namely, a shifted-charge scheme and three eliminated-charge schemes. The comparisons were carried out for a series of molecules in terms of proton affinities and geometries. Point charges derived from various charge models were tested. The results demonstrate that it is critical to preserve charge and bond dipole and that it is important to use accurate MM point charges in QM/MM boundary treatments. The RCD scheme was further applied to study the H atom transfer reaction CH3 + CH3CH2CH2OH --> CH4 + CH2CH2CH2OH. Various QM levels of theory were tested to demonstrate the generality of the methodology. It is encouraging to find that the QM/MM calculations obtained a reaction energy, barrier height, saddle-point geometry, and imaginary frequency at the saddle point in quite good agreement with full QM calculations at the same level. Furthermore, analysis based on energy decomposition revealed the quantitatively similar interaction energies between the QM and the MM subsystems for the reactant, for the saddle point, and for the product. These interaction energies almost cancel each other energetically, resulting in negligibly small net effects on the reaction energy and barrier height. However, the charge distribution of the QM atoms is greatly affected by the polarization effect of the MM point charges. The QM/MM charge distribution agrees much better with full QM results than does the unpolarized charge distribution of the capped primary subsystem.

RM1: A reparameterization of AM1 for H, C, N, O, P, S, F, Cl, Br, and I

Twenty years ago, the landmark AM1 was introduced, and has since had an increasingly wide following among chemists due to its consistently good results and time-tested reliability--being presently available in countless computational quantum chemistry programs. However, semiempirical molecular orbital models still are of limited accuracy and need to be improved if the full potential of new linear scaling techniques, such as MOZYME and LocalSCF, is to be realized. Accordingly, in this article we present RM1 (Recife Model 1): a reparameterization of AM1. As before, the properties used in the parameterization procedure were: heats of formation, dipole moments, ionization potentials and geometric variables (bond lengths and angles). Considering that the vast majority of molecules of importance to life can be assembled by using only six elements: C, H, N, O, P, and S, and that by adding the halogens we can now build most molecules of importance to pharmaceutical research, our training set consisted of 1736 molecules, representative of organic and biochemistry, containing C, H, N, O, P, S, F, Cl, Br, and I atoms. Unlike AM1, and similar to PM3, all RM1 parameters have been optimized. For enthalpies of formation, dipole moments, ionization potentials, and interatomic distances, the average errors in RM1, for the 1736 molecules, are less than those for AM1, PM3, and PM5. Indeed, the average errors in kcal x mol(-1) of the enthalpies of formation for AM1, PM3, and PM5 are 11.15, 7.98, and 6.03, whereas for RM1 this value is 5.77. The errors, in Debye, of the dipole moments for AM1, PM3, PM5, and RM1 are, respectively, 0.37, 0.38, 0.50, and 0.34. Likewise, the respective errors for the ionization potentials, in eV, are 0.60, 0.55, 0.48, and 0.45, and the respective errors, in angstroms, for the interatomic distances are 0.036, 0.029, 0.037, and 0.027. The RM1 average error in bond angles of 6.82 degrees is only slightly higher than the AM1 figure of 5.88 degrees, and both are much smaller than the PM3 and PM5 figures of 6.98 degrees and 9.83 degrees, respectively. Moreover, a known error in PM3 nitrogen charges is corrected in RM1. Therefore, RM1 represents an improvement over AM1 and its similar successor PM3, and is probably very competitive with PM5, which is a somewhat different model, and not fully disclosed. RM1 possesses the same analytical construct and the same number of parameters for each atom as AM1, and, therefore, can be easily implemented in any software that already has AM1, not requiring any change in any line of code, with the sole exception of the values of the parameters themselves.