Enzyme mechanisms from molecular modeling and isotope effects

The application of kinetic isotope effects and molecular modeling to characterize three enzyme-catalyzed reactions is presented; the mechanism of the chloroacid dehalogenase catalyzed reaction is approached using chlorine kinetic isotope effects and solvent kinetic isotope effects. The pre-steady-state phase of the reaction catalyzed by methylmalonyl-CoA mutase is approached by different QM/MM schemes and the results are validated by comparison with the experimental value of the deuterium kinetic isotope effect. Finally, a procedure for improving QM/MM calculations is illustrated by analysis of the trihydroxynaphthalene reductase-catalyzed reaction.     

Use of multistate Bennett acceptance ratio method for free-energy calculations from enhanced sampling and free-energy perturbation

Multistate Bennett acceptance ratio (MBAR) works as a method to analyze molecular dynamics (MD) simulation data after the simulations have been finished. It is widely used to estimate free-energy changes between different states and averaged properties at the states of interest. MBAR allows us to treat a wide range of states from those at different temperature/pressure to those with different model parameters. Due to the broad applicability, the MBAR equations are rather difficult to apply for free-energy calculations using different types of MD simulations including enhanced conformational sampling methods and free-energy perturbation. In this review, we first summarize the basic theory of the MBAR equations and categorize the representative usages into the following four: (i) perturbation, (ii) scaling, (iii) accumulation, and (iv) full potential energy. For each, we explain how to prepare input data using MD simulation trajectories for solving the MBAR equations. MBAR is also useful to estimate reliable free-energy differences using MD trajectories based on a semi-empirical quantum mechanics/molecular mechanics (QM/MM) model and ab initio QM/MM energy calculations on the MD snapshots. We also explain how to use the MBAR software in the GENESIS package, which we call mbar_analysis, for the four representative cases. The proposed estimations of free-energy changes and thermodynamic averages are effective and useful for various biomolecular systems.     

Proton transfer in carbonic anhydrase is controlled by electrostatics rather than the orientation of the acceptor

Combined quantum mechanical/molecular mechanical (QM/MM) simulations are carried out to analyze factors that dictate the proton transfer in carbonic anhydrase II (CAII), an enzyme that has been used as a prototypical example of long-range proton transfers in biomolecules. In contrast to the long-held conjecture in the experimental literature, the computed potentials of mean force (PMF) suggest that the proton transfer in CAII is not very sensitive to the orientation of the acceptor group (His 64) and, therefore, the number of water molecules that bridge the donor (zinc-water) and acceptor groups. Perturbative analysis indicates that a series of polar and charged residues close to the transfer pathways make the dominant contribution to the barrier and exothermicity of the proton transfer reaction, thus supporting the proposal from previous studies of Warshel and co-workers using a somewhat simpler QM/MM model that electrostatic interactions play a major role in the proton transfer in CAII. The PMF results are in striking contrast to previous analysis using the same QM/MM method but an ensemble of minimum energy path (MEP) calculations, which found a steep dependence of the barrier height on the number of bridging water molecules. Analysis of the configurations sampled in the PMF and MEP simulations suggests that this difference arises because the PMF simulations sample a largely stepwise mechanism while the local MEP calculations artificially favored concerted transfers due to the specific protocol used to generate the initial configurations. Therefore, this study presents a compelling argument for carrying out proper conformational sampling in the study of long-range proton transfers. Finally, we illustrate that Phi analysis, which has been widely used in protein folding studies, can potentially generate new mechanistic information for long-range proton transfers regarding the sequence of events. The results of the perturbation analysis and the Phi analysis provide opportunities for experimentally testing the mechanistic proposals from this study and our recent work in which a stepwise "proton hole" transfer pathway has been proposed.     

Computational and experimental studies on the catalytic mechanism of biliverdin-IXbeta reductase

BVR-B (biliverdin-IXbeta reductase) also known as FR (flavin reductase) is a promiscuous enzyme catalysing the pyridine-nucleotide-dependent reduction of a variety of flavins, biliverdins, PQQ (pyrroloquinoline quinone) and ferric ion. Mechanistically it is a good model for BVR-A (biliverdin-IXalpha reductase), a potential pharmacological target for neonatal jaundice and also a potential target for adjunct therapy to maintain protective levels of biliverdin-IXalpha during organ transplantation. In a commentary on the structure of BVR-B it was noted that one outstanding issue remained: whether the mechanism was a concerted hydride transfer followed by protonation of a pyrrolic anion or protonation of the pyrrole followed by hydride transfer. In the present study we have attempted to address this question using QM/MM (quantum mechanics/molecular mechanics) calculations. QM/MM potential energy surfaces show that the lowest energy pathway proceeds with a positively charged pyrrole intermediate via two transition states. These initial calculations were performed with His(153) as the source of the proton. However site-directed mutagenesis studies with both the H153A and the H153N mutant reveal that His(153) is not required for catalytic activity. We have repeated the calculation with a solvent hydroxonium donor and obtain a similar energy landscape indicating that protonation of the pyrrole is the most likely first step followed by hydride transfer and that the required proton may come from bulk solvent. The implications of the present study for the design of inhibitors of BVR-A are discussed.     

