Do Quantum Mechanical Energies Calculated for Small Models of Protein-Active Sites Converge?

Quick-and-Easy Validation of Protein-Ligand Binding Models Using Fragment-Based Semiempirical Quantum Chemistry

Electronic structure calculations in enzymes converge very slowly with respect to the size of the model region that is described using quantum mechanics (QM), requiring hundreds of atoms to obtain converged results and exhibiting substantial sensitivity (at least in smaller models) to which amino acids are included in the QM region. As such, there is considerable interest in developing automated procedures to construct a QM model region based on well-defined criteria. However, testing such procedures is burdensome due to the cost of large-scale electronic structure calculations. Here, we show that semiempirical methods can be used as alternatives to density functional theory (DFT) to assess convergence in sequences of models generated by various automated protocols. The cost of these convergence tests is reduced even further by means of a many-body expansion. We use this approach to examine convergence (with respect to model size) of protein-ligand binding energies. Fragment-based semiempirical calculations afford well-converged interaction energies in a tiny fraction of the cost required for DFT calculations. Two-body interactions between the ligand and single-residue amino acid fragments afford a low-cost way to construct a "QM-informed" enzyme model of reduced size, furnishing an automatable active-site model-building procedure. This provides a streamlined, user-friendly approach for constructing ligand binding-site models that requires neither a priori information nor manual adjustments. Extension to model-building for thermochemical calculations should be straightforward.

Eliminating Imaginary Vibrational Frequencies in Quantum-Chemical Cluster Models of Enzymatic Active Sites

In constructing finite models of enzyme active sites for quantum-chemical calculations, atoms at the periphery of the model must be constrained to prevent unphysical rearrangements during geometry relaxation. A simple fixed-atom or "coordinate-lock" approach is commonly employed but leads to undesirable artifacts in the form of small imaginary frequencies. These preclude evaluation of finite-temperature free-energy corrections, limiting thermochemical calculations to enthalpies only. Full-dimensional vibrational frequency calculations are possible by replacing the fixed-atom constraints with harmonic confining potentials. Here, we compare that approach to an alternative strategy in which fixed-atom contributions to the Hessian are simply omitted. While the latter strategy does eliminate imaginary frequencies, it tends to underestimate both the zero-point energy and the vibrational entropy while introducing artificial rigidity. Harmonic confining potentials eliminate imaginary frequencies and provide a flexible means to construct active-site models that can be used in unconstrained geometry relaxations, affording better convergence of reaction energies and barrier heights with respect to the model size, as compared to models with fixed-atom constraints.

The influence of model building schemes and molecular dynamics sampling on QM-cluster models: the chorismate mutase case study

Most QM-cluster models of enzymes are constructed based on X-ray crystal structures, which limits comparison to in vivo structure and mechanism. The active site of chorismate mutase from Bacillus subtilis and the enzymatic transformation of chorismate to prephenate is used as a case study to guide construction of QM-cluster models built first from the X-ray crystal structure, then from molecular dynamics (MD) simulation snapshots. The Residue Interaction Network ResidUe Selector (RINRUS) software toolkit, developed by our group to simplify and automate the construction of QM-cluster models, is expanded to handle MD to QM-cluster model workflows. Several options, some employing novel topological clustering from residue interaction network (RIN) information, are evaluated for generating conformational clustering from MD simulation. RINRUS then generates a statistical thermodynamic framework for QM-cluster modeling of the chorismate mutase mechanism via refining 250 MD frames with density functional theory (DFT). The 250 QM-cluster models sampled provide a mean ΔG‡ of 10.3 ± 2.6 kcal mol-1 compared to the experimental value of 15.4 kcal mol-1 at 25 °C. While the difference between theory and experiment is consequential, the level of theory used is modest and therefore "chemical" accuracy is unexpected. More important are the comparisons made between QM-cluster models designed from the X-ray crystal structure versus those from MD frames. The large variations in kinetic and thermodynamic properties arise from geometric changes in the ensemble of QM-cluster models, rather from the composition of the QM-cluster models or from the active site-solvent interface. The findings open the way for further quantitative and reproducible calibration in the field of computational enzymology using the model construction framework afforded with the RINRUS software toolkit.

Mechanism of AGT-Mediated Repair of dG-dC Cross-Links in the Drug Resistance to Chloroethylnitrosoureas: Molecular Docking, MD Simulation, and ONIOM (QM/MM) Investigation

Chloroethylnitrosoureas (CENUs) are important chemotherapies applied in the treatment of cancer. They exert anticancer activity by inducing DNA interstrand cross-links (ICLs) via the formation of two O6-alkylguanine intermediates, O6-chloroethylguanine (O6-ClEtG) and N1,O6-ethanoguanine (N1,O6-EtG). However, O6-alkylguanine-DNA alkyltransferase (AGT), a DNA-repair enzyme, can restore the O6-alkylguanine damages and thereby obstruct the formation of ICLs (dG-dC cross-link). In this study, the inhibitory mechanism of ICL formation was investigated to elucidate the drug resistance of CENUs mediated by AGT in detail. Based on the structures of the substrate-enzyme complexes obtained from docking and MD simulations, two ONIOM (QM/MM) models with different sizes of the QM region were constructed. The model with a larger QM region, which included the substrate (O6-ClEtG or N1,O6-EtG), a water molecule, and five residues (Tyr114, Cys145, His146, Lys165, and Glu172) in the active pocket of AGT, accurately described the repairing reaction and generated the results coinciding with the experimental outcomes. The repair process consists of two sequential steps: hydrogen transfer to form a thiolate anion on Cys145 and alkyl transfer from the O6 site of guanine (the rate-limiting step). The repair of N1,O6-EtG was more favorable than that of O6-ClEtG from both kinetics and thermodynamics aspects. Moreover, the comparison of the repairing process with the formation of dG-dC cross-link and the inhibition of AGT by O6-benzylguanine (O6-BG) showed that the presence of AGT could effectively interrupt the formation of ICLs leading to drug resistance, and the inhibition of AGT by O6-BG that was energetically more favorable than the repair of O6-ClEtG could not prevent the repair of N1,O6-EtG. Therefore, it is necessary to completely eliminate AGT activity before CENUs medication to enhance the chemotherapeutic effectiveness. This work provides reasonable explanations for the supposed mechanism of AGT-mediated drug resistance of CENUs and will assist in the development of novel CENU chemotherapies and their medication strategies.

Dynamic Protonation States Underlie Carbene Formation in ThDP-Dependent Enzymes: A Theoretical Study

The activation mechanism of thiamine diphosphate (ThDP) in enzymes has long been the subject of intense research and controversial discussion. Particularly contentious is the formation of a carbene intermediate, the first one observed in an enzyme. For the formation of the carbene to take place, both intramolecular and intermolecular proton transfer pathways have been proposed. However, the physiologically relevant pH of ThDP-dependent enzymes around neutrality does not seem to be suitable for the formation of such reactive chemical species. Herein, we investigate the general mechanism of activation of the ThDP cofactor in human transketolase (TKT), by means of electronic structure methods. We show that in the case of the human TKT, the carbene species is accessible through a pKa shift induced by the electrostatics of a neighboring histidine residue (H110), whose protonation state change modulates the pKa of ThDP and suppresses the latter by more than 6 pH units. Our findings highlight that ThDP enzymes activate the cofactor beyond simple geometric constraints and the canonical glutamate. Such observations in nature can pave the way for the design of biomimetic carbene catalysts and the engineering of tailored enzymatic carbenes.

Environmental biotransformation mechanisms by flavin-dependent monooxygenase: A computational study

The enzyme-catalyzed metabolic biotransformation of xenobiotics plays a significant role in toxicology evolution and subsequently environmental health risk assessment. Recent studies noted that the phase I human flavin-dependent monooxygenase (e.g., FMO3) can catalyze xenobiotics into more toxic metabolites. However, details of the metabolic mechanisms are insufficient. To fill the mechanism in the gaps, the systemic density functional theory calculations were performed to elucidate diverse FMO-catalyzed oxidation reactions toward environmental pollutants, including denitrification (e.g., nitrophenol), N-oxidation (e.g., nicotine), desulfurization (e.g., fonofos), and dehalogenation (e.g., pentachlorophenol). Similar to the active center compound 0 of cytochrome P450, FMO mainly catalyzed reactions with the structure of the tricyclic isoalloxazine C-4a-hydroperoxide (FADHOOH). As will be shown, FMO-catalyzed pathways are more favorable with a concerted than stepwise mechanism; Deprotonation is necessary to initiate the oxidation reactions for phenolic substrates; The regioselectivity of nicotine by FMO prefers the N-oxidation other than N-demethylation pathway; Formation of the P-S-O triangle ring is the key step for desulfurization of fonofos by FMO. We envision that these fundamental mechanisms catalyzed by FMO with a computational method can be extended to other xenobiotics of similar structures, which may aid the high-throughput screening and provide theoretical predictions in the future.

Catalytic Reaction Mechanism of Glyoxalase II: A Quantum Mechanics/Molecular Mechanics Study

Methylglyoxal (MG) is a reactive and toxic compound produced in carbohydrate, lipid, and amino acid metabolism. The glyoxalase system is the main detoxifying route for MG and consists of two enzymes, glyoxalase I (GlxI) and glyoxalase II (GlxII). GlxI catalyzes the formation of S-d-lactoylglutathione from hemithioacetal, and GlxII converts this intermediate to d-lactate. A relationship between the glyoxalase system and some diseases like diabetes has been shown, and inhibiting enzymes of this system may be an effective means of controlling certain diseases. A detailed understanding of the reaction mechanism of an enzyme is essential to the rational design of competitive inhibitors. In this work, we use quantum mechanics/molecular mechanics (QM/MM) calculations and energy refinement utilizing the big-QM and QM/MM thermodynamic cycle perturbation methods to propose a mechanism for the GlxII reaction that starts with a nucleophilic attack of the bridging OH^- group on the substrate. The coordination of the substrate to the Zn ions places its electrophilic center close to the hydroxide group, enabling the reaction to proceed. Our estimated reaction energies are in excellent agreement with experimental data, thus demonstrating the reliability of our approach and the proposed mechanism. Additionally, we examined alternative protonation states of Asp-29, Asp-58, Asp-134, and the bridging hydroxide ion in the catalytic process. However, these give less favorable reactions, a poorer reproduction of the crystal structure geometry of the active site, and higher root-mean-squared deviations of the active site residues in molecular dynamics simulations.

Fragment-Based Calculations of Enzymatic Thermochemistry Require Dielectric Boundary Conditions

Electronic structure calculations on enzymes require hundreds of atoms to obtain converged results, but fragment-based approximations offer a cost-effective solution. We present calculations on enzyme models containing 500-600 atoms using the many-body expansion, comparing to benchmarks in which the entire enzyme-substrate complex is described at the same level of density functional theory. When the amino acid fragments contain ionic side chains, the many-body expansion oscillates under vacuum boundary conditions but rapid convergence is restored using low-dielectric boundary conditions. This implies that full-system calculations in the gas phase are inappropriate benchmarks for assessing errors in fragment-based approximations. A three-body protocol retains sub-kilocalorie per mole fidelity with respect to a supersystem calculation, as does a two-body calculation combined with a full-system correction at a low-cost level of theory. These protocols pave the way for application of high-level quantum chemistry to large systems via rigorous, ab initio treatment of many-body polarization.

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.

Harder, better, faster, stronger: Large-scale QM and QM/MM for predictive modeling in enzymes and proteins

Computational prediction of enzyme mechanism and protein function requires accurate physics-based models and suitable sampling. We discuss recent advances in large-scale quantum mechanical (QM) modeling of biochemical systems that have reduced the cost of high-accuracy models. Tradeoffs between sampling and accuracy have motivated modeling with molecular mechanics (MM) in a multiscale QM/MM or iterative approach. Limitations to both conventional density-functional theory and classical MM force fields remain for describing noncovalent interactions in comparison to experiment or wavefunction theory. Because predictions of enzyme action (i.e. electrostatics), free energy barriers, and mechanisms are sensitive to the protocol and embedding method in QM/MM, convergence tests and systematic methods for quantifying QM-level interactions are a needed, active area of development.

Cheminformatic quantum mechanical enzyme model design: A catechol-O-methyltransferase case study

To accurately simulate the inner workings of an enzyme active site with quantum mechanics (QM), not only must the reactive species be included in the model but also important surrounding residues, solvent, or coenzymes involved in crafting the microenvironment. Our lab has been developing the Residue Interaction Network Residue Selector (RINRUS) toolkit to utilize interatomic contact network information for automated, rational residue selection and QM-cluster model generation. Starting from an x-ray crystal structure of catechol-O-methyltransferase, RINRUS was used to construct a series of QM-cluster models. The reactant, product, and transition state of the methyl transfer reaction were computed for a total of 550 models, and the resulting free energies of activation and reaction were used to evaluate model convergence. RINRUS-designed models with only 200-300 atoms are shown to converge. RINRUS will serve as a cornerstone for improved and automated cheminformatics-based enzyme model design.

Thermochemical and Quantum Descriptor Calculations for Gaining Insight into Ricin Toxin A (RTA) Inhibitors

In this work, we performed a study to assess the interactions between the ricin toxin A (RTA) subunit of ricin and some of its inhibitors using modern semiempirical quantum chemistry and ONIOM quantum mechanics/molecular mechanics (QM/MM) methods. Two approaches were followed (calculation of binding enthalpies, ΔH _bind, and reactivity quantum chemical descriptors) and compared with the respective half-maximal inhibitory concentration (IC₅₀) experimental data, to gain insight into RTA inhibitors and verify which quantum chemical method would better describe RTA-ligand interactions. The geometries for all RTA-ligand complexes were obtained after running classical molecular dynamics simulations in aqueous media. We found that single-point energy calculations of ΔH _bind with the PM6-DH+, PM6-D3H4, and PM7 semiempirical methods and ONIOM QM/MM presented a good correlation with the IC₅₀ data. We also observed, however, that the correlation decreased significantly when we calculated ΔH _bind after full-atom geometry optimization with all semiempirical methods. Based on the results from reactivity descriptors calculations for the cases studied, we noted that both types of interactions, molecular overlap and electrostatic interactions, play significant roles in the overall affinity of these ligands for the RTA binding pocket.

QM/MM Study of the Catalytic Reaction of Myrosinase; Importance of Assigning Proper Protonation States of Active-Site Residues

Myrosinase from Sinapis alba hydrolyzes glycosidic bonds of β-d-S-glucosides. The enzyme shows an enhanced activity in the presence of l-ascorbic acid. In this work, we employed combined quantum mechanical and molecular mechanical (QM/MM) calculations and molecular dynamics simulations to study the catalytic reaction of wild-type myrosinase and its E464A, Q187A, and Q187E mutants. Test calculations show that a proper QM region to study the myrosinase reaction must contain the whole substrate, models of Gln-187, Glu-409, Gln-39, His-141, Asn-186, Tyr-330, Glu-464, Arg-259, and a water molecule. Furthermore, to make the deglycosylation step possible, Arg-259 must be charged, Glu-464 must be protonated on OE2, and His-141 must be protonated on the NE2 atom. The results indicate that assigning proper protonation states of the residues is more important than the size of the model QM system. Our model reproduces the anomeric retaining characteristic of myrosinase and also reproduces the experimental fact that ascorbate increases the rate of the reaction. A water molecule in the active site, positioned by Gln-187, helps the aglycon moiety of the substrate to stabilize the buildup of negative charge during the glycosylation reaction and this in turn makes the moiety a better leaving group. The water molecule also lowers the glycosylation barrier by ∼9 kcal/mol. The results indicate that the Q187E and E464A mutants but not the Q187A mutant can perform the glycosylation step. However, the energy profiles for the deglycosylation step of the mutants are not similar to that of the wild-type enzyme. The Glu-464 residue lowers the barriers of the glycosylation and deglycosylation steps. The ascorbate ion can act as a general base in the reaction of the wild-type enzyme only if the Glu-464 and His-141 residues are properly protonated.

QM/MM study of the binding of H2 to MoCu CO dehydrogenase: development and applications of improved H2 van der Waals parameters

The MoCu CO dehydrogenase enzyme not only transforms CO into CO₂ but it can also oxidise H₂. Even if its hydrogenase activity has been known for decades, a debate is ongoing on the most plausible mode for the binding of H₂ to the enzyme active site and the hydrogen oxidation mechanism. In the present work, we provide a new perspective on the MoCu-CODH hydrogenase activity by improving the in silico description of the enzyme. Energy refinement—by means of the BigQM approach—was performed on the intermediates involved in the dihydrogen oxidation catalysis reported in our previously published work (Rovaletti, et al. “Theoretical Insights into the Aerobic Hydrogenase Activity of Molybdenum–Copper CO Dehydrogenase.” Inorganics 7 (2019) 135). A suboptimal description of the H₂–HN(backbone) interaction was observed when the van der Waals parameters described in previous literature for H₂ were employed. Therefore, a new set of van der Waals parameters is developed here in order to better describe the hydrogen–backbone interaction. They give rise to improved binding modes of H₂ in the active site of MoCu CO dehydrogenase. Implications of the resulting outcomes for a better understanding of hydrogen oxidation catalysis mechanisms are proposed and discussed.

Two-Substrate Glyoxalase I Mechanism: A Quantum Mechanics/Molecular Mechanics Study

Glyoxalase I (GlxI) is an important enzyme that catalyzes the detoxification of methylglyoxal (MG) with the help of glutathione (H-SG). It is currently unclear whether MG and H-SG are substrates of GlxI or whether the enzyme processes hemithioacetal (HTA), which is nonenzymatically formed from MG and H-SG. Most previous studies have concentrated on the latter mechanism. Here, we study the two-substrate reaction mechanism of GlxI from humans (HuGlxI) and corn (ZmGlxI), which are Zn(II)-active and -inactive, respectively. Hybrid quantum mechanics/molecular mechanics calculations were used to obtain geometrical structures of the stationary points along reaction paths, and big quantum mechanical systems with more than 1000 atoms and free-energy perturbations were used to improve the quality of the calculated energies. We studied, on an equal footing, all reasonable reaction paths to the S- and R-enantiomers of HTA from MG and H-SG (the latter was considered in two different binding modes). The results indicate that the MG and H-SG reaction in both enzymes can follow the same path to reach S-HTA. However, the respective overall barriers and reaction energies are different for the two enzymes (6.1 and -9.8 kcal/mol for HuGlxI and 15.7 and -2.2 kcal/mol for ZmGlxI). The first reaction step to produce S-HTA is facilitated by a crystal water molecule that forms hydrogen bonds with a Glu and a Thr residue in the active site. The two enzymes also follow similar paths to R-HTA. However, the reactions reach a deprotonated and protonated R-HTA in the human and corn enzymes, respectively. The production of deprotonated R-HTA in HuGlxI is consistent with other theoretical and experimental works. However, our calculations show a different behavior for ZmGlxI (both S- and R-HTA can be formed in the enzyme with the alcoholic proton on HTA). This implies that Glu-144 of corn GlxI is not basic enough to keep the alcoholic proton. In HuGlxI, the two binding modes of H-SG that lead to S- and R-HTA are degenerate, but the barrier leading to R-HTA is lower than the barrier to S-HTA. On the other hand, ZmGlxI prefers the binding mode, which produces S-HTA; this observation is consistent with experiments. Based on the results, we present a modification for a previously proposed two-substrate reaction mechanism for ZmGlxI.

Kohn-Sham Density Functional Calculations Reveal Proton Wires in the Enolization and Carboxylase Reactions Catalyzed by Rubisco

Ribulose 1,5-bisphosphate (RuBP) carboxylase-oxygenase (Rubisco) plays a fundamental role in the carbon cycle by fixing the atmospheric CO₂ used in photosynthesis. Rubisco is all the more remarkable because it must catalyze some difficult multistep reaction chemistry involving proton transfers within the one active site. In the present study, we have used Kohn-Sham density functional theory at the B3LYP/6-31G* level with basis set superposition error and dispersion corrections (B3LYP-gCP-D3) to examine the possibility that the proton transfers can take place through molecular wires (including active-site water molecules) via the classical Grotthuss proton-shuttle mechanism. The results support an essential role for water molecules found in the crystal structures of Rubisco complexes as facilitators of proton transport in all the rate-limiting (catalytic) reaction steps through a network of short proton wires within the Rubisco active site. We suggest that completion of the initial product turnover (cycle) requires two excess protons produced in the initial carbamylation that is required for Rubisco activation. By use of proton wires, a large number of reaction steps may be accommodated within a single active site without necessitating the input of excessive conformational strain energy arising from the movement of residue side chains into positions where direct protonation of substrates can occur. The involvement of the identified types of proton wires in the kinetic mechanism is capable of providing a unique explanation for various experimental observations, including deuterium isotope effects and the results of site-directed mutagenesis experiments, and may thus provide a realistic solution to the problem of Rubisco's challenging chemistry.

Quantum Mechanics/Molecular Mechanics Study of the Reaction Mechanism of Glyoxalase I

Glyoxalase I (GlxI) is a member of the glyoxalase system, which is important in cell detoxification and converts hemithioacetals of methylglyoxal (a cytotoxic byproduct of sugar metabolism that may react with DNA or proteins and introduce nucleic acid strand breaks, elevated mutation frequencies, and structural or functional changes of the proteins) and glutathione into d-lactate. GlxI accepts both the S and R enantiomers of hemithioacetal, but converts them to only the S-d enantiomer of lactoylglutathione. Interestingly, the enzyme shows this unusual specificity with a rather symmetric active site (a Zn ion coordinated to two glutamate residues; Glu-99 and Glu-172), making the investigation of its reaction mechanism challenging. Herein, we have performed a series of combined quantum mechanics and molecular mechanics calculations to study the reaction mechanism of GlxI. The substrate can bind to the enzyme in two different modes, depending on the direction of its alcoholic proton (H2; toward Glu-99 or Glu-172). Our results show that the S substrate can react only if H2 is directed toward Glu-99 and the R substrate only if H2 is directed toward Glu-172. In both cases, the reactions lead to the experimentally observed S-d enantiomer of the product. In addition, the results do not show any low-energy paths to the wrong enantiomer of the product from neither the S nor the R substrate. Previous studies have presented several opposing mechanisms for the conversion of R and S enantiomers of the substrate to the correct enantiomer of the product. Our results confirm one of them for the S substrate, but propose a new one for the R substrate.

Using Atomic Confining Potentials for Geometry Optimization and Vibrational Frequency Calculations in Quantum-Chemical Models of Enzyme Active Sites

Quantum-chemical studies of enzymatic reaction mechanisms sometimes use truncated active-site models as simplified alternatives to mixed quantum mechanics molecular mechanics (QM/MM) procedures. Eliminating the MM degrees of freedom reduces the complexity of the sampling problem, but the trade-off is the need to introduce geometric constraints in order to prevent structural collapse of the model system during geometry optimizations that do not contain a full protein backbone. These constraints may impair the efficiency of the optimization, and care must be taken to avoid artifacts such as imaginary vibrational frequencies. We introduce a simple alternative in which terminal atoms of the model system are placed in soft harmonic confining potentials rather than being rigidly constrained. This modification is simple to implement and straightforward to use in vibrational frequency calculations, unlike iterative constraint-satisfaction algorithms, and allows the optimization to proceed without constraint even though the practical result is to fix the anchor atoms in space. The new approach is more efficient for optimizing minima and transition states, as compared to the use of fixed-atom constraints, and also more robust against unwanted imaginary frequencies. We illustrate the method by application to several enzymatic reaction pathways where entropy makes a significant contribution to the relevant reaction barriers. The use of confining potentials correctly describes reaction paths and facilitates calculation of both vibrational zero-point and finite-temperature entropic corrections to barrier heights.

Reaction Mechanism of [NiFe] Hydrogenase Studied by Computational Methods

[NiFe] hydrogenases catalyze the reversible conversion of molecular hydrogen to protons and electrons. This seemingly simple reaction has attracted much attention because of the prospective use of H₂ as a clean fuel. In this paper, we have studied the full reaction mechanism of this enzyme with various computational methods. Geometries were obtained with combined quantum mechanical and molecular mechanics (QM/MM) calculations. To get more accurate energies and obtain a detailed account of the surroundings, we performed big-QM calculations with 819 atoms in the QM region. Moreover, QM/MM thermodynamic cycle perturbation calculations were performed to obtain free energies. Finally, density matrix renormalisation group complete active space self-consistent field calculations were carried out to study the electronic structures of the various states in the reaction mechanism. Our calculations indicate that the Ni-L state is not involved in the reaction mechanism. Instead, the Ni-C state is reduced by one electron and then the bridging hydride ion is transferred to the sulfur atom of Cys546 as a proton and the two electrons transfer to the Ni ion. This step turned out to be rate-determining with an energy barrier of 58 kJ/mol, which is consistent with the experimental rate of 750 ± 90 s^-1 (corresponding to ∼52 kJ/mol). The cleavage of the H-H bond is facile with an energy barrier of 33 kJ/mol, according to our calculations. We also find that the reaction energies are sensitive to the size of the QM system, the basis set, and the density functional theory method, in agreement with previous studies.

Large-scale QM/MM free energy simulations of enzyme catalysis reveal the influence of charge transfer

Hybrid quantum mechanical-molecular mechanical (QM/MM) simulations provide key insights into enzyme structure-function relationships. Numerous studies have demonstrated that large QM regions are needed to systematically converge ground state, zero temperature properties with electrostatic embedding QM/MM. However, it is not well known if ab initio QM/MM free energy simulations have this same dependence, in part due to the hundreds of thousands of energy evaluations required for free energy estimations that in turn limit QM region size. Here, we leverage recent advances in electronic structure efficiency and accuracy to carry out range-separated hybrid density functional theory free energy simulations in a representative methyltransferase. By studying 200 ps of ab initio QM/MM dynamics for each of five QM regions from minimal (64 atoms) to one-sixth of the protein (544 atoms), we identify critical differences between large and small QM region QM/MM in charge transfer between substrates and active site residues as well as in geometric structure and dynamics that coincide with differences in predicted free energy barriers. Distinct geometric and electronic structure features in the largest QM region indicate that important aspects of enzymatic rate enhancement in methyltransferases are identified with large-scale electronic structure.

QM/MM study of the reaction mechanism of sulfite oxidase

Sulfite oxidase is a mononuclear molybdenum enzyme that oxidises sulfite to sulfate in many organisms, including man. Three different reaction mechanisms have been suggested, based on experimental and computational studies. Here, we study all three with combined quantum mechanical (QM) and molecular mechanical (QM/MM) methods, including calculations with large basis sets, very large QM regions (803 atoms) and QM/MM free-energy perturbations. Our results show that the enzyme is set up to follow a mechanism in which the sulfur atom of the sulfite substrate reacts directly with the equatorial oxo ligand of the Mo ion, forming a Mo-bound sulfate product, which dissociates in the second step. The first step is rate limiting, with a barrier of 39–49 kJ/mol. The low barrier is obtained by an intricate hydrogen-bond network around the substrate, which is preserved during the reaction. This network favours the deprotonated substrate and disfavours the other two reaction mechanisms. We have studied the reaction with both an oxidised and a reduced form of the molybdopterin ligand and quantum-refinement calculations indicate that it is in the normal reduced tetrahydro form in this protein.

Quantum chemical approaches to [NiFe] hydrogenase

The mechanism by which [NiFe] hydrogenase catalyses the oxidation of molecular hydrogen is a significant yet challenging topic in bioinorganic chemistry. With far-reaching applications in renewable energy and carbon mitigation, significant effort has been invested in the study of these complexes. In particular, computational approaches offer a unique perspective on how this enzyme functions at an electronic and atomistic level. In this article, we discuss state-of-the art quantum chemical methods and how they have helped deepen our comprehension of [NiFe] hydrogenase. We outline the key strategies that can be used to compute the (i) geometry, (ii) electronic structure, (iii) thermodynamics and (iv) kinetic properties associated with the enzymatic activity of [NiFe] hydrogenase and other bioinorganic complexes.

Quantum-Chemical Study of the Discrimination against dNTP in the Nucleotide Addition Reaction in the Active Site of RNA Polymerase II

Eukaryotic RNA polymerase II catalyzes the transcription of DNA into mRNA very efficiently and with an extremely low error rate with regard to matching base and sugar moiety. Despite its importance, little is known about how it discriminates against 2'-deoxy NTPs during the chemical reaction. To investigate the differences in the addition reactions of ATP and dATP, we used FF-MD and QM/MM calculations within a nudged elastic band approach, which allowed us to find the energetically accessible reaction coordinates. By converging the QM size, we found that 800 QM atoms are necessary to properly describe the active site. We show how the absence of a single hydrogen bond between the enzyme and the NTP 2'-OH group leads to an increase of the reaction barrier by 16 kcal/mol and therefore conclude that Arg446 is the key residue in the discrimination process.

Multiscale Modelling of Lytic Polysaccharide Monooxygenases

Lytic polysaccharide monooxygenase (LPMO) enzymes have attracted considerable attention owing to their ability to enhance polysaccharide depolymerization, making them interesting with respect to production of biofuel from cellulose. LPMOs are metalloenzymes that contain a mononuclear copper active site, capable of activating dioxygen. However, many details of this activation are unclear. Some aspects of the mechanism have previously been investigated from a computational angle. Yet, either these studies have employed only molecular mechanics (MM), which are inaccurate for metal active sites, or they have described only the active site with quantum mechanics (QM) and neglected the effect of the protein. Here, we employ hybrid QM and MM (QM/MM) methods to investigate the first steps of the LPMO mechanism, which is reduction of Cu^II to Cu^I and the formation of a Cu^II-superoxide complex. In the latter complex, the superoxide can bind either in an equatorial or an axial position. For both steps, we obtain structures that are markedly different from previous suggestions, based on small QM-cluster calculations. Our calculations show that the equatorial isomer of the superoxide complex is over 60 kJ/mol more stable than the axial isomer because it is stabilized by interactions with a second-coordination-sphere glutamine residue, suggesting a possible role for this residue. The coordination of superoxide in this manner agrees with recent experimental suggestions.

Influence of Coupling and Embedding Schemes on QM Size Convergence in QM/MM Approaches for the Example of a Proton Transfer in DNA

The influence of embedding and coupling schemes on the convergence of the QM size in the QM/MM approach is investigated for the transfer of a proton in a DNA base pair. We find that the embedding scheme (mechanical or electrostatic) has a much greater impact on the convergence behavior than the coupling scheme (additive QM/MM or subtractive ONIOM). To achieve size convergence, QM regions with up to 6000 atoms are necessary for pure QM or mechanical embedding. In contrast, electrostatic embedding converges faster: for the example of the transfer of a proton between DNA base pairs, we recommend including at least five base pairs and 5 Å of solvent (including counterions) into the QM region, i.e., a total of 1150 atoms.

H2binding to the active site of [NiFe] hydrogenase studied by multiconfigurational and coupled-cluster methods

CCSD(T) and DMRG-CASPT2 calculations show that H₂prefers to bind to Ni rather than to Fe in [NiFe] hydrogenase.

How Large Should the QM Region Be in QM/MM Calculations? The Case of Catechol O-Methyltransferase

Hybrid quantum mechanical-molecular mechanical (QM/MM) simulations are widely used in studies of enzymatic catalysis. Until recently, it has been cost prohibitive to determine the asymptotic limit of key energetic and structural properties with respect to increasingly large QM regions. Leveraging recent advances in electronic structure efficiency and accuracy, we investigate catalytic properties in catechol O-methyltransferase, a prototypical methyltransferase critical to human health. Using QM regions ranging in size from reactants-only (64 atoms) to nearly one-third of the entire protein (940 atoms), we show that properties such as the activation energy approach within chemical accuracy of the large-QM asymptotic limits rather slowly, requiring approximately 500-600 atoms if the QM residues are chosen simply by distance from the substrate. This slow approach to asymptotic limit is due to charge transfer from protein residues to the reacting substrates. Our large QM/MM calculations enable identification of charge separation for fragments in the transition state as a key component of enzymatic methyl transfer rate enhancement. We introduce charge shift analysis that reveals the minimum number of protein residues (approximately 11-16 residues or 200-300 atoms for COMT) needed for quantitative agreement with large-QM simulations. The identified residues are not those that would be typically selected using criteria such as chemical intuition or proximity. These results provide a recipe for a more careful determination of QM region sizes in future QM/MM studies of enzymes.

Binding-affinity predictions of HSP90 in the D3R Grand Challenge 2015 with docking, MM/GBSA, QM/MM, and free-energy simulations

We have estimated the binding affinity of three sets of ligands of the heat-shock protein 90 in the D3R grand challenge blind test competition. We have employed four different methods, based on five different crystal structures: first, we docked the ligands to the proteins with induced-fit docking with the Glide software and calculated binding affinities with three energy functions. Second, the docked structures were minimised in a continuum solvent and binding affinities were calculated with the MM/GBSA method (molecular mechanics combined with generalised Born and solvent-accessible surface area solvation). Third, the docked structures were re-optimised by combined quantum mechanics and molecular mechanics (QM/MM) calculations. Then, interaction energies were calculated with quantum mechanical calculations employing 970-1160 atoms in a continuum solvent, combined with energy corrections for dispersion, zero-point energy and entropy, ligand distortion, ligand solvation, and an increase of the basis set to quadruple-zeta quality. Fourth, relative binding affinities were estimated by free-energy simulations, using the multi-state Bennett acceptance-ratio approach. Unfortunately, the results were varying and rather poor, with only one calculation giving a correlation to the experimental affinities larger than 0.7, and with no consistent difference in the quality of the predictions from the various methods. For one set of ligands, the results could be strongly improved (after experimental data were revealed) if it was recognised that one of the ligands displaced one or two water molecules. For the other two sets, the problem is probably that the ligands bind in different modes than in the crystal structures employed or that the conformation of the ligand-binding site or the whole protein changes.

Converging ligand-binding free energies obtained with free-energy perturbations at the quantum mechanical level

In this article, the convergence of quantum mechanical (QM) free-energy simulations based on molecular dynamics simulations at the molecular mechanics (MM) level has been investigated. We have estimated relative free energies for the binding of nine cyclic carboxylate ligands to the octa-acid deep-cavity host, including the host, the ligand, and all water molecules within 4.5 Å of the ligand in the QM calculations (158-224 atoms). We use single-step exponential averaging (ssEA) and the non-Boltzmann Bennett acceptance ratio (NBB) methods to estimate QM/MM free energy with the semi-empirical PM6-DH2X method, both based on interaction energies. We show that ssEA with cumulant expansion gives a better convergence and uses half as many QM calculations as NBB, although the two methods give consistent results. With 720,000 QM calculations per transformation, QM/MM free-energy estimates with a precision of 1 kJ/mol can be obtained for all eight relative energies with ssEA, showing that this approach can be used to calculate converged QM/MM binding free energies for realistic systems and large QM partitions. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Asn47 and Phe114 modulate the inner sphere reorganization energies of type zero copper proteins

The geometric structures and electron transfer properties of type 1 Cu proteins are reasonably understood at the molecular level (E. I. Solomon and R. G. Hadt, Coord. Chem. Rev., 2011, 255, 774-789, J. J. Warren, K. M. Lancaster, J. H. Richards and H. B. Gray, J. Inorg. Biochem., 2012, 115, 119-126). Much understanding of type 1 copper electron transfer reactivity has come from site directed mutagenesis studies. For example, artificial "type zero" Cu-centres constructed in cupredoxin-azurin have showcased the capacity of outer-sphere hydrogen bonding networks to enhance Cu II/I electron transfer reactivity. In this paper, we have elaborated on earlier kinetics and electronic structural studies of type zero Cu by calculating the inner sphere reorganization energies of type 1, type 2, and type zero Cu proteins using density functional theory (DFT). Although the choice of density functionals for copper systems is not straightforward, we have benchmarked the density functionals against the recently reported ESI-PES data for two synthetic copper models (S. Niu, D.-L. Huang, P. D. Dau, H.-T. Liu, L.-S. Wang and T. J. Ichiye, Chem. Theory Comput., 2014, 10, 1283). For the Cu proteins, our calculations predict that changes in the coordination number upon metal reduction lead to large inner sphere reorganization energies for type 2 Cu sites, whereas retention in the coordination number is observed for type zero Cu sites. These variations in the coordination number are modulated by the outer-sphere coordinating residues Asn47 and Phe114, which are involved in hydrogen bonding with the Asp112 side chain.

A machine learning correction for DFT non-covalent interactions based on the S22, S66 and X40 benchmark databases

Background

Non-covalent interactions (NCIs) play critical roles in supramolecular chemistries; however, they are difficult to measure. Currently, reliable computational methods are being pursued to meet this challenge, but the accuracy of calculations based on low levels of theory is not satisfactory and calculations based on high levels of theory are often too costly. Accordingly, to reduce the cost and increase the accuracy of low-level theoretical calculations to describe NCIs, an efficient approach is proposed to correct NCI calculations based on the benchmark databases S22, S66 and X40 (Hobza in Acc Chem Rev 45: 663–672, 2012; Řezáč et al. in J Chem Theory Comput 8:4285, 2012).

Results

A novel type of NCI correction is presented for density functional theory (DFT) methods. In this approach, the general regression neural network machine learning method is used to perform the correction for DFT methods on the basis of DFT calculations. Various DFT methods, including M06-2X, B3LYP, B3LYP-D3, PBE, PBE-D3 and ωB97XD, with two small basis sets (i.e., 6-31G* and 6-31+G*) were investigated. Moreover, the conductor-like polarizable continuum model with two types of solvents (i.e., water and pentylamine, which mimics a protein environment with ε = 4.2) were considered in the DFT calculations. With the correction, the root mean square errors of all DFT calculations were improved by at least 70 %. Relative to CCSD(T)/CBS benchmark values (used as experimental NCI values because of its high accuracy), the mean absolute error of the best result was 0.33 kcal/mol, which is comparable to high-level ab initio methods or DFT methods with fairly large basis sets. Notably, this level of accuracy is achieved within a fraction of the time required by other methods. For all of the correction models based on various DFT approaches, the validation parameters according to OECD principles (i.e., the correlation coefficient R², the predictive squared correlation coefficient q² and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{cv}^{2}$$\end{document}qcv2 from cross-validation) were >0.92, which suggests that the correction model has good stability, robustness and predictive power.

Conclusions

The correction can be added following DFT calculations. With the obtained molecular descriptors, the NCIs produced by DFT methods can be improved to achieve high-level accuracy. Moreover, only one parameter is introduced into the correction model, which makes it easily applicable. Overall, this work demonstrates that the correction model may be an alternative to the traditional means of correcting for NCIs.Graphical abstract

Noncovalent interactions between the second coordination sphere and the active site of [NiFeSe] hydrogenase

QM/MM studies on seven truncated models of the oxidized as-isolated state of the [NiFeSe] Hases reveal the influence of the residues in the second coordination sphere on the active site.

Large-Scale First-Principles Molecular Dynamics Simulations with Electrostatic Embedding: Application to Acetylcholinesterase Catalysis

Enzymes are complicated solvated systems that typically require many atoms to simulate their function with any degree of accuracy. We have recently developed numerical techniques for large scale first-principles molecular dynamics simulations and applied them to the study of the enzymatic reaction catalyzed by acetylcholinesterase. We carried out density functional theory calculations for a quantum-mechanical (QM) subsystem consisting of 612 atoms with an O(N) complexity finite-difference approach. The QM subsystem is embedded inside an external potential field representing the electrostatic effect due to the environment. We obtained finite-temperature sampling by first-principles molecular dynamics for the acylation reaction of acetylcholine catalyzed by acetylcholinesterase. Our calculations show two energy barriers along the reaction coordinate for the enzyme-catalyzed acylation of acetylcholine. The second barrier (8.5 kcal/mol) is rate-limiting for the acylation reaction and in good agreement with experiment.

Carbohydrate-protein interactions: molecular modeling insights

The article reviews the significant contributions to, and the present status of, applications of computational methods for the characterization and prediction of protein-carbohydrate interactions. After a presentation of the specific features of carbohydrate modeling, along with a brief description of the experimental data and general features of carbohydrate-protein interactions, the survey provides a thorough coverage of the available computational methods and tools. At the quantum-mechanical level, the use of both molecular orbitals and density-functional theory is critically assessed. These are followed by a presentation and critical evaluation of the applications of semiempirical and empirical methods: QM/MM, molecular dynamics, free-energy calculations, metadynamics, molecular robotics, and others. The usefulness of molecular docking in structural glycobiology is evaluated by considering recent docking- validation studies on a range of protein targets. The range of applications of these theoretical methods provides insights into the structural, energetic, and mechanistic facets that occur in the course of the recognition processes. Selected examples are provided to exemplify the usefulness and the present limitations of these computational methods in their ability to assist in elucidation of the structural basis underlying the diverse function and biological roles of carbohydrates in their dialogue with proteins. These test cases cover the field of both carbohydrate biosynthesis and glycosyltransferases, as well as glycoside hydrolases. The phenomenon of (macro)molecular recognition is illustrated for the interactions of carbohydrates with such proteins as lectins, monoclonal antibodies, GAG-binding proteins, porins, and viruses.

The Importance of the MM Environment and the Selection of the QM Method in QM/MM Calculations: Applications to Enzymatic Reactions

In this chapter, we discuss the influence of an anisotropic protein environment on the reaction mechanisms of saccharopine reductase and uroporphyrinogen decarboxylase, respectively, via the use of a quantum mechanical and molecular mechanical (QM/MM) approach. In addition, we discuss the importance of selecting a suitable DFT functional to be used in a QM/MM study of a key intermediate in the mechanism of 8R-lipoxygenase, a nonheme iron enzyme. In the case of saccharopine reductase, while the enzyme utilizes a substrate-assisted catalytic pathway, it was found that only through treating the polarizing effect of the active site, via the use of an electronic embedding formalism, was agreement with experimental kinetic data obtained. Similarly, in the case of uroporphyrinogen decarboxylase, the effect of the protein environment on the catalytic mechanism was found to be such that the calculated rate-limiting barrier is in good agreement with related experimentally determined values for the first decarboxylation of the substrate. For 8R-lipoxygenase, it was found that the geometries and energies of the multicentered open-shell intermediate complexes formed during the mechanism are quite sensitive to the choice of the density functional theory method. Thus, while density functional theory has become the method of choice in QM/MM studies, care must be taken in the selection of a particular high-level method.