High-accuracy computation of reaction barriers in enzymes

CHARMM-GUI QM/MM Interfacer for a Quantum Mechanical and Molecular Mechanical (QM/MM) Simulation Setup: 1. Semiempirical Methods

Quantum mechanical (QM) treatments, when combined with molecular mechanical (MM) force fields, can effectively handle enzyme-catalyzed reactions without significantly increasing the computational cost. In this context, we present CHARMM-GUI QM/MM Interfacer, a web-based cyberinfrastructure designed to streamline the preparation of various QM/MM simulation inputs with ligand modification. The development of QM/MM Interfacer has been achieved through integration with existing CHARMM-GUI modules, such as PDB Reader and Manipulator, Solution Builder, and Membrane Builder. In addition, new functionalities have been developed to facilitate the one-stop preparation of QM/MM systems and enable interactive and intuitive ligand modifications and QM atom selections. QM/MM Interfacer offers support for a range of semiempirical QM methods, including AM1(+/d), PM3(+/PDDG), MNDO(+/d, +/PDDG), PM6, RM1, and SCC-DFTB, tailored for both AMBER and CHARMM. A nontrivial setup related to ligand modification, link-atom insertion, and charge distribution is automatized through intuitive user interfaces. To illustrate the robustness of QM/MM Interfacer, we conducted QM/MM simulations of three enzyme-substrate systems: dihydrofolate reductase, insulin receptor kinase, and oligosaccharyltransferase. In addition, we have created three tutorial videos about building these systems, which can be found at https://www.charmm-gui.org/demo/qmi. QM/MM Interfacer is expected to be a valuable and accessible web-based tool that simplifies and accelerates the setup process for hybrid QM/MM simulations.

Multilevel quantum mechanical calculations show the role of promoter molecules in the dehydration of methanol to dimethyl ether in H-ZSM-5

Methyl carboxylate esters promote the formation of dimethyl ether (DME) from the dehydration of methanol in H-ZSM-5 zeolite. We employ a multilevel quantum method to explore the possible associative and dissociative mechanisms in the presence, and absence, of six methyl ester promoters. This hybrid method combines density functional theory, with dispersion corrections (DFT-D3), for the full periodic system, with second-order Møller-Plesset perturbation theory (MP2) for small clusters representing the reaction site, and coupled cluster with single, double, and perturbative triple substitution (CCSD(T)) for the reacting molecules. The calculated adsorption enthalpy of methanol, and reaction enthalpies of the dehydration of methanol to DME within H-ZSM-5, agree with experiment to within chemical accuracy (∼4 kJ mol-1). For the promoters, a reaction pathway via an associative mechanism gives lower overall reaction enthalpies and barriers compared to the reaction with methanol only. Each stage of this mechanism is explored and related to experimental data. We provide evidence that suggests the promoter's adsorption to the Brønsted acid site is the most important factor dictating its efficiency.

Nonequilibrium Wet-Dry Cycling Acts as a Catalyst for Chemical Reactions

Recent experimental studies suggest that wet-dry cycles and coexisting phases can each strongly alter chemical processes. The mechanisms of why and to what degree chemical processes are altered when subjected to evaporation and condensation are unclear. To close this gap, we developed a theoretical framework for nondilute chemical reactions subject to nonequilibrium conditions of evaporation and condensation. We find that such conditions can change the half-time of the product's yield by more than an order of magnitude, depending on the substrate-solvent interaction. We show that the cycle frequency strongly affects the chemical turnover when the system is maintained out of equilibrium by wet-dry cycles. There exists a resonance behavior in the cycle frequency where the turnover is maximal. This resonance behavior enables wet-dry cycles to select specific chemical reactions, suggesting a potential mechanism for chemical evolution in prebiotic soups at early Earth.

Exploration of biochemical reactivity with a QM/MM growing string method

In this work, we have implemented the single-ended growing string method using a hybrid internal/Cartesian coordinate scheme within our in-house QM/MM package, QoMMMa, representing the first implementation of the growing string method in the QM/MM framework. The goal of the implementation was to facilitate generation of QM/MM reaction pathways with minimal user input, and also to improve the quality of the pathways generated as compared to the widely used adiabatic mapping approach. We have validated the algorithm against a reaction which has been studied extensively in previous computational investigations - the Claisen rearrangement catalysed by chorismate mutase. The nature of the transition state and the height of the barrier was predicted well using our algorithm, where more than 88% of the pathways generated were deemed to be of production quality. Directly compared to using adiabatic mapping, we found that while our QM/MM single-ended growing string method is slightly less efficient, it readily produces reaction pathways with fewer discontinuites and thus minimises the need for involved remapping of unsatisfactory energy profiles.

Accurate Calculation of Isomerization and Conformational Energies of Larger Molecules Using Explicitly Correlated Local Coupled Cluster Methods in Molpro and ORCA

An overview of the approximations in the explicitly correlated local coupled cluster methods PNO-LCCSD(T)-F12 in Molpro and DLPNO-CCSD(T)F12 in ORCA is given. Options to select the domains of projected atomic orbitals (PAOs), pair natural orbitals (PNOs), and triples natural orbitals (TNOs) in both programs are described and compared in detail. The two programs are applied to compute isomerization and conformational energies of the ISOL24 and ACONFL test sets, where the former is part of the GMTKN55 benchmark suite. Thorough studies of basis set effects are presented for selected systems. These revealed large intramolecular basis set superposition effects that make it practically impossible to reliably determine the complete basis set (CBS) limits without including explicitly correlated terms. The latter strongly reduce the basis set dependence and at the same time also errors caused by the local domain approximations. On the basis of these studies, the PNO-LCCSD(T)-F12 method is applied to determine new reference energies for the above-mentioned benchmark sets. We are confident that our results should agree within a few tenths of a kcal mol-1 with the (unknown) CCSD(T)/CBS values, which therefore allowed us to define computational settings for accurate explicitly correlated local coupled cluster methods with moderate computational effort. With these protocols, especially PNO-LCCSD(T)-F12b/AVTZ', reliable reference values for comprehensive benchmark sets can be generated efficiently. This can significantly advance the development and evaluation of the performance of approximate electronic structure methods, especially improved density functional approximations or machine learning approaches.

Enzymatic reconstitution of salicylate formation in promysalin biosynthesis

Promysalin is an amphipathic antibiotic isolated from Pseudomonas promysalinigenes (previously Pseudomonas putida RW10S1) which shows potent antibacterial activities against Gram-negative pathogens by inactivating succinate dehydrogenase. Based on the in-vivo studies, promysalin is hypothesized to be assembled from three building blocks: salicylic acid, proline, and myristic acid via a proposed but uncharacterized hybrid NRPS-PKS biosynthetic pathway. So far, no in-vitro biosynthetic studies have been reported for this promising antibiotic. Here, we report the first in-vitro reconstitution and biochemical characterization of two early enzymes on the pathway: PpgH, an isochorismate synthase (IS), and PpgG, an isochorismate pyruvate lyase (IPL) which are involved in the biosynthesis of salicylic acid, the polar fragment of promysalin. We also report a secondary chorismate mutase (CM) activity for PpgG. Based on our biochemical experiments, preliminary mechanistic proposals have been postulated for PpgH and PpgG. We believe this study will lay a strong foundation for elucidating the functions and mechanisms of other intriguing enzymes of the promysalin biosynthesis pathway, which may potentially unravel interesting enzyme chemistries and promote pathway engineering in the future.

Challenges in modelling dynamic processes in realistic nanostructured materials at operating conditions

The question is addressed in how far current modelling strategies are capable of modelling dynamic phenomena in realistic nanostructured materials at operating conditions. Nanostructured materials used in applications are far from perfect; they possess a broad range of heterogeneities in space and time extending over several orders of magnitude. Spatial heterogeneities from the subnanometre to the micrometre scale in crystal particles with a finite size and specific morphology, impact the material's dynamics. Furthermore, the material's functional behaviour is largely determined by the operating conditions. Currently, there exists a huge length-time scale gap between attainable theoretical length-time scales and experimentally relevant scales. Within this perspective, three key challenges are highlighted within the molecular modelling chain to bridge this length-time scale gap. Methods are needed that enable (i) building structural models for realistic crystal particles having mesoscale dimensions with isolated defects, correlated nanoregions, mesoporosity, internal and external surfaces; (ii) the evaluation of interatomic forces with quantum mechanical accuracy albeit at much lower computational cost than the currently used density functional theory methods and (iii) derivation of the kinetics of phenomena taking place in a multi-length-time scale window to obtain an overall view of the dynamics of the process. This article is part of a discussion meeting issue 'Supercomputing simulations of advanced materials'.

Best Practices on QM/MM Simulations of Biological Systems

During the second half of the 20th century, following structural biology hallmark works on DNA and proteins, biochemists shifted their questions from "what does this molecule look like?" to "how does this process work?". Prompted by the theoretical and practical developments in computational chemistry, this led to the emergence of biomolecular simulations and, along with the 2013 Nobel Prize in Chemistry, to the development of hybrid QM/MM methods. QM/MM methods are necessary whenever the problem we want to address involves chemical reactivity and/or a change in the system's electronic structure, with archetypal examples being the studies of an enzyme's reaction mechanism and a metalloprotein's active site. In the last decades QM/MM methods have seen an increasing adoption driven by their incorporation in widely used biomolecular simulation software. However, properly setting up a QM/MM simulation is not an easy task, and several issues need to be properly addressed to obtain meaningful results. In the present work, we describe both the theoretical concepts and practical issues that need to be considered when performing QM/MM simulations. We start with a brief historical perspective on the development of these methods and describe when and why QM/MM methods are mandatory. Then we show how to properly select and analyze the performance of the QM level of theory, the QM system size, and the position and type of the boundaries. We show the relevance of performing prior QM model system (or QM cluster) calculations in a vacuum and how to use the corresponding results to adequately calibrate those derived from QM/MM. We also discuss how to prepare the starting structure and how to select an adequate simulation strategy, including those based on geometry optimizations as well as free energy methods. In particular, we focus on the determination of free energy profiles using multiple steered molecular dynamics (MSMD) combined with Jarzynski's equation. Finally, we describe the results for two illustrative and complementary examples: the reaction performed by chorismate mutase and the study of ligand binding to hemoglobins. Overall, we provide many practical recommendations (or shortcuts) together with important conceptualizations that we hope will encourage more and more researchers to incorporate QM/MM studies into their research projects.

Modeling Ensembles of Enzyme Reaction Pathways with Hi-MSM Reveals the Importance of Accounting for Pathway Diversity

Conventional quantum mechanical-molecular mechanics (QM/MM) simulation approaches for modeling enzyme reactions often assume that there is one dominant reaction pathway and that this pathway can be sampled starting from an X-ray structure of the enzyme. These assumptions reduce computational cost; however, their validity has not been extensively tested. This is due in part to the lack of a rigorous formalism for integrating disparate pathway information from dynamical QM/MM calculations. Here, we present a way to model ensembles of reaction pathways efficiently using a divide-and-conquer strategy through Hierarchical Markov State Modeling (Hi-MSM). This approach allows information on multiple, distinct pathways to be incorporated into a chemical kinetic model, and it allows us to test these two assumptions. Applying Hi-MSM to the reaction carried out by dihydrofolate reductase (DHFR) we find (i) there are multiple, distinct pathways significantly contributing to the overall flux of the reaction that the conventional approach does not identify and (ii) that the conventional approach does not identify the dominant reaction pathway. Thus, both assumptions underpinning the conventional approach are violated. Since DHFR is a relatively small enzyme, and configuration space scales exponentially with protein size, accounting for multiple reaction pathways is likely to be necessary for most 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.

Multistructural microiteration combined with QM/MM-ONIOM electrostatic embedding

Multistructural microiteration (MSM) is a method to take account of contributions of multiple surrounding structures in a geometrical optimization or reaction path calculation using the quantum mechanics/molecular mechanics (QM/MM) ONIOM method. In this study, we combined MSM with the electrostatic embedding (EE) scheme of the QM/MM-ONIOM method by extending its original formulation for mechanical embedding (ME). MSM-EE takes account of the polarization in the QM region induced by point charges assigned to atoms in the multiple surrounding structures, where the point charges are scaled by the weight factor of each surrounding structure determined through MSM. The performance of MSM-EE was compared with that of the other methods, i.e., ONIOM-ME, ONIOM-EE, and MSM-ME, by applying them to three chemical processes: (1) chorismate-to-prephenate transformation in aqueous solution, (2) the same transformation as (1) in an enzyme, and (3) hydroxylation in p-hydroxybenzoate hydroxylase. These numerical tests of MSM-EE yielded barriers and reaction energies close to experimental values with computational costs comparable to those of the other three methods.

Factors That Determine the Variation of Equilibrium and Kinetic Properties of QM/MM Enzyme Simulations: QM Region, Conformation, and Boundary Condition

To analyze the impact of various technical details on the results of quantum mechanical (QM)/molecular mechanical (MM) enzyme simulations, including the QM region size, catechol-O-methyltransferase (COMT) is studied as a model system using an approximate QM/MM method (DFTB3/CHARMM). The results show that key equilibrium and kinetic properties for methyl transfer in COMT exhibit limited variations with respect to the size of the QM region, which ranges from ∼100 to ∼500 atoms in this study. With extensive sampling, local and global structural characteristics of the enzyme are largely conserved across the studied QM regions, while the nature of the transition state (e.g., secondary kinetic isotope effect) and reaction exergonicity are largely maintained. Deviations in the free energy profile with different QM region sizes are similar in magnitude to those observed with changes in other simulation protocols, such as different initial enzyme conformations and boundary conditions. Electronic structural properties, such as the covariance matrix of residual charge fluctuations, appear to exhibit rather long-range correlations, especially when the peptide backbone is included in the QM region; this observation holds when a range-separated DFT approach is used as the QM region, suggesting that delocalization error is unlikely the origin. Overall, the analyses suggest that multiple simulation details determine the results of QM/MM enzyme simulations with comparable contributions.

Wavefunction-Based Electrostatic-Embedding QM/MM Using CFOUR through MiMiC

We present an interface of the wavefunction-based quantum chemical software CFOUR to the multiscale modeling framework MiMiC. Electrostatic embedding of the quantum mechanical (QM) part is achieved by analytic evaluation of one-electron integrals in CFOUR, while the rest of the QM/molecular mechanical (MM) operations are treated according to the previous MiMiC-based QM/MM implementation. Long-range electrostatic interactions are treated by a multipole expansion of the potential from the QM electron density to reduce the computational cost without loss of accuracy. Testing on model water/water systems, we verified that the CFOUR interface to MiMiC is robust, guaranteeing fast convergence of the self-consistent field cycles and optimal conservation of the energy during the integration of the equations of motion. Finally, we verified that the CFOUR interface to MiMiC is compatible with the use of a QM/QM multiple time-step algorithm, which effectively reduces the cost of ab initio MD (AIMD) or QM/MM-MD simulations using higher level wavefunction-based approaches compared to cheaper density functional theory-based ones. The new wavefunction-based AIMD and QM/MM-MD implementations were tested and validated for a large number of wavefunction approaches, including Hartree-Fock and post-Hartree-Fock methods like Møller-Plesset, coupled-cluster, and complete active space self-consistent field.

High throughput and quantitative enzymology in the genomic era

Accurate predictions from models based on physical principles are the ultimate metric of our biophysical understanding. Although there has been stunning progress toward structure prediction, quantitative prediction of enzyme function has remained challenging. Realizing this goal will require large numbers of quantitative measurements of rate and binding constants and the use of these ground-truth data sets to guide the development and testing of these quantitative models. Ground truth data more closely linked to the underlying physical forces are also desired. Here, we describe technological advances that enable both types of ground truth measurements. These advances allow classic models to be tested, provide novel mechanistic insights, and place us on the path toward a predictive understanding of enzyme structure and function.

Multiscale Workflow for Modeling Ligand Complexes of Zinc Metalloproteins

Zinc metalloproteins are ubiquitous, with protein zinc centers of structural and functional importance, involved in interactions with ligands and substrates and often of pharmacological interest. Biomolecular simulations are increasingly prominent in investigations of protein structure, dynamics, ligand interactions, and catalysis, but zinc poses a particular challenge, in part because of its versatile, flexible coordination. A computational workflow generating reliable models of ligand complexes of biological zinc centers would find broad application. Here, we evaluate the ability of alternative treatments, using (nonbonded) molecular mechanics (MM) and quantum mechanics/molecular mechanics (QM/MM) at semiempirical (DFTB3) and density functional theory (DFT) levels of theory, to describe the zinc centers of ligand complexes of six metalloenzyme systems differing in coordination geometries, zinc stoichiometries (mono- and dinuclear), and the nature of interacting groups (specifically the presence of zinc-sulfur interactions). MM molecular dynamics (MD) simulations can overfavor octahedral geometries, introducing additional water molecules to the zinc coordination shell, but this can be rectified by subsequent semiempirical (DFTB3) QM/MM MD simulations. B3LYP/MM geometry optimization further improved the accuracy of the description of coordination distances, with the overall effectiveness of the approach depending upon factors, including the presence of zinc-sulfur interactions that are less well described by semiempirical methods. We describe a workflow comprising QM/MM MD using DFTB3 followed by QM/MM geometry optimization using DFT (e.g., B3LYP) that well describes our set of zinc metalloenzyme complexes and is likely to be suitable for creating accurate models of zinc protein complexes when structural information is more limited.

QM/MM and molecular dynamics investigation of the mechanism of covalent inhibition of TAK1 kinase

TAK1 is a serine/threonine kinase which is involved in the moderation of cell survival and death via the TNFα signalling pathway. It is also implicated in a range of cancer and anti-inflammatory diseases. Drug discovery efforts on this target have focused on both traditional reversible ATP-binding site inhibitors and increasingly popular irreversible covalent binding inhibitors. Irreversible inhibitors can offer benefits in terms of potency, selectivity and PK/PD meaning they are increasingly pursued where the strategy exists. TAK1 kinase differs from the better-known kinase EGFR in that the reactive cysteine nucleophile targeted by electrophilic inhibitors is located towards the back of the ATP binding site, not at its mouth. While a wealth of structural and computational effort has been spent exploring EGFR, only limited studies on TAK1 have been reported. In this work we report the first QM/MM study on TAK1 aiming to better understand aspects of covalent adduct formation. Our goal is to identify the general base in the catalytic reaction, whether the process proceeds via a stepwise or concerted pathway, and how the highly flexible G-loop and A-loop affect the catalytic cysteine located nearby.

Gold-Catalyzed Cycloisomerization of Sulfur Ylides to Dihydrobenzothiepines

The metal-promoted nucleophilic addition of sulfur ylides to π-systems is a well-established reactivity. However, the driving force of such transformations, elimination of a sulfide moiety, entails stoichiometric byproducts making them unfavorable in terms of atom economy. In this work, a new take on sulfur ylide chemistry is reported, an atom-economical gold(I)-catalyzed synthesis of dihydrobenzo[b]thiepines. The reaction proceeds under mild conditions at room temperature.

Multi-level free energy simulation with a staged transformation approach

Combining multiple levels of theory in free energy simulations to balance computational accuracy and efficiency is a promising approach for studying processes in the condensed phase. While the basic idea has been proposed and explored for quite some time, it remains challenging to achieve convergence for such multi-level free energy simulations as it requires a favorable distribution overlap between different levels of theory. Previous efforts focused on improving the distribution overlap by either altering the low-level of theory for the specific system of interest or ignoring certain degrees of freedom. Here, we propose an alternative strategy that first identifies the degrees of freedom that lead to gaps in the distributions of different levels of theory and then treats them separately with either constraints or restraints or by introducing an intermediate model that better connects the low and high levels of theory. As a result, the conversion from the low level to the high level model is done in a staged fashion that ensures a favorable distribution overlap along the way. Free energy components associated with different steps are mostly evaluated explicitly, and thus, the final result can be meaningfully compared to the rigorous free energy difference between the two levels of theory with limited and well-defined approximations. The additional free energy component calculations involve simulations at the low level of theory and therefore do not incur high computational costs. The approach is illustrated with two simple but non-trivial solution examples, and factors that dictate the reliability of the result are discussed.

Benchmarking quantum mechanical methods for calculating reaction energies of reactions catalyzed by enzymes

To assess the accuracy of different quantum mechanical methods for biochemical modeling, the reaction energies of 20 small model reactions (chosen to represent chemical steps catalyzed by commonly studied enzymes) were calculated. The methods tested included several popular Density Functional Theory (DFT) functionals, second-order Møller Plesset perturbation theory (MP2) and its spin-component scaled variant (SCS-MP2), and coupled cluster singles and doubles and perturbative triples (CCSD(T)). Different basis sets were tested. CCSD(T)/aug-cc-pVTZ results for all 20 reactions were used to benchmark the other methods. It was found that MP2 and SCS-MP2 reaction energy calculation results are similar in quality to CCSD(T) (mean absolute error (MAE) of 1.2 and 1.3 kcal mol−1, respectively). MP2 calculations gave a large error in one case, and are more subject to basis set effects, so in general SCS-MP2 calculations are a good choice when CCSD(T) calculations are not feasible. Results with different DFT functionals were of reasonably good quality (MAEs of 2.5–5.1 kcal mol−1), whereas popular semi-empirical methods (AM1, PM3, SCC-DFTB) gave much larger errors (MAEs of 11.6–14.6 kcal mol−1). These results should be useful in guiding methodological choices and assessing the accuracy of QM/MM calculations on enzyme-catalyzed reactions.

Temperature, Dynamics, and Enzyme-Catalyzed Reaction Rates

We review the adaptations of enzyme activity to different temperatures. Psychrophilic (cold-adapted) enzymes show significantly different activation parameters (lower activation enthalpies and entropies) from their mesophilic counterparts. Furthermore, there is increasing evidence that the temperature dependence of many enzyme-catalyzed reactions is more complex than is widely believed. Many enzymes show curvature in plots of activity versus temperature that is not accounted for by denaturation or unfolding. This is explained by macromolecular rate theory: A negative activation heat capacity for the rate-limiting chemical step leads directly to predictions of temperature optima; both entropy and enthalpy are temperature dependent. Fluctuations in the transition state ensemble are reduced compared to the ground state. We show how investigations combining experiment with molecular simulation are revealing fundamental details of enzyme thermoadaptation that are relevant for understanding aspects of enzyme evolution. Simulations can calculate relevant thermodynamic properties (such as activation enthalpies, entropies, and heat capacities) and reveal the molecular mechanisms underlying experimentally observed behavior.

The Molpro quantum chemistry package

Molpro is a general purpose quantum chemistry software package with a long development history. It was originally focused on accurate wavefunction calculations for small molecules but now has many additional distinctive capabilities that include, inter alia, local correlation approximations combined with explicit correlation, highly efficient implementations of single-reference correlation methods, robust and efficient multireference methods for large molecules, projection embedding, and anharmonic vibrational spectra. In addition to conventional input-file specification of calculations, Molpro calculations can now be specified and analyzed via a new graphical user interface and through a Python framework.

Design and SAR Analysis of Covalent Inhibitors Driven by Hybrid QM/MM Simulations

Quantum mechanics/molecular mechanics (QM/MM) hybrid technique is emerging as a reliable computational method to investigate and characterize chemical reactions occurring in enzymes. From a drug discovery perspective, a thorough understanding of enzyme catalysis appears pivotal to assist the design of inhibitors able to covalently bind one of the residues belonging to the enzyme catalytic machinery. Thanks to the current advances in computer power, and the availability of more efficient algorithms for QM-based simulations, the use of QM/MM methodology is becoming a viable option in the field of covalent inhibitor design. In the present review, we summarized our experience in the field of QM/MM simulations applied to drug design problems which involved the optimization of agents working on two well-known drug targets, namely fatty acid amide hydrolase (FAAH) and epidermal growth factor receptor (EGFR). In this context, QM/MM simulations gave valuable information in terms of geometry (i.e., of transition states and metastable intermediates) and reaction energetics that allowed to correctly predict inhibitor binding orientation and substituent effect on enzyme inhibition. What is more, enzyme reaction modelling with QM/MM provided insights that were translated into the synthesis of new covalent inhibitor featured by a unique combination of intrinsic reactivity, on-target activity, and selectivity.

Thermal activation of methane by MgO+: temperature dependent kinetics, reactive molecular dynamics simulations and statistical modeling

The kinetics methane activation (MgO⁺ + CH₄) was studied experimentally and computationally by running and analyzing reactive atomistic simulations.

Hydroxyl Radical-Coupled Electron-Transfer Mechanism of Flavin-Dependent Hydroxylases

Class A flavin-dependent hydroxylases (FdHs) catalyze the hydroxylation of organic compounds in a site- and stereoselective manner. In stark contrast, conventional synthetic routes require environmentally hazardous reagents and give modest yields. Thus, understanding the detailed mechanism of this class of enzymes is essential to their rational manipulation for applications in green chemistry and pharmaceutical production. Both electrophilic substitution and radical intermediate mechanisms have been proposed as interpretations of FdH hydroxylation rates and optical spectra. While radical mechanistic steps are often difficult to examine directly, modern quantum chemistry calculations combined with statistical mechanical approaches can yield detailed mechanistic models providing insights that can be used to differentiate reaction pathways. In the current work, we report quantum mechanical/molecular mechanical (QM/MM) calculations on the fungal TropB enzyme that shows an alternative reaction pathway in which hydroxylation through a hydroxyl radical-coupled electron-transfer mechanism is significantly favored over electrophilic substitution. Furthermore, QM/MM calculations on several modified flavins provide a more consistent interpretation of the experimental trends in the reaction rates seen experimentally for a related enzyme, para-hydroxybenzoate hydroxylase. These calculations should guide future enzyme and substrate design strategies and broaden the scope of biological spin chemistry.

The reaction mechanism of Zika virus NS2B/NS3 serine protease inhibition by dipeptidyl aldehyde: a QM/MM study

Zika virus (ZIKV) infection has become a global public health problem, associated with microcephaly in newborns and Guillain-Barré syndrome in adults. Currently, there are no commercially available anti-ZIKV drugs. The viral protease NS2B/NS3, which is involved in viral replication and maturation, is a potential drug target. Peptidomimetic aldehyde inhibitors bind covalently to the catalytic S135 of the NS3 protease. Here, we apply hybrid quantum mechanics/molecular mechanics (QM/MM) free-energy simulations at the PDDG-PM3/ff14SB level to investigate the inhibition mechanism of the ZIKV protease by a dipeptidyl aldehyde inhibitor (acyl-KR-aldehyde). The results show that proton transfer from the catalytic S135 to H51 occurs in concert with nucleophilic addition on the aldehyde warhead by S135. The anionic covalent complex between the dipeptidyl aldehyde and the ZIKV protease is analogous to the tetrahedral intermediate for substrate hydrolysis. Spontaneous protonation by H51 forms the hemiacetal. In addition, we use correlated ab initio QM/MM potential energy path calculations at levels up to LCCSD(T)/(aug)-cc-pVTZ to obtain accurate potential energy profiles of the reaction, which also support a concerted mechanism. These results provide detailed insight into the mechanism of ZIKV protease inhibition by a peptidyl aldehyde inhibitor, which will guide in the design of inhibitors.

Projector-Based Embedding Eliminates Density Functional Dependence for QM/MM Calculations of Reactions in Enzymes and Solution

Combined quantum mechanics/molecular mechanics (QM/MM) methods are increasingly widely utilized in studies of reactions in enzymes and other large systems. Here, we apply a range of QM/MM methods to investigate the Claisen rearrangement of chorismate to prephenate, in solution, and in the enzyme chorismate mutase. Using projector-based embedding in a QM/MM framework, we apply treatments up to the CCSD(T) level. We test a range of density functional QM/MM methods and QM region sizes. The results show that the calculated reaction energetics are significantly more sensitive to the choice of density functional than they are to the size of the QM region in these systems. Projector-based embedding of a wave function method in DFT reduced the 13 kcal/mol spread in barrier heights calculated at the DFT/MM level to a spread of just 0.3 kcal/mol, essentially eliminating dependence on the functional. Projector-based embedding of correlated ab initio methods provides a practical method for achieving high accuracy for energy profiles derived from DFT and DFT/MM calculations for reactions in condensed phases.

Comparing Hydrolysis and Transglycosylation Reactions Catalyzed by Thermus thermophilus β-Glycosidase. A Combined MD and QM/MM Study

The synthesis of oligosaccharides and other carbohydrate derivatives is of relevance for the advancement of glycosciences both at the fundamental and applied level. For many years, glycosyl hydrolases (GHs) have been explored to catalyze the synthesis of glycosidic bonds. In particular, retaining GHs can catalyze a transglycosylation (T) reaction that competes with hydrolysis (H). This has been done either employing controlled conditions in wild type GHs or by engineering new mutants. The goal, which is to increase the T/H ratio, has been achieved with moderate success in several cases despite the fact that the molecular basis for T/H modulation are unclear. Here we have used QM(DFT)/MM calculations to compare the glycosylation, hydrolysis and transglycosylation steps catalyzed by wild type Thermus thermophilus β-glycosidase (family GH1), a retaining glycosyl hydrolase for which a transglycosylation yield of 36% has been determined experimentally. The three transition states have a strong oxocarbenium character and ring conformations between ⁴H₃ and ⁴E. The atomic charges at the transition states for hydrolysis and transglycosylation are very similar, except for the more negative charge of the oxygen atom of water when compared to that of the acceptor Glc. The glycosylation transition state has a stronger S_N2 character than the deglycosylation ones and the proton transfer is less advanced. At the QM(PBE0/TZVP)/MM level, the TS for transglycosylation has shorter O4_GLC-C1_FUC (forming bond) distance and longer OE2_GLU338-C1_FUC (breaking) distance than the hydrolysis one, although the H_ACC proton is closer to the Glu164 base in the hydrolysis TS. The QM(SCC-DFTB)/MM free energy maxima show the inverted situation, although the hydrolysis TS presents significant structural fluctuations. The 3-OH_GLC group of the acceptor Glc (transglycosylation) and WAT432 (neighbor water in hydrolysis) are identified to stabilize the oxocarbenium transition states through interaction with O5_FUC and O4_FUC. The analysis of interaction suggests that perturbing the Glu392-Fuc interaction could increase the T/H ratio, either by direct mutation of this residue or indirectly as reported experimentally in the Asn390I and Phe401S cases. The molecular understanding of similarities and differences between hydrolysis and transglycosylation steps may be of help in the design of new biocatalysts for glycan synthesis.

Quantum Mechanics/Molecular Mechanics (QM/MM) Calculations Support a Concerted Reaction Mechanism for the Zika Virus NS2B/NS3 Serine Protease with Its Substrate

Zika virus (ZIKV) is mainly transmitted to humans by Aedes species mosquitoes and is associated with serious pathological disorders including microcephaly in newborns and Guillain-Barré syndrome in adults. Currently, there is no vaccine or anti-ZIKV drug available for preventing or controlling ZIKV infection. An attractive drug target for ZIKV treatment is a two-compartment (NS2B/NS3) serine protease that processes viral polyprotein during infection. Here, conventional molecular dynamics simulations of the ZIKV protease in complex with peptide substrate (TGKRS) sequence at the C-terminus of NS2B show that the substrate is in the active conformation for the cleavage reaction by ZIKV protease. Hybrid quantum mechanics/molecular mechanics (QM/MM) umbrella sampling simulations (PM6/ff14SB) of acylation results reveal that proton transfer from S135 to H51 and nucleophilic attack on the substrate by S135 are concerted. The rate-limiting step involves the formation of a tetrahedral intermediate. In addition, the single-point energy QM/MM calculations, precisely at the level of coupled cluster theory (LCCSD(T)/(aug)-cc-pVTZ), were performed to correct the potential energy profiles for the first step of the acylation process. The average computed activation barrier at this level of theory is 16.3 kcal mol^-1. Therefore, the computational approaches presented here are helpful for further designing of NS2B/NS3 inhibitors based on transition-state analogues.

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.

Theoretical Evaluation of the Reaction Mechanism of Serine Hydroxymethyltransferase

Serine hydroxymethyltransferase (SHMT) is a pyridoxal phosphate (PLP)-dependent enzyme that catalyzes the reversible conversion of serine and tetrahydrofolate (THF) to glycine and 5,10-methylene THF. SHMT is a folate pathway enzyme and is therefore of considerable medical interest due to its role as an important intervention point for antimalarial, anticancer, and antibacterial treatments. Despite considerable experimental effort, the precise reaction mechanism of SHMT remains unclear. In this study, we explore the mechanism of SHMT to determine the roles of active site residues and the nature and the sequence of chemical steps. Molecular dynamics (MD) methods were employed to generate a suitable starting structure which then underwent analysis using hybrid quantum mechanical/molecular mechanical (QM/MM) simulations. The QM region consisted of 12 key residues, two substrates, and explicit solvent. Our results show that the catalytic reaction proceeds according to a retro-aldol synthetic process with His129 acting as the general base in the reaction. The rate-determining step involves the cleavage of the PLP-serine aldimine C_α-C_β bond and the formation of formaldehyde in line with experimental evidence. The pyridyl ring of the PLP-serine aldimine substrate exists in deprotonated form, being stabilized directly by Asp208 via a strong H-bond, as well as through interactions with Arg371, Lys237, and His211, and with the surrounding protein which was electrostatically embedded. This knowledge has the potential to impact the design and development of new inhibitors.

Oxidative dehalogenation and denitration by a flavin-dependent monooxygenase is controlled by substrate deprotonation

Enzymes that are capable of detoxifying halogenated phenols (HPs) and nitrophenols (NPs) are valuable for bioremediation and waste biorefining. HadA monooxygenase was found to perform dual functions of oxidative dehalogenation (hydroxylation plus halide elimination) and denitration (hydroxylation plus nitro elimination). Rate constants associated with individual steps of HadA reactions with phenol, halogenated phenols and nitrophenols were measured using combined transient kinetic approaches of stopped-flow absorbance/fluorescence and rapid-quench flow techniques. Density functional theory was used to calculate the thermodynamic and electronic parameters associated with hydroxylation and group elimination steps. These parameters were correlated with the rate constants of hydroxylation, group elimination, and overall product formation to identify factors controlling individual steps. The results indicated that the hydroxylation rate constant is higher when the pK a of the phenolic group is lower, i.e. it is more easily deprotonated, but not higher when the energy gap between the E LUMO of the C4a-hydroperoxy-FAD intermediate and the E HOMO of the phenolate substrate is lower. These data suggest that the substrate deprotonation has a higher energy barrier than the -OH transfer, and thus controls the hydroxylation step. For the group elimination, the process is controlled by the ability of the C-X bond to break. For the overall product formation (hydroxylation and group elimination combined), this analysis showed that the rate constant of product formation is dependent on the pK a value of the substrate, indicating that the overall reaction is controlled by substrate deprotonation. This step also likely has the highest energy barrier and thus controls the overall process of oxidative dehalogenation and denitration by HadA. This report is the first to identify a key mechanistic factor controlling the enzymatic processes of oxidative dehalogenation and denitration.

Error-Controlled Exploration of Chemical Reaction Networks with Gaussian Processes

For a theoretical understanding of the reactivity of complex chemical systems, relative energies of stationary points on potential energy hypersurfaces need to be calculated to high accuracy. Due to the large number of intermediates present in all but the simplest chemical processes, approximate quantum chemical methods are required that allow for fast evaluations of the relative energies but at the expense of accuracy. Despite the plethora of benchmark studies, the accuracy of a quantum chemical method is often difficult to assess. Moreover, a significant improvement of a method's accuracy (e.g., through reparameterization or systematic model extension) is rarely possible. Here, we present a new approach that allows for the systematic, problem-oriented, and rolling improvement of quantum chemical results through the application of Gaussian processes. Due to its Bayesian nature, reliable error estimates are provided for each prediction. A reference method of high accuracy can be employed if the uncertainty associated with a particular calculation is above a given threshold. The new data point is then added to a growing data set in order to continuously improve the model and, as a result, all subsequent predictions. Previous predictions are validated by the updated model to ensure that uncertainties remain within the given confidence bound, which we call backtracking. We demonstrate our approach with the example of a complex chemical reaction network.

Analysis of Phosphoryl-Transfer Enzymes with QM/MM Free Energy Simulations

We discuss the application of quantum mechanics/molecular mechanics (QM/MM) free energy simulations to the analysis of phosphoryl transfers catalyzed by two enzymes: alkaline phosphatase and myosin. We focus on the nature of the transition state and the issue of mechanochemical coupling, respectively, in the two enzymes. The results illustrate unique insights that emerged from the QM/MM simulations, especially concerning the interpretation of experimental data regarding the nature of enzymatic transition states and coupling between global structural transition and catalysis in the active site. We also highlight a number of technical issues worthy of attention when applying QM/MM free energy simulations, and comment on a number of technical and mechanistic issues that require further studies.

Rationalising pKa shifts in Bacillus circulans xylanase with computational studies

Bacillus circulans xylanase (BcX), a family 11 glycoside hydrolase, catalyses the hydrolysis of xylose polymers with a net retention of stereochemistry. Glu172 in BcX is believed to act as a general acid by protonating the aglycone during glycosylation, and then as a general base to facilitate the deglycosylation step. The key to the dual role of this general acid/base lies in its protonation states, which depend on its intrinsic pK_a value and the specific environment which it resides within. To fully understand the detailed molecular features in BcX to establish the dual role of Glu172, we present a combined study based on both atomistic simulations and empirical models to calculate pK_a shifts for the general acid/base Glu172 in BcX at different functional states. Its pK_a values and those of nearby residues, obtained based on QM/MM free energy calculations, MCCE and PROPKA, show a good agreement with available experimental data. Additionally, our study provides additional insights into the effects of structural and electrostatic perturbations caused by mutations and chemical modifications, suggesting that the local solvation environment and mutagenesis of the residues adjacent to Glu172 establish its dual role during hydrolysis. The strengths and limitations of various methods for calculating pK_as and pK_a shifts have also been discussed.

Multiscale analysis of enantioselectivity in enzyme-catalysed 'lethal synthesis' using projector-based embedding

The action of fluoroacetate as a broad-spectrum mammalian pesticide depends on the 'lethal synthesis' of fluorocitrate by citrate synthase, through a subtle enantioselective enolization of fluoroacetyl-coenzyme A. In this work, we demonstrate how a projection-based embedding method can be applied to calculate coupled cluster (CCSD(T)) reaction profiles from quantum mechanics/molecular mechanics optimized pathways for this enzyme reaction. Comparison of pro-R and pro-S proton abstraction in citrate synthase at the CCSD(T)-in-DFT//MM level yields the correct enantioselectivity. We thus demonstrate the potential of projection-based embedding for determining stereoselectivity in enzymatic systems. We further show that the method is simple to apply, eliminates variability due to the choice of density functional theory functional and allows the efficient calculation of CCSD(T) quality enzyme reaction barriers.

Ab Initio QM/MM Modeling of the Rate-Limiting Proton Transfer Step in the Deamination of Tryptamine by Aromatic Amine Dehydrogenase

Aromatic amine dehydrogenase (AADH) and related enzymes are at the heart of debates on the roles of quantum tunneling and protein dynamics in catalysis. The reaction of tryptamine in AADH involves significant quantum tunneling in the rate-limiting proton transfer step, shown by large H/D primary kinetic isotope effects (KIEs), with unusual temperature dependence. We apply correlated ab initio combined quantum mechanics/molecular mechanics (QM/MM) methods, at levels up to local coupled cluster theory (LCCSD(T)/(aug)-cc-pVTZ), to calculate accurate potential energy surfaces for this reaction, which are necessary for quantitative analysis of tunneling contributions and reaction dynamics. Different levels of QM/MM treatment are tested. Multiple pathways are calculated with fully flexible transition state optimization by the climbing-image nudged elastic band method at the density functional QM/MM level. The average LCCSD(T) potential energy barriers to proton transfer are 16.7 and 14.0 kcal/mol for proton transfer to the two carboxylate atoms of the catalytic base, Asp128β. The results show that two similar, but distinct pathways are energetically accessible. These two pathways have different barriers, exothermicity and curvature, and should be considered in analyses of the temperature dependence of reaction and KIEs in AADH and other enzymes. These results provide a benchmark for this prototypical enzyme reaction and will be useful for developing empirical models, and analyzing experimental data, to distinguish between different conceptual models of enzyme catalysis.

New Insights into the Catalytic Mechanism of Aldose Reductase: A QM/MM Study

Aldose reductase is the first enzyme of the polyol pathway in which glucose is converted to fructose via sorbitol. The understanding of this key enzyme is important as it has been linked to some diabetes mellitus complications. The mechanism of the enzyme was investigated using a hybrid quantum mechanics/molecular mechanics (QM/MM) method. It was found that depending on the protonation state of His110 the mechanism can be concerted or stepwise and the proton donor can be either Tyr48 or His110. These findings are different from the previous theoretical studies based on QM/MM calculations using either AM1 or HF/4-31G, in which the reduction is, respectively, a stepwise or one-step process. The QM/MM energy barriers for the reduction of d-glyceraldehyde were evaluated at a B3LYP/6-31G* level for both HIP and HIE protonation states of His110. These were, respectively, 6.5 ± 2.2 and 16.7 ± 1.0 kcal/mol, which makes only the HIE protonation state consistent with the experimental value of 14.8 kcal/mol derived from kinetics experiments and makes Tyr48 the most probable proton donor.

Unraveling the Addition-Elimination Mechanism of EPSP Synthase through Computer Modeling

Enolpyruvyl transfer from phosphoenolpyruvate (PEP) to the hydroxyl group of shikimate-5-OH-3-phosphate (S3P) is catalyzed by 5-enolpyruvylshikimate 3-phosphate (EPSP) synthase in a reaction that involves breaking the C-O bond of PEP. Catalysis involves an addition-elimination mechanism with the formation of a tetrahedral intermediate (THI). Experiments have elucidated the mechanism of THI formation and breakdown. However, the catalytic action of EPSP synthase and the individual roles of catalytic residues Asp313 and Glu341 remains unclear. We have used a hybrid quantum mechanical/molecular mechanical (QM/MM) approach to explore the free energy surface in a reaction catalyzed by EPSP synthase. The Glu341 was the most favorable acid/base catalyst. Our results indicate that the protonation of PEP C3 precedes the nucleophilic attack on PEP C2 in the addition mechanism. Also, the breaking of the C-O bond of THI to form an EPSP cation intermediate must occur before proton transfer from PEP C3 to Glu341 in the elimination mechanism. Analysis of the FES supports cationic intermediate formation during the reaction catalyzed by EPSP synthase. Finally, the computational model indicates a proton transfer shift (Hammond shift) from Glu341 to C3 for an enzyme-based reaction with the shifted transition state, earlier than in the reference reaction in water.

Quantum chemical modeling of the reaction path of chorismate mutase based on the experimental substrate/product complex

Chorismate mutase is a well-known model enzyme, catalyzing the Claisen rearrangement of chorismate to prephenate. Recent high-resolution crystal structures along the reaction coordinate of this enzyme enabled computational analyses at unprecedented detail. Using quantum chemical simulations, we investigated how the catalytic reaction mechanism is affected by electrostatic and hydrogen-bond interactions. Our calculations showed that the transition state (TS) was mainly stabilized electrostatically, with Arg90 playing the leading role. The effect was augmented by selective hydrogen-bond formation to the TS in the wild-type enzyme, facilitated by a small-scale local induced fit. We further identified a previously underappreciated water molecule, which separates the negative charges during the reaction. The analysis includes the wild-type enzyme and a non-natural enzyme variant, where the catalytic arginine was replaced with an isosteric citrulline residue.