Lys169 of Human Glucokinase Is a Determinant for Glucose Phosphorylation: Implication for the Atomic Mechanism of Glucokinase Catalysis

Glucokinase (GK), a glucose sensor, maintains plasma glucose homeostasis via phosphorylation of glucose and is a potential therapeutic target for treating maturity-onset diabetes of the young (MODY) and persistent hyperinsulinemic hypoglycemia of infancy (PHHI). To characterize the catalytic mechanism of glucose phosphorylation by GK, we combined molecular modeling, molecular dynamics (MD) simulations, quantum mechanics/molecular mechanics (QM/MM) calculations, experimental mutagenesis and enzymatic kinetic analysis on both wild-type and mutated GK. Our three-dimensional (3D) model of the GK-Mg²⁺-ATP-glucose (GMAG) complex, is in agreement with a large number of mutagenesis data, and elucidates atomic information of the catalytic site in GK for glucose phosphorylation. A 10-ns MD simulation of the GMAG complex revealed that Lys169 plays a dominant role in glucose phosphorylation. This prediction was verified by experimental mutagenesis of GK (K169A) and enzymatic kinetic analyses of glucose phosphorylation. QM/MM calculations were further used to study the role of Lys169 in the catalytic mechanism of the glucose phosphorylation and we found that Lys169 enhances the binding of GK with both ATP and glucose by serving as a bridge between ATP and glucose. More importantly, Lys169 directly participates in the glucose phosphorylation as a general acid catalyst. Our findings provide mechanistic details of glucose phorphorylation catalyzed by GK, and are important for understanding the pathogenic mechanism of MODY.     

Multiscale methods for macromolecular simulations

In this article we review the key modeling tools available for simulating biomolecular systems. We consider recent developments and representative applications of mixed quantum mechanics/molecular mechanics (QM/MM), elastic network models (ENMs), coarse-grained molecular dynamics, and grid-based tools for calculating interactions between essentially rigid protein assemblies. We consider how the different length scales can be coupled, both in a sequential fashion (e.g. a coarse-grained or grid model using parameterization from MD simulations), and via concurrent approaches, where the calculations are performed together and together control the progression of the simulation. We suggest how the concurrent coupling approach familiar in the context of QM/MM calculations can be generalized, and describe how this has been done in the CHARMM macromolecular simulation package.     

Reaction-path energetics and kinetics of the hydride transfer reaction catalyzed by dihydrofolate reductase

We have studied the hydride transfer reaction catalyzed by the enzyme dihydrofolate reductase (DHFR) and the coenzyme nicotinamide adenine dinucleotide phosphate (NADPH); the substrate is 5-protonated 7,8-dihydrofolate, and the product is tetrahydrofolate. The potential energy surface is modeled by a combined quantum mechanical-molecular mechanical (QM/MM) method employing Austin model 1 (AM1) and a simple valence bond potential for 69 QM atoms and employing the CHARMM22 and TIP3P molecular mechanics force fields for the other 21 399 atoms; the QM and MM regions are joined by two boundary atoms treated by the generalized hybrid orbital (GHO) method. All simulations are carried out using periodic boundary conditions at neutral pH and 298 K. In stage 1, a reaction coordinate is defined as the difference between the breaking and forming bond distances to the hydride ion, and a quasithermodynamic free energy profile is calculated along this reaction coordinate. This calculation includes quantization effects on bound vibrations but not on the reaction coordinate, and it is used to locate the variational transition state that defines a transition state ensemble. Then, the key interactions at the reactant, variational transition state, and product are analyzed in terms of both bond distances and electrostatic energies. The results of both analyses support the conclusion derived from previous mutational studies that the M20 loop of DHFR makes an important contribution to the electrostatic stabilization of the hydride transfer transition state. Third, transmission coefficients (including recrossing factors and multidimensional tunneling) are calculated and averaged over the transition state ensemble. These averaged transmission coefficients, combined with the quasithermodynamic free energy profile determined in stage 1, allow us to calculate rate constants, phenomenological free energies of activation, and primary and secondary kinetic isotope effects. A primary kinetic isotope effect (KIE) of 2.8 has been obtained, in good agreement with the experimentally determined value of 3.0 and with the value 3.2 calculated previously. The primary KIE is mainly a consequence of the quantization of bound vibrations. In contrast, the secondary KIE, with a value of 1.13, is almost entirely due to dynamical effects on the reaction coordinate, especially tunneling.     

Binding free energy calculation with QM/MM hybrid methods for Abl-Kinase inhibitor

We report a Quantum mechanics/Molecular Mechanics–Poisson-Boltzmann/ Surface Area (QM/MM-PB/SA) method to calculate the binding free energy of c-Abl human tyrosine kinase by combining the QM and MM principles where the ligand is treated quantum mechanically and the rest of the receptor by classical molecular mechanics. To study the role of entropy and the flexibility of the protein ligand complex in a solvated environment, molecular dynamics calculations are performed using a hybrid QM/MM approach. This work shows that the results of the QM/MM approach are strongly correlated with the binding affinity. The QM/MM interaction energy in our reported study confirms the importance of electronic and polarization contributions, which are often neglected in classical MM-PB/SA calculations. Moreover, a comparison of semi-empirical methods like DFTB-SCC, PM3, MNDO, MNDO-PDDG, and PDDG-PM3 is also performed. The results of the study show that the implementation of a DFTB-SCC semi-empirical Hamiltonian that is derived from DFT gives better results than other methods. We have performed such studies using the AMBER molecular dynamic package for the first time. The calculated binding free energy is also in agreement with the experimentally determined binding affinity for c-Abl tyrosine kinase complex with Imatinib.     

Simulated <Superscript>18</Superscript>O kinetic isotope effects in enzymatic hydrolysis of guanosine triphosphate

We compare the computed on the base of quantum mechanical-molecular mechanical (QM/MM) modeling kinetic isotope effects (KIEs) for a series of the ¹⁸O-labeled substrates in enzymatic hydrolysis of guanosine triphosphate (GTP) with those measured experimentally. Following the quantitative structure-activity relationship concept, we introduce the correlation between KIEs and structure of substrates with the help of a labeling index, which also aids better imaging of presentation of both experimental and theoretical data. An evident correlation of the computed and measured KIEs provides support to the predominantly dissociative-type reaction mechanism of enzymatic GTP hydrolysis predicted in QM/MM simulations.     

Determinants of reactivity and selectivity in soluble epoxide hydrolase from quantum mechanics/molecular mechanics modeling

Soluble epoxide hydrolase (sEH) is an enzyme involved in drug metabolism that catalyzes the hydrolysis of epoxides to form their corresponding diols. sEH has a broad substrate range and shows high regio- and enantioselectivity for nucleophilic ring opening by Asp333. Epoxide hydrolases therefore have potential synthetic applications. We have used combined quantum mechanics/molecular mechanics (QM/MM) umbrella sampling molecular dynamics (MD) simulations (at the AM1/CHARMM22 level) and high-level ab initio (SCS-MP2) QM/MM calculations to analyze the reactions, and determinants of selectivity, for two substrates: trans-stilbene oxide (t-SO) and trans-diphenylpropene oxide (t-DPPO). The calculated free energy barriers from the QM/MM (AM1/CHARMM22) umbrella sampling MD simulations show a lower barrier for phenyl attack in t-DPPO, compared with that for benzylic attack, in agreement with experiment. Activation barriers in agreement with experimental rate constants are obtained only with the highest level of QM theory (SCS-MP2) used. Our results show that the selectivity of the ring-opening reaction is influenced by several factors, including proximity to the nucleophile, electronic stabilization of the transition state, and hydrogen bonding to two active site tyrosine residues. The protonation state of His523 during nucleophilic attack has also been investigated, and our results show that the protonated form is most consistent with experimental findings. The work presented here illustrates how determinants of selectivity can be identified from QM/MM simulations. These insights may also provide useful information for the design of novel catalysts for use in the synthesis of enantiopure compounds.     

Catalytic mechanism and energy barriers for butyrylcholinesterase-catalyzed hydrolysis of cocaine

The geometries of the transition states, intermediates, and prereactive enzyme-substrate complex and the corresponding energy barriers have been determined by performing hybrid quantum mechanical/molecular mechanical (QM/MM) calculations on butyrylcholinesterase (BChE)-catalyzed hydrolysis of (-)- and (+)-cocaine. The energy barriers were evaluated by performing QM/MM calculations with the QM method at the MP2/6-31+G* level and the MM method using the AMBER force field. These calculations allow us to account for the protein environmental effects on the transition states and energy barriers of these enzymatic reactions, showing remarkable effects of the protein environment on intermolecular hydrogen bonding (with an oxyanion hole), which is crucial for the transition state stabilization and, therefore, on the energy barriers. The calculated energy barriers are consistent with available experimental kinetic data. The highest barrier calculated for BChE-catalyzed hydrolysis of (-)- and (+)-cocaine is associated with the third reaction step, but the energy barrier calculated for the first step is close to the highest and is so sensitive to the protein environment that the first reaction step can be rate determining for (-)-cocaine hydrolysis catalyzed by a BChE mutant. The computational results provide valuable insights into future design of BChE mutants with a higher catalytic activity for (-)-cocaine.     

How does the protein environment optimize the thermodynamics of thiol sulfenylation? Insights from model systems to QM/MM calculations on human 2-Cys peroxiredoxin

Protein thiol/sulfenic acid oxidation potentials provide a tool to select specific oxidation agents, but are experimentally difficult to obtain. Here, insights into the thiol sulfenylation thermodynamics are obtained from model calculations on small systems and from a quantum mechanics/molecular mechanics (QM/MM) analysis on human 2-Cys peroxiredoxin thioredoxin peroxidase B (Tpx-B). To study thiol sulfenylation in Tpx-B, our recently developed computational method to determine reduction potentials relatively compared to a reference system and based on reaction energies reduction potential from electronic energies is updated. Tpx-B forms a sulfenic acid (R-SO^?) on one of its active site cysteines during reactive oxygen scavenging. The observed effect of the conserved active site residues is consistent with the observed hydrogen bond interactions in the QM/MM optimized Tpx-B structures and with free energy calculations on small model systems. The ligand effect could be linked to the complexation energies of ligand L with CH₃S^? and CH₃SO^?. Compared to QM only calculations on Tpx-B’s active site, the QM/MM calculations give an improved understanding of sulfenylation thermodynamics by showing that other residues from the protein environment other than the active site residues can play an important role.     

Effects of a distal mutation on active site chemistry

Previous studies of Escherichia coli dihydrofolate reductase (ecDHFR) have demonstrated that residue G121, which is 19 A from the catalytic center, is involved in catalysis, and long distance dynamical motions were implied. Specifically, the ecDHFR mutant G121V has been extensively studied by various experimental and theoretical tools, and the mutation's effect on kinetic, structural, and dynamical features of the enzyme has been explored. This work examined the effect of this mutation on the physical nature of the catalyzed hydride transfer step by means of intrinsic kinetic isotope effects (KIEs), their temperature dependence, and activation parameters as described previously for wild type ecDHFR [Sikorski, R. S., et al. (2004) J. Am. Chem. Soc. 126, 4778-4779]. The temperature dependence of initial velocities was used to estimate activation parameters. Isotope effects on the preexponential Arrhenius factors, and the activation energy, could be rationalized by an environmentally coupled hydrogen tunneling model, similar to the one used for the wild-type enzyme. Yet, in contrast to that in the wild type, fluctuations of the donor-acceptor distance were now required. Secondary (2 degrees ) KIEs were also measured for both H- and D-transfer, and as in the case of the wild-type enzyme, no coupled motion was detected. Despite these similarities, the reduced rates, the slightly inflated primary (1 degrees ) KIEs, and their temperature dependence, together with relatively deflated 2 degrees KIEs, indicate that the potential surface prearrangement was not as ideal as for the wild-type enzyme. These findings support theoretical studies suggesting that the G121V mutation led to a different conformational ensemble of reactive states and less effective rearrangement of the potential surface but has an only weak effect on H-tunneling.     

Exploring the RING-Catalyzed Ubiquitin Transfer Mechanism by MD and QM/MM Calculations

Ubiquitylation is a universal mechanism for controlling cellular functions. A large family of ubiquitin E3 ligases (E3) mediates Ubiquitin (Ub) modification. To facilitate Ub transfer, RING E3 ligases bind both the substrate and ubiquitin E2 conjugating enzyme (E2) linked to Ub via a thioester bond to form a catalytic complex. The mechanism of Ub transfer catalyzed by RING E3 remains elusive. By employing a combined computational approach including molecular modeling, molecular dynamics (MD) simulations, and quantum mechanics/molecular mechanics (QM/MM) calculations, we characterized this catalytic mechanism in detail. The three-dimensional model of dimeric RING E3 ligase RNF4 RING, E2 ligase UbcH5A, Ub and the substrate SUMO2 shows close contact between the substrate and Ub transfer catalytic center. Deprotonation of the substrate lysine by D117 on UbcH5A occurs with almost no energy barrier as calculated by MD and QM/MM calculations. Then, the side chain of the activated lysine gets close to the thioester bond via a conformation change. The Ub transfer pathway begins with a nucleophilic addition that forms an oxyanion intermediate of a 4.23 kcal/mol energy barrier followed by nucleophilic elimination, resulting in a Ub modified substrate by a 5.65 kcal/mol energy barrier. These results provide insight into the mechanism of RING-catalyzed Ub transfer guiding the discovery of Ub system inhibitors.     

An Alternative Mechanism for the Methylation of Phosphoethanolamine Catalyzed by Plasmodium falciparum Phosphoethanolamine Methyltransferase

The phosphobase methylation pathway catalyzed by the phosphoethanolamine methyltransferase in Plasmodium falciparum (PfPMT), the malaria parasite, offers an attractive target for anti-parasitic drug development. PfPMT methylates phosphoethanolamine (pEA) to phosphocholine for use in membrane biogenesis. Quantum mechanics and molecular mechanics (QM/MM) calculations tested the proposed reaction mechanism for methylation of pEA involving the previously identified Tyr-19–His-132 dyad, which indicated an energetically unfavorable mechanism. Instead, the QM/MM calculations suggested an alternative mechanism involving Asp-128. The reaction coordinate involves the stepwise transfer of a proton to Asp-128 via a bridging water molecule followed by a typical S_n2-type methyl transfer from S-adenosylmethionine to pEA. Functional analysis of the D128A, D128E, D128Q, and D128N PfPMT mutants shows a loss of activity with pEA but not with the final substrate of the methylation pathway. X-ray crystal structures of the PfPMT-D128A mutant in complex with S-adenosylhomocysteine and either pEA or phosphocholine reveal how mutation of Asp-128 disrupts a hydrogen bond network in the active site. The combined QM/MM, biochemical, and structural studies identify a key role for Asp-128 in the initial step of the phosphobase methylation pathway in Plasmodium and provide molecular insight on the evolution of multiple activities in the active site of the PMT.     

Application of a SCC-DFTB QM/MM approach to the investigation of the catalytic mechanism of fatty acid amide hydrolase

Self-consistent charge density functional tight binding (SCC-DFTB) is a promising method for hybrid quantum mechanics/molecular mechanics (QM/MM) simulations of enzyme-catalyzed reactions. The acylation reaction of fatty acid amide hydrolase (FAAH), a promising drug target, was investigated by applying a SCC-DFTB/CHARMM27 scheme. Calculated potential energy barriers resulted in reasonable agreement with experiments for oleamide (OA) and oleoylmethyl ester (OME) substrates, outperforming previous calculations performed at the PM3/CHARMM22 level. Furthermore, the experimental preference of FAAH in hydrolyzing OA faster than OME was adequately reproduced by calculations. All these findings indicate that the SCC-DFTB/CHARMM27 approach can be successfully applied to mechanistic investigations of FAAH-catalyzed reactions.     

Environmentally coupled hydrogen tunneling. Linking catalysis to dynamics.

Many biological C-H activation reactions exhibit nonclassical kinetic isotope effects (KIEs). These nonclassical KIEs are too large (kH/kD > 7) and/or exhibit unusual temperature dependence such that the Arrhenius prefactor KIEs (AH/AD) fall outside of the semiclassical range near unity. The focus of this minireview is to discuss such KIEs within the context of the environmentally coupled hydrogen tunneling model. Full tunneling models of hydrogen transfer assume that protein or solvent fluctuations generate a reactive configuration along the classical, heavy-atom coordinate, from which the hydrogen transfers via nuclear tunneling. Environmentally coupled tunneling also invokes an environmental vibration (gating) that modulates the tunneling barrier, leading to a temperature-dependent KIE. These properties directly link enzyme fluctuations to the reaction coordinate for hydrogen transfer, making a quantum view of hydrogen transfer necessarily a dynamic view of catalysis. The environmentally coupled hydrogen tunneling model leads to a range of magnitudes of KIEs, which reflect the tunneling barrier, and a range of AH/AD values, which reflect the extent to which gating modulates hydrogen transfer. Gating is the primary determinant of the temperature dependence of the KIE within this model, providing insight into the importance of this motion in modulating the reaction coordinate. The potential use of variable temperature KIEs as a direct probe of coupling between environmental dynamics and the reaction coordinate is described. The evolution from application of a tunneling correction to a full tunneling model in enzymatic H transfer reactions is discussed in the context of a thermophilic alcohol dehydrogenase and soybean lipoxygenase-1.