High-Accuracy Computation of Reaction Barriers in Enzymes

The Mechanism of Biochemical NO‐Sensing: Insights from Computational Chemistry

The binding of small gas molecules such as NO and CO plays a major role in the signaling routes of the human body. The sole NO‐receptor in humans is soluble guanylyl cyclase (sGC) – a histidine‐ligated heme protein, which, upon NO binding, activates a downstream signaling cascade. Impairment of NO‐signaling is linked, among others, to cardiovascular and inflammatory diseases. In the present work, we use a combination of theoretical tools such as MD simulations, high‐level quantum chemical calculations and hybrid QM/MM methods to address various aspects of NO binding and to elucidate the most likely reaction paths and the potential intermediates of the reaction. As a model system, the H‐NOX protein from Shewanella oneidensis (So H‐NOX) homologous to the NO‐binding domain of sGC is used. The signaling route is predicted to involve NO binding to form a six‐coordinate intermediate heme‐NO complex, followed by relatively facile His decoordination yielding a five‐coordinate adduct with NO on the distal side with possible isomerization to the proximal side through binding of a second NO and release of the first one. MD simulations show that the His sidechain can quite easily rotate outward into solvent, with this motion being accompanied in our simulations by shifts in helix positions that are consistent with this decoordination leading to significant conformational change in the protein.

Computational Study of the Peroxygenase Mechanism Catalyzed by Hemoglobin Dehaloperoxidase Involved in the Degradation of Chlorophenols

The biochemical evidence showed that hemoglobin dehaloperoxidase (DHP B) from Amphitrite Ornata is a multifunctional hemoprotein that catalyzes both dehalogenation and hydroxylation of halophenols via the peroxidase and peroxygenase mechanism, respectively, which sets the basis for the degradation of halophenols. In the peroxygenase mechanism, the reaction was previously suggested to be triggered either by the hydrogen atom abstraction by the Fe═O center or by the proton abstraction by His55. To illuminate the peroxygenase mechanism of DHP B at the atomistic level, on the basis of the high-resolution crystal structure, computational models were constructed, and a series of quantum mechanical/molecular mechanical calculations have been performed. According to the calculation results, the pathway (Path a) initiated by the H-abstraction by the Fe═O center is feasible. In another pathway (Path b), His55 can abstract the proton from the hydroxyl group of the substrate (4-Cl-o-cresol) to initiate the reaction; however, its feasibility depends on the prior electron transfer from the substrate to the porphyrin group. The rate-limiting step of Path a is the OH-rebound, which corresponds to an energy barrier of 14.7 kcal/mol at the quartet state. His55 acts as an acid-base catalyst and directly involves in the catalysis. Our mutant study indicates that His55 can be replaced by other titratable residues. These findings may provide useful information for further understanding of the catalytic reaction of DHP B and for the design of enzymes in the degradation of pollutants, in particular, halophenols.

Comparative QM/MM Molecular Dynamics and Umbrella Sampling Simulations: Interaction of the Zinc-Bound Intermediate Gem-Diolate Trapoxin A Inhibitor and Acetyl-l-lysine Substrate with Histone Deacetylase 8

Targeting the genetic material without destruction is a priority to develop safe anticancer drugs. Histone deacetylase 8 (HDAC8), which is proved to be involved in carcinogenesis, is an enzyme associated with the chromatin for post-translational deacetylation of acetylated lysine. In this study, HDAC8 co-crystallized with the intermediate state tetrapeptide Trapoxin A (TA) inhibitor and the holoenzyme are utilized to find their conformational ensembles. Furthermore, the co-crystallized intermediate gem-diolate TA was used to find optimum interactions with the active site residues by conventional molecular dynamics (MD) simulation and QM/MM umbrella sampling. Finally, the intermediate state of the acetyl-l-lysine substrate was explored by QM/MM steered MD and compared to the binding of the intermediate state of the inhibitor. This research showed that HDAC8 is flexible and exists in conformational ensembles in its holoenzyme state. Binding of the intermediate state TA stabilizes its conformation. The optimum binding to the active site of HDAC8 for structures of gem-diolate TA (intermediate state) and acetyl-l-lysine (intermediate state) was determined according to the corresponding energy profiles. The use of these models will aid in the design of potentially reversible, potent, and selective inhibitors of HDAC8 for cancer treatment.

Accurate Quantum Chemical Free Energies at Affordable Cost

Free energy sampling methods allow studying the full dynamics of activated processes. Unfortunately, the affordable accuracy of the potential describing the energy and forces of the system is usually rather low. Here we introduce a new method that by combining metadynamics and free energy perturbation allows calculating accurate quantum chemical free energies for chemical reactions. To prove the effectiveness of this new approach we study the S_N2 reaction of CH₃F + Cl^- → CH₃Cl + F^- in vacuo and solvated by water. Comparisons are made with harmonic transition-state theory to show how this method could provide accurate equilibrium and rate constants for complex systems.

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.

Uncertainty quantification for quantum chemical models of complex reaction networks

For the quantitative understanding of complex chemical reaction mechanisms, it is, in general, necessary to accurately determine the corresponding free energy surface and to solve the resulting continuous-time reaction rate equations for a continuous state space. For a general (complex) reaction network, it is computationally hard to fulfill these two requirements. However, it is possible to approximately address these challenges in a physically consistent way. On the one hand, it may be sufficient to consider approximate free energies if a reliable uncertainty measure can be provided. On the other hand, a highly resolved time evolution may not be necessary to still determine quantitative fluxes in a reaction network if one is interested in specific time scales. In this paper, we present discrete-time kinetic simulations in discrete state space taking free energy uncertainties into account. The method builds upon thermo-chemical data obtained from electronic structure calculations in a condensed-phase model. Our kinetic approach supports the analysis of general reaction networks spanning multiple time scales, which is here demonstrated for the example of the formose reaction. An important application of our approach is the detection of regions in a reaction network which require further investigation, given the uncertainties introduced by both approximate electronic structure methods and kinetic models. Such cases can then be studied in greater detail with more sophisticated first-principles calculations and kinetic simulations.

A combined density functional theory and numerical simulation investigation of levels of chirality transfer and regioselectivity for the radical cyclizations of N-methyl-, N-ethyl- and N-isopropyl-substituted ortho-halo-N-acryloylanilides

There is a high chirality transfer ratio when the N-alkyl substituent is larger in volume than isopropyl.

QM/MM Study of the Conversion of Oxophlorin into Verdoheme by Heme Oxygenase

Heme oxygenase is an enzyme that degrades heme, thereby recycling iron in most organisms, including humans. Pervious density functional theory (DFT) calculations have suggested that iron(III) hydroxyheme, an intermediate generated in the first step of heme degradation by heme oxygenase, is converted to iron(III) superoxo oxophlorin in the presence of dioxygen. In this article, we have studied the detailed mechanism of conversion of iron(III) superoxo oxophlorin to verdoheme by using combined quantum mechanics and molecular mechanics (QM/MM) calculations. The calculations employed the B3LYP method and the def2-QZVP basis set, considering dispersion effects with the DFT-D3 approach, obtaining accurate energies with large QM regions of almost 1000 atoms. The reaction was found to be exothermic by -35 kcal/mol, with a rate-determining barrier of 19 kcal/mol in the doublet state. The protein environment and especially water in the enzyme pocket significantly affects the reaction by decreasing the reaction activation energies and changing the structures by providing strategic hydrogen bonds.

Tryptophan lyase (NosL): mechanistic insights into amine dehydrogenation and carboxyl fragment migration by QM/MM calculations

QM/MM calculations suggest two feasible pathways for the breaking of the C–C bond of the substrate. The breaking of the Cα–Cβ bond leads to the final product, whereas the cleavage of the Cα–C bond will terminate in the EPR-trapped radical intermediate.

What's My Substrate? Computational Function Assignment of Candida parapsilosis ADH5 by Genome Database Search, Virtual Screening, and QM/MM Calculations

Zinc-dependent medium chain reductase from Candida parapsilosis can be used in the reduction of carbonyl compounds to pharmacologically important chiral secondary alcohols. To date, the nomenclature of cpADH5 is differing (CPCR2/RCR/SADH) in the literature, and its natural substrate is not known. In this study, we utilized a substrate docking based virtual screening method combined with KEGG, MetaCyc pathway, and Candida genome databases search for the discovery of natural substrates of cpADH5. The virtual screening of 7834 carbonyl compounds from the ZINC database provided 94 aldehydes or methyl/ethyl ketones as putative carbonyl substrates. Out of which, 52 carbonyl substrates of cpADH5 with catalytically active docking pose were identified by employing mechanism based substrate docking protocol. Comparison of the virtual screening results with KEGG, MetaCyc database search, and Candida genome pathway analysis suggest that cpADH5 might be involved in the Ehrlich pathway (reduction of fusel aldehydes in leucine, isoleucine, and valine degradation). Our QM/MM calculations and experimental activity measurements affirmed that butyraldehyde substrates are the potential natural substrates of cpADH5, suggesting a carbonyl reductase role for this enzyme in butyraldehyde reduction in aliphatic amino acid degradation pathways. Phylogenetic tree analysis of known ADHs from Candida albicans shows that cpADH5 is close to caADH5. We therefore propose, according to the experimental substrate identification and sequence similarity, the common name butyraldehyde dehydrogenase cpADH5 for Candida parapsilosis CPCR2/RCR/SADH.

Ab Initio Calculation of Rate Constants for Molecule-Surface Reactions with Chemical Accuracy

The ab initio prediction of reaction rate constants for systems with hundreds of atoms with an accuracy that is comparable to experiment is a challenge for computational quantum chemistry. We present a divide-and-conquer strategy that departs from the potential energy surfaces obtained by standard density functional theory with inclusion of dispersion. The energies of the reactant and transition structures are refined by wavefunction-type calculations for the reaction site. Thermal effects and entropies are calculated from vibrational partition functions, and the anharmonic frequencies are calculated separately for each vibrational mode. This method is applied to a key reaction of an industrially relevant catalytic process, the methylation of small alkenes over zeolites. The calculated reaction rate constants (free energies), pre-exponential factors (entropies), and enthalpy barriers show that our computational strategy yields results that agree with experiment within chemical accuracy limits (less than one order of magnitude).

Acceleration of Ab Initio QM/MM Calculations under Periodic Boundary Conditions by Multiscale and Multiple Time Step Approaches

Development of multiscale ab initio quantum mechanical and molecular mechanical (AI-QM/MM) method for periodic boundary molecular dynamics (MD) simulations and their acceleration by multiple time step approach are described. The developed method achieves accuracy and efficiency by integrating the AI-QM/MM level of theory and the previously developed semiempirical (SE) QM/MM-Ewald sum method [J. Chem. Theory Comput. 2005, 1, 2] extended to the smooth particle-mesh Ewald (PME) summation method. In the developed methods, the total energy of the simulated system is evaluated at the SE-QM/MM-PME level of theory to include long-range QM/MM electrostatic interactions, which is then corrected on the fly using the AI-QM/MM level of theory within the real space cutoff. The resulting energy expression enables decomposition of total forces applied to each atom into forces determined at the low-level SE-QM/MM method and correction forces at the AI-QM/MM level, to integrate the system using the reversible reference system propagator algorithm. The resulting method achieves a substantial speed-up of the entire calculation by minimizing the number of time-consuming energy and gradient evaluations at the AI-QM/MM level. Test calculations show that the developed multiple time step AI-QM/MM method yields MD trajectories and potential of mean force profiles comparable to single time step QM/MM results. The developed method, together with message passing interface (MPI) parallelization, accelerates the present AI-QM/MM MD simulations about 30-fold relative to the speed of single-core AI-QM/MM simulations for the molecular systems tested in the present work, making the method less than one order slower than the SE-QM/MM methods under periodic boundary conditions.

Quantum mechanical modeling: a tool for the understanding of enzyme reactions

Most enzyme reactions involve formation and cleavage of covalent bonds, while electrostatic effects, as well as dynamics of the active site and surrounding protein regions, may also be crucial. Accordingly, special computational methods are needed to provide an adequate description, which combine quantum mechanics for the reactive region with molecular mechanics and molecular dynamics describing the environment and dynamic effects, respectively. In this review we intend to give an overview to non-specialists on various enzyme models as well as established computational methods and describe applications to some specific cases. For the treatment of various enzyme mechanisms, special approaches are often needed to obtain results, which adequately refer to experimental data. As a result of the spectacular progress in the last two decades, most enzyme reactions can be quite precisely treated by various computational methods.

In silico screening of 393 mutants facilitates enzyme engineering of amidase activity in CalB

Our previously presented method for high throughput computational screening of mutant activity (Hediger et al., 2012) is benchmarked against experimentally measured amidase activity for 22 mutants of Candida antarctica lipase B (CalB). Using an appropriate cutoff criterion for the computed barriers, the qualitative activity of 15 out of 22 mutants is correctly predicted. The method identifies four of the six most active mutants with ≥3-fold wild type activity and seven out of the eight least active mutants with ≤0.5-fold wild type activity. The method is further used to screen all sterically possible (386) double-, triple- and quadruple-mutants constructed from the most active single mutants. Based on the benchmark test at least 20 new promising mutants are identified.

Mapping enzymatic catalysis using the effective fragment molecular orbital method: towards all ab initio biochemistry

We extend the Effective Fragment Molecular Orbital (EFMO) method to the frozen domain approach where only the geometry of an active part is optimized, while the many-body polarization effects are considered for the whole system. The new approach efficiently mapped out the entire reaction path of chorismate mutase in less than four days using 80 cores on 20 nodes, where the whole system containing 2398 atoms is treated in the ab initio fashion without using any force fields. The reaction path is constructed automatically with the only assumption of defining the reaction coordinate a priori. We determine the reaction barrier of chorismate mutase to be kcal mol⁻¹ for MP2/cc-pVDZ and for MP2/cc-pVTZ in an ONIOM approach using EFMO-RHF/6-31G(d) for the high and low layers, respectively.

Benchmark Study for the Cysteine-Histidine Proton Transfer Reaction in a Protein Environment: Gas Phase, COSMO, QM/MM Approaches

Proton transfer reactions are of crucial interest for the investigation of proteins. We have investigated the accuracy of commonly used quantum chemical methods for the description of proton transfer reactions in different environments (gas phase, COSMO, QM/MM) using the proton transfer between the catalytic dyad residues cysteine 145 and histidine 41 of SARS coronavirus main protease as a case study. The test includes thermodynamic, kinetic, and structural properties. The study comprises computationally demanding ab initio approaches (HF, CC2, MP2, SCS-CC2, SCS-MP2, CCSD(T)), popular density functional theories (BLYP, B3LYP, M06-2X), and semiempirical methods (MNDO/d, AM1, RM1, PM3, PM6). The approximated coupled cluster approach LCCSD(T) is taken as a reference method. We find that the robustness of the tested methods with respect to the environment correlates well with the level of theory. As an example HF, CC2, MP2, and their SCS variants show similar errors for gas phase, COSMO, or QM/MM computations. In contrast for semiempirical methods, the errors strongly diversify if one goes from gas phase to COSMO or QM/MM. Particular problems are observed for the recent semiempirical methods PM6 and RM1, which show the best performance for gas phase calculations but possess larger errors in conjunction with COSMO. Finally, a combination of SCS-MP2 and B3LYP or M06-2X allows reliable estimates about remaining errors.

A computational methodology to screen activities of enzyme variants

We present a fast computational method to efficiently screen enzyme activity. In the presented method, the effect of mutations on the barrier height of an enzyme-catalysed reaction can be computed within 24 hours on roughly 10 processors. The methodology is based on the PM6 and MOZYME methods as implemented in MOPAC2009, and is tested on the first step of the amide hydrolysis reaction catalyzed by the Candida Antarctica lipase B (CalB) enzyme. The barrier heights are estimated using adiabatic mapping and shown to give barrier heights to within 3 kcal/mol of B3LYP/6-31G(d)//RHF/3-21G results for a small model system. Relatively strict convergence criteria (0.5 kcal/(molÅ)), long NDDO cutoff distances within the MOZYME method (15 Å) and single point evaluations using conventional PM6 are needed for reliable results. The generation of mutant structures and subsequent setup of the semiempirical calculations are automated so that the effect on barrier heights can be estimated for hundreds of mutants in a matter of weeks using high performance computing.

Optimization of orbital-specific virtuals in local Møller-Plesset perturbation theory

We present an orbital-optimized version of our orbital-specific-virtuals second-order Møller-Plesset perturbation theory (OSV-MP2). The OSV model is a local correlation ansatz with a small basis of virtual functions for each occupied orbital. It is related to the Pulay-Saebø approach, in which domains of virtual orbitals are drawn from a single set of projected atomic orbitals; but here the virtual functions associated with a particular occupied orbital are specifically tailored to the correlation effects in which that orbital participates. In this study, the shapes of the OSVs are optimized simultaneously with the OSV-MP2 amplitudes by minimizing the Hylleraas functional or approximations to it. It is found that optimized OSVs are considerably more accurate than the OSVs obtained through singular value decomposition of diagonal blocks of MP2 amplitudes, as used in our earlier work. Orbital-optimized OSV-MP2 recovers smooth potential energy surfaces regardless of the number of virtuals. Full optimization is still computationally demanding, but orbital optimization in a diagonal or Kapuy-type MP2 approximation provides an attractive scheme for determining accurate OSVs.

Correlating DFT-calculated energy barriers to experiments in nonheme octahedral Fe(IV)O species

The experimentally measured bimolecular reaction rate constant, k(2), should in principle correlate with the theoretically calculated rate-limiting free energy barrier, ΔG(≠), through the Eyring equation, but it fails quite often to do so due to the inability of current computational methods to account in a precise manner for all the factors contributing to ΔG(≠). This is further aggravated by the exponential sensitivity of the Eyring equation to these factors. We have taken herein a pragmatic approach for C-H activation reactions of 1,4-cyclohexadiene with a variety of octahedral nonheme Fe(IV)O complexes. The approach consists of empirically determining two constants that would aid in predicting experimental k(2) values uniformly from theoretically calculated electronic energy (ΔE(≠)) values. Shown in this study is the predictive power as well as insights into energy relationships in Fe(IV)O C-H activation reactions. We also find that the difference between ΔG(≠) and ΔE(≠) converges at slow reactions, in a manner suggestive of changes in the importance of the triplet spin state weight in the overall reaction.

QM/MM methods for biomolecular systems

Combined quantum-mechanics/molecular-mechanics (QM/MM) approaches have become the method of choice for modeling reactions in biomolecular systems. Quantum-mechanical (QM) methods are required for describing chemical reactions and other electronic processes, such as charge transfer or electronic excitation. However, QM methods are restricted to systems of up to a few hundred atoms. However, the size and conformational complexity of biopolymers calls for methods capable of treating up to several 100,000 atoms and allowing for simulations over time scales of tens of nanoseconds. This is achieved by highly efficient, force-field-based molecular mechanics (MM) methods. Thus to model large biomolecules the logical approach is to combine the two techniques and to use a QM method for the chemically active region (e.g., substrates and co-factors in an enzymatic reaction) and an MM treatment for the surroundings (e.g., protein and solvent). The resulting schemes are commonly referred to as combined or hybrid QM/MM methods. They enable the modeling of reactive biomolecular systems at a reasonable computational effort while providing the necessary accuracy.