Molecular dynamics study of taxadiene synthase catalysis

Molecular dynamics (MD) simulations have been performed to study the dynamic behavior of noncovalent enzyme carbocation complexes involved in the cyclization of geranylgeranyl diphosphate to taxadiene catalyzed by taxadiene synthase (TXS). Taxadiene and the observed four side products originate from the deprotonation of carbocation intermediates. The MD simulations of the TXS carbocation complexes provide insights into potential deprotonation mechanisms of such carbocations. The MD results do not support a previous hypothesis that carbocation tumbling is a key factor in the deprotonation of the carbocations by pyrophosphate. Instead water bridges are identified which may allow the formation of side products via multiple proton transfer reactions. A novel reaction path for taxadiene formation is proposed on the basis of the simulations. © 2018 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.     

Studies on taxadiene synthase: interception of the cyclization cascade at the isocembrene stage with GGPP analogues

[reaction: see text] The cyclization of GGPP to taxadiene, catalyzed by taxadiene synthase, has been suggested to proceed through a series of monocyclic isocembrenyl- and bicyclic verticillyl-carbocationic intermediary stages. A set of GGPP analogues with abolished or perturbed pi-nucleophilicity at the delta10 double bond (GGPP numbering) was synthesized and incubated with taxadiene synthase to intercept the cyclization cascade at the monocyclic stage. Each analogue was transformed by taxadiene synthase in vitro to hydrocarbon products in varying yields, and the structures of the major product in each reaction were solved by GCEIMS and one- and two-dimensional (1H and 13C) NMR and found to be 14-membered monocyclic isocembrenyl diterpenes, indicating that the first C-C bond formation catalyzed by taxadiene synthase could be uncoupled from the other subsequent bond formation events by using suitably designed substrate analogues. The formation and isolation of these isocembrenyl diterpene products using taxadiene synthase supports proposals that the isocembrenyl cation is an intermediate in the cyclization of GGPP to taxadiene.     

Quantum Mechanics/Molecular Mechanics Simulations Show Saccharide Distortion is Required for Reaction in Hen Egg-White Lysozyme

Hybrid quantum mechanics/molecular mechanics (QM/MM) calculations on lysozyme show significant distortion of the bound saccharide is required to facilitate the catalytic reaction.     

Quantum Mechanics/Molecular Mechanics Modeling of Enzymatic Processes: Caveats and Breakthroughs

Nature has developed large groups of enzymatic catalysts with the aim to transfer substrates into useful products, which enables biosystems to perform all their natural functions. As such, all biochemical processes in our body (we drink, we eat, we breath, we sleep, etc.) are governed by enzymes. One of the problems associated with research on biocatalysts is that they react so fast that details of their reaction mechanisms cannot be obtained with experimental work. In recent years, major advances in computational hardware and software have been made and now large (bio)chemical systems can be studied using accurate computational techniques. One such technique is the quantum mechanics/molecular mechanics (QM/MM) technique, which has gained major momentum in recent years. Unfortunately, it is not a black‐box method that is easily applied, but requires careful set‐up procedures. In this work we give an overview on the technical difficulties and caveats of QM/MM and discuss work‐protocols developed in our groups for running successful QM/MM calculations.     

QM/MM Studies of the Matrix Metalloproteinase 2 (MMP2) Inhibition Mechanism of (S)-SB-3CT and its Oxirane Analogue

SB-3CT, (4-phenoxyphenylsulfonyl)methylthiirane, is a potent, mechanism-based inhibitor of the gelatinase sub-class of the matrix metalloproteinase (MMP) family of zinc proteases. The gelatinase MMPs are unusual in that there are several examples where both enantiomers of a racemic inhibitor have comparable inhibitory abilities. SB-3CT is one such example. Here, the inhibition mechanism of the MMP2 gelatinase by the (S)-SB-3CT enantiomer and its oxirane analogue is examined computationally, and compared to the mechanism of (R)-SB-3CT. Inhibition of MMP2 by (R)-SB-3CT was shown previously to involve enzyme-catalyzed C-H deprotonation adjacent to the sulfone, with concomitant opening by β-elimination of the sulfur of the three-membered thiirane ring. Similarly to the R enantiomer, (S)-SB-3CT was docked into the active site of MMP2, followed by molecular dynamics simulation to prepare the complex for combined quantum mechanics and molecular mechanics (QM/MM) calculations. QM/MM calculations with B3LYP/6-311+G(d,p) for the QM part (46 atoms) and the AMBER force field for the MM part were used to compare the reaction of (S)-SB-3CT and its oxirane analogue in the active site of MMP2 (9208 atoms). These calculations show that the barrier for the proton abstraction coupled ring opening reaction of (S)-SB-3CT in the MMP2 active site is 4.4 kcal/mol lower than its oxirane analogue, and the ring opening reaction energy of (S)-SB-3CT is only 1.6 kcal/mol less exothermic than its oxirane analogue. Calculations also show that the protonation of the ring-opened products by water is thermodynamically much more favorable for the alkoxide obtained from the oxirane, than for the thiolate obtained from the thiirane. In contrast to (R)-SB-3CT and the R-oxirane analogue, the double bonds of the ring-opened products of (S)-SB-3CT and its S-oxirane analogue have the cis-configuration. Vibrational frequency and intrinsic reaction path calculations on a reduced size QM/MM model (2747 atoms) provide additional insight into the mechanism. These calculations yield 5.9 and 6.7 for the deuterium kinetic isotope effect for C-H bond cleavage in the transition state for the R and S enantiomers of SB-3CT, in good agreement with the experimental results.     

QM/MM modeling of benzene hydroxylation in human cytochrome P450 2C9

The mechanism of benzene hydroxylation was investigated in the realistic enzyme environment of the human CYP 2C9 by using quantum mechanical/molecular mechanical (QM/MM) calculations of the whole reaction profile using the B3LYP method to describe the QM region. The calculated QM/MM barriers for addition of the active species Compound I to benzene are consistent with experimental rate constants for benzene metabolism in CYP 2E1. In contrast to gas-phase model calculations, our results suggest that competing side-on and face-on geometries of arene addition may both occur in the case of aromatic ring oxidation in cytochrome P450s. QM/MM profiles for three different rearrangement pathways of the initially formed sigma-adduct, leading to formation of epoxide, ketone, and an N-protonated porphyrin species, were calculated. Our results suggest that epoxide and ketone products form with comparable ease in the face-on pathway, whereas epoxide formation is preferred in the side-on pathway. Additionally, rearrangement to the N-protonated porphyrin species was found to be competitive with side-on epoxide formation. This suggests that overall, the competition between formation of epoxide and phenol final products in P450 oxidation of aromatic substrates is quite finely balanced.     

Taxadiene synthase-catalyzed cyclization of 6-fluorogeranylgeranyl diphosphate to 7-fluoroverticillenes

The mechanism of the taxadiene synthase-catalyzed cyclization of (E,E,E)-geranylgeranyl diphosphate (GGPP, 7) to taxadiene (5) is proposed to proceed through a verticillen-12-yl carbocation intermediate (8) that undergoes an 11 --> 7 proton transfer leading to formation of the C ring. The substrate analogue 6-fluoroGGPP (17) was synthesized to elucidate the stereochemistry of the putative verticillenyl intermediate. It was expected that the inductive electron-withdrawing effect of the fluoro substituent would prevent the critical proton transfer to the Delta(7) double bond and thereby derail the cyclization at the bicyclic stage. Incubation of the fluoro analogue with recombinant taxadiene synthase yielded a mixture of three major and two minor fluoro diterpenes according to GC/MS analyses. The three major products were identified as the exocyclic, endocyclic, and 4(20)-methylene 7-fluoroverticillenes, i.e., Delta(3,7,12 (18)), Delta(3,7,12), and Delta(4(20),7,11) isomers (22, 23, and 24) on the basis of (1)H NMR analyses and comparisons with the parent bicyclic diterpenes. The H1beta, H11alpha (1S,11R) configurations at the bridgehead positions of 22 were established by means of NOE experiments and CD spectra. The absolute configuration of (+)-verticillol (4) was revised after the anomalous dispersion X-ray analysis of (+)-verticillol p-iodobenzoate. Of particular note, all absolute configurations of verticillane diterpenes in the literature should be reversed. This work affords compelling evidence supporting the H11alpha (11R) stereochemistry of the verticillen-12-yl(+) ion intermediate in the taxadiene synthase-catalyzed reaction and illustrates the capability of vinyl fluoro analogues to intercept complex cyclization cascades.     

Assessing the validity of DLPNO-CCSD(T) in the calculation of activation and reaction energies of ubiquitous enzymatic reactions

The domain-based local pair natural orbital coupled-cluster with single, double, and perturbative triples excitation (DLPNO-CCSD(T)) method was employed to portray the activation and reaction energies of four ubiquitous enzymatic reactions, and its performance was confronted to CCSD(T)/complete basis set (CBS) to assess its accuracy and robustness in this specific field. The DLPNO-CCSD(T) results were also confronted to those of a set of density functionals (DFs) to understand the benefit of implementing this technique in enzymatic quantum mechanics/molecular mechanics (QM/MM) calculations as a second QM component, which is often treated with DF theory (DFT). On average, the DLPNO-CCSD(T)/aug-cc-pVTZ results were 0.51 kcal·mol⁻¹ apart from the canonic CCSD(T)/CBS, without noticeable biases toward any of the reactions under study. All DFs fell short to the DLPNO-CCSD(T), both in terms of accuracy and robustness, which suggests that this method is advantageous to characterize enzymatic reactions and that its use in QM/MM calculations, either alone or in conjugation with DFT, in a two-region QM layer (DLPNO-CCSD(T):DFT), should enhance the quality and faithfulness of the results.     

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.     

YinYang atom: a simple combined ab initio quantum mechanical molecular mechanical model

A simple interface is proposed for combined quantum mechanical (QM) molecular mechanical (MM) calculations for the systems where the QM and MM regions are connected through covalent bonds. Within this model, the atom that connects the two regions, called YinYang atom here, serves as an ordinary MM atom to other MM atoms and as a hydrogen-like atom to other QM atoms. Only one new empirical parameter is introduced to adjust the length of the connecting bond and is calibrated with the molecule propanol. This model is tested with the computation of equilibrium geometries and protonation energies for dozens of molecules. Special attention is paid on the influence of MM point charges on optimized geometry and protonation energy, and it is found that it is important to maintain local charge-neutrality in the MM region in order for the accurate calculation of the protonation and deprotonation energies. Overall the simple YinYang atom model yields comparable results to some other QM/MM models.     

Geometry optimization with QM/MM,ONIOM, and other combined methods. I. Microiterations and constraints

Hybrid energy methods such as QM/MM and ONIOM, that combine different levels of theory into one calculation, have been very successful in describing large systems. Geometry optimization methods can take advantage of the partitioning of these calculations into a region treated at a quantum mechanical (QM) level of theory and the larger, remaining region treated by an inexpensive method such as molecular mechanics (MM). A series of microiterations can be employed to fully optimize the MM region for each optimization step in the QM region. Cartesian coordinates are used for the MM region and are chosen so that the internal coordinates of the QM region remain constant during the microiterations. The coordinates of the MM region are augmented to permit rigid body translation and rotation of the QM region. This is essential if any atoms in the MM region are constrained, but it also improves the efficiency of unconstrained optimizations. Because of the microiterations, special care is needed for the optimization step in the QM region so that the system remains in the same local valley during the course of the optimization. The optimization methodology with microiterations, constraints, and step-size control are illustrated by calculations on bacteriorhodopsin and other systems.     

Quantum mechanical/molecular mechanical study on the mechanism of the enzymatic Baeyer-Villiger reaction

We report a combined quantum mechanical/molecular mechanical (QM/MM) study on the mechanism of the enzymatic Baeyer-Villiger reaction catalyzed by cyclohexanone monooxygenase (CHMO). In QM/MM geometry optimizations and reaction path calculations, density functional theory (B3LYP/TZVP) is used to describe the QM region consisting of the substrate (cyclohexanone), the isoalloxazine ring of C4a-peroxyflavin, the side chain of Arg-329, and the nicotinamide ring and the adjacent ribose of NADP(+), while the remainder of the enzyme is represented by the CHARMM force field. QM/MM molecular dynamics simulations and free energy calculations at the semiempirical OM3/CHARMM level employ the same QM/MM partitioning. According to the QM/MM calculations, the enzyme-reactant complex contains an anionic deprotonated C4a-peroxyflavin that is stabilized by strong hydrogen bonds with the Arg-329 residue and the NADP(+) cofactor. The CHMO-catalyzed reaction proceeds via a Criegee intermediate having pronounced anionic character. The initial addition reaction has to overcome an energy barrier of about 9 kcal/mol. The formed Criegee intermediate occupies a shallow minimum on the QM/MM potential energy surface and can undergo fragmentation to the lactone product by surmounting a second energy barrier of about 7 kcal/mol. The transition state for the latter migration step is the highest point on the QM/MM energy profile. Gas-phase reoptimizations of the QM region lead to higher barriers and confirm the crucial role of the Arg-329 residue and the NADP(+) cofactor for the catalytic efficiency of CHMO. QM/MM calculations for the CHMO-catalyzed oxidation of 4-methylcyclohexanone reproduce and rationalize the experimentally observed (S)-enantioselectivity for this substrate, which is governed by the conformational preferences of the corresponding Criegee intermediate and the subsequent transition state for the migration step.     

Toward a consistent treatment of polarization in model QM/MM calculations

The concept of model chemistries within hybrid QM/MM calculations has been addressed through analysis of the polarization energy determined by two distinct approaches based on (i) induced charges and (ii) induced dipoles. The quantum mechanical polarization energy for four configurations of the water dimer has been determined for a range of basis sets using Morokuma energy decomposition analysis. This benchmark value has been compared to the fully classical polarization energy determined using the induced dipole approach, and the molecular mechanics polarization energy calculated using induced charges within the MM region of hybrid QM/MM calculations. From the water dimer calculations, it is concluded that the induced charge approach is consistent with medium sized basis set calculations whereas the induced dipole approach is consistent with large basis set calculations. This result is highly relevant to the concept of QM/MM model chemistries.     

Spin-projected QM/MM Free Energy Simulations for Oxidation Reaction of Guanine in B−DNA by Singlet Oxygen

Guanine is the most susceptible base to oxidation damage induced by reactive oxygen species including singlet oxygen (¹O₂, ¹Δ_g). We clarify whether the first step of guanine oxidation in B−DNA proceeds via either a zwitterionic or a diradical intermediate. The free energy profiles are calculated by means of a combined quantum mechanical and molecular mechanical (QM/MM) method coupled with the adaptive biasing force (ABF) method. To describe the open-shell electronic structure of ¹O₂ correctly, the broken-symmetry spin-unrestricted density functional theory (BS−UDFT) with an approximate spin projection (AP) correction is applied to the QM region. We find that the effect of spin contamination on the activation and reaction free energies is up to ∼8 kcal mol⁻¹, which is too large to be neglected. The QM(AP−ULC−BLYP)/MM-based free energy calculations also reveal that the reaction proceeds through a diradical transition state, followed by a conversion to a zwitterionic intermediate. Our computed activation energy of 5.2 kcal mol⁻¹ matches experimentally observed range (0∼6 kcal mol⁻¹).     

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.     

Glutathione transferase: new model for glutathione activation

Glutathione transferases are enzymes of the cellular detoxification system that metabolize a vast spectrum of xenobiotic and endobiotic toxic compounds. They are homodimers or heterodimers and each monomer has an active center composed of a G-site in which glutathione (GSH) binds and an H-site for the electrophilic substrate. When GSH binds to the G-site, the pKa value of its thiol group drops by 2.5 units; this promotes its deprotonation and, therefore, produces a strong nucleophilic thiolate that is able to react with the electrophilic substrate. The mechanism behind the deprotonation of the thiol group is still unknown. Some studies point to the fact that the GSH glutamyl alpha-carboxylate group is essential for GSH activation, whereas others indicate the importance of the active-center water molecules. On the basis of QM/MM calculations, we propose a mechanism of GSH activation in which a water molecule, acting as a bridge, is able to assist in the transfer of the proton from the GSH thiol group to the GSH glutamyl alpha-carboxylate group, after an initial GSH conformational rearrangement. We calculated the potential of mean force of this GSH structural rearrangement that would be necessary for the approach of both groups and we then performed a QM/MM ONIOM scan of water-assisted proton transfer. The overall free-energy barrier for the process is consistent with experimental studies of the enzyme kinetics.     

Calculating distribution coefficients based on multi-scale free energy simulations: an evaluation of MM and QM/MM explicit solvent simulations of water-cyclohexane transfer in the SAMPL5 challenge

One of the central aspects of biomolecular recognition is the hydrophobic effect, which is experimentally evaluated by measuring the distribution coefficients of compounds between polar and apolar phases. We use our predictions of the distribution coefficients between water and cyclohexane from the SAMPL5 challenge to estimate the hydrophobicity of different explicit solvent simulation techniques. Based on molecular dynamics trajectories with the CHARMM General Force Field, we compare pure molecular mechanics (MM) with quantum-mechanical (QM) calculations based on QM/MM schemes that treat the solvent at the MM level. We perform QM/MM with both density functional theory (BLYP) and semi-empirical methods (OM1, OM2, OM3, PM3). The calculations also serve to test the sensitivity of partition coefficients to solute polarizability as well as the interplay of the quantum-mechanical region with the fixed-charge molecular mechanics environment. Our results indicate that QM/MM with both BLYP and OM2 outperforms pure MM. However, this observation is limited to a subset of cases where convergence of the free energy can be achieved.