Quantum Chemical Modeling of Enzymatic Reactions: The Case of Decarboxylation

Predicting Carbonic Anhydrase Binding Affinity: Insights from QM Cluster Models

A systematic series of QM cluster models has been developed to predict the trend in the carbonic anhydrase binding affinity of a structurally diverse dataset of ligands. Reference DLPNO-CCSD(T)/CBS binding energies were generated for a cluster model and used to evaluate the performance of contemporary density functional theory methods, including Grimme's "3c" DFT composite methods (r2SCAN-3c and ωB97X-3c). It is demonstrated that when validated QM methods are used, the predictive power of the cluster models improves systematically with the size of the cluster models. This provided valuable insights into the key interactions that need to be modeled quantum mechanically and could inform how the QM region should be defined in hybrid quantum mechanics/molecular mechanics (QM/MM) models. The use of r2SCAN-3c on the largest cluster model composed of 16 residues appears to be an economical approach to predicting binding trends compared with using more robust DFT methods such as ωB97M-V and provides a significant improvement compared with docking.

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.

Theoretical Analysis of Selectivity Differences in Ketoreductases toward Aldehyde and Ketone Carbonyl Groups

Lactobacillus kefir alcohol dehydrogenase (LkADH) and ketoreductase from Chryseobacterium sp. CA49 (ChKRED12) exhibit different chemoselectivity and stereoselectivity toward a substrate with both keto and aldehyde carbonyl groups. LkADH selectively reduces the keto carbonyl group while retaining the aldehyde carbonyl group, producing optically pure R-alcohols. In contrast, ChKRED12 selectively reduces the aldehyde group and exhibits low reactivity toward ketone carbonyls. This study investigated the structural basis for these differences and the role of specific residues in the active site. Molecular dynamics (MD) simulations and quantum chemical calculations were used to investigate the interactions between the substrate and the enzymes and the essential cause of this phenomenon. The present study has revealed that LkADH and ChKRED12 exhibit significant differences in the structure of their respective active pockets, which is a crucial determinant of their distinct chemoselectivity toward the same substrate. Moreover, residues N89, N113, and E144 within LkADH as well as Q151 and D190 within ChKRED12 have been identified as key contributors to substrate stabilization within the active pocket through electrostatic interactions and van der Waals forces, followed by hydride transfer utilizing the coenzyme NADPH. Furthermore, the enantioselectivity mechanism of LkADH has been elucidated using quantum chemical methods. Overall, these findings not only provide fundamental insights into the underlying reasons for the observed differences in selectivity but also offer a detailed mechanistic understanding of the catalytic reaction.

How to Compute Atomistic Insight in DFT Clusters: The REG-IQA Approach

The relative energy gradient (REG) method is paired with the topological energy partitioning method interacting quantum atoms (IQA), as REG-IQA, to provide detailed and unbiased knowledge on the intra- and interatomic interactions. REG operates on a sequence of geometries representing a dynamical change of a system. Its recent application to peptide hydrolysis of the human immunodeficiency virus-1 (HIV-1) protease (PDB code: 4HVP) has demonstrated its full potential in recovering reaction mechanisms and through-space electrostatic and exchange-correlation effects, making it a compelling tool for analyzing enzymatic reactions. In this study, the computational efficiency of the REG-IQA method for the 133-atom HIV-1 protease quantum mechanical system is analyzed in every detail and substantially improved by means of three different approaches. The first approach of smaller integration grids for IQA integrations reduces the computational overhead by about a factor of 3. The second approach uses the line-simplification Ramer-Douglas-Peucker (RDP) algorithm, which outputs the minimal number of geometries necessary for the REG-IQA analysis for a predetermined root mean squared error (RMSE) tolerance. This cuts the computational time of the whole REG analysis by a factor of 2 if an RMSE of 0.5 kJ/mol is considered. The third approach consists of a "biased" or "unbiased" selection of a specific subset of atoms of the whole initial quantum mechanical model wave-function, which results in more than a 10-fold speed-up per geometry for the IQA calculation, without deterioration of the outcome of the REG-IQA analysis. Finally, to show the capability of these approaches, the findings gathered from the HIV-1 protease system are also applied to a different system named haloalcohol dehalogenase (HheC). In summary, this study takes the REG-IQA method to a computationally feasible and highly accurate level, making it viable for the analysis of a multitude of enzymatic systems.

Uncovering the Mechanism of Azepino-Indole Skeleton Formation via Pictet-Spengler Reaction by Strictosidine Synthase: A Quantum Chemical Investigation

Strictosidine synthase (STR) catalyzes the Pictet-Spengler (PS) reaction of tryptamine and secologanin to produce strictosidine. Recent studies demonstrated that the enzyme can also catalyze the reaction of non-natural substrates to form new alkaloid skeletons. For example, the PS condensation of 1H-indole-4-ethanamine with secologanin could be promoted by the STR from Rauvolfia serpentina (RsSTR) to generate a rare class of skeletons with a seven-membered ring, namely azepino-[3,4,5-cd]-indoles, which are precursors for the synthesis of new compounds displaying antimalarial activity. In the present study, the detailed reaction mechanism of RsSTR-catalyzed formation of the rare seven-membered azepino-indole skeleton through the PS reaction was revealed at the atomic level by quantum chemical calculations. The structures of the transition states and intermediates involved in the reaction pathway were optimized, and the energetics of the complete reaction were analyzed. Based on our calculation results, the most likely pathway of the enzyme-catalyzed reaction was determined, and the rate-determining step of the reaction was clarified. The mechanistic details obtained in the present study are important in understanding the promiscuous activity of RsSTR in the formation of the rare azepino-indole skeleton molecule and are also helpful in designing STR enzymes for the synthesis of other new alkaloid skeleton molecules.

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.

QM Calculations Revealed that Outer-Sphere Electron Transfer Boosted O-O Bond Cleavage in the Multiheme-Dependent Cytochrome bd Oxygen Reductase

The cytochrome bd oxygen reductase catalyzes the four-electron reduction of dioxygen to two water molecules. The structure of this enzyme reveals three heme molecules in the active site, which differs from that of heme-copper cytochrome c oxidase. The quantum chemical cluster approach was used to uncover the reaction mechanism of this intriguing metalloenzyme. The calculations suggested that a proton-coupled electron transfer reduction occurs first to generate a ferrous heme b₅₉₅. This is followed by the dioxygen binding at the heme d center coupled with an outer-sphere electron transfer from the ferrous heme b₅₉₅ to the dioxygen moiety, affording a ferric ion superoxide intermediate. A second proton-coupled electron transfer produces a heme d ferric hydroperoxide, which undergoes efficient O-O bond cleavage facilitated by an outer-sphere electron transfer from the ferrous heme b₅₉₅ to the O-O σ* orbital and an inner-sphere proton transfer from the heme d hydroxyl group to the leaving hydroxide. The synergistic benefits of the two types of hemes rationalize the highly efficient oxygen reduction repertoire for the multi-heme-dependent cytochrome bd oxygen reductase family.

Mechanism for the Halogenation and Azidation of Lysine Catalyzed by Non-heme Iron BesD Enzyme

Selective halogenation is important in synthetic chemistry. BesD, a new member of the non-heme Fe(II)/α-ketoglutarate (αKG)-dependent halogenase family, can activate the sp3 C-H bond and halogenate lysine, in particular without a carrier protein. Using the density functional calculations, a chlorination mechanism in BesD has been proposed, mainly including the formation of Cl-Fe(IV)=O through the αKG decarboxylation, the isomerization of Cl-Fe(IV)=O, the substrate hydrogen abstraction by Fe(IV)=O, and the rebound of chloro to the substrate carbon radical. The hydrogen abstraction is rate-limiting. The isomerization of Cl-Fe(IV)=O is essential for the hydrogen abstraction and the chiral selectivity. The BesD-catalyzed bromination and azidation of lysine adopt the same mechanism as the chlorination. The hardly-changed overall barriers indicate that the introduced ligands (X) do not affect the reaction rate significantly, implying that the X-introduced reactions catalyzed by BesD may be extended to other X anions.

Computational Survey of Recent Experimental Developments in the Hydroxylation Mechanism of Kynurenine 3-Monooxygenase

Recently, two new mechanistic proposals for the kynurenine 3-monooxygenase (KMO) catalyzed hydroxylation reaction of l-Kynurenine (l-Kyn) have been proposed. According to the first proposal, instead of the distal oxygen, the proximal oxygen of the hydroperoxide intermediate of flavin adenine dinucleotide (FAD) is transferred to the substrate ring. The second study proposes that l-Kyn participates in its base form in the reaction. To address these proposals, the reaction was reconsidered with a 386 atom quantum cluster model that is based on a recent X-ray structure (PDB id: 6FOX). The computations were carried out at the UB3LYP/6-311+G(2d,2p)//UB3LYP/6-31G(d,p) level with solvation (polarizable continuum model) and dispersion (DFT-D3(BJ)) corrections. To supplement the results of the density functional theory (DFT) calculations, molecular dynamics (MD) simulations of the protein-substrate complex were employed. The comparison of a proximal oxygen transfer mechanism to the distal oxygen transfer mechanism revealed that the former requires too high of a barrier energy while the latter validated our previous results. According to the MD simulations, the hydroperoxy moiety does not favor an alignment that might promote the proximal oxygen transfer mechanism. In the second part of the study, hydroxylation reaction with the base form of l-Kyn was sought. Although DFT calculations confirmed a much more facile reaction with the base form of l-Kyn, a mechanism which would allow the deprotonation of the l-Kyn before the oxygen transfer could not be determined with the X-ray-based positions. A concerted mechanism with both the oxygen transfer and the deprotonation required a high barrier energy. A stepwise mechanism involving the deprotonation of l-Kyn was found, starting from an MD frame. The overall barrier of the oxygen transfer step of this model was found to be in the range of that of with neutral l-Kyn. MD simulations supported the idea of ineffectiveness of the nearby shell surrounding the utilized active site core on the deprotonation of l-Kyn.

The Se-S Bond Formation in the Covalent Inhibition Mechanism of SARS-CoV-2 Main Protease by Ebselen-like Inhibitors: A Computational Study

The inhibition mechanism of the main protease (Mpro) of SARS-CoV-2 by ebselen (EBS) and its analog with a hydroxyl group at position 2 of the benzisoselenazol-3(2H)-one ring (EBS-OH) was studied by using a density functional level of theory. Preliminary molecular dynamics simulations on the apo form of Mpro were performed taking into account both the hydrogen donor and acceptor natures of the Nδ and Nε of His41, a member of the catalytic dyad. The potential energy surfaces for the formation of the Se-S covalent bond mediated by EBS and EBS-OH on Mpro are discussed in detail. The EBS-OH shows a distinctive behavior with respect to EBS in the formation of the noncovalent complex. Due to the presence of canonical H-bonds and noncanonical ones involving less electronegative atoms, such as sulfur and selenium, the influence on the energy barriers and reaction energy of the Minnesota hybrid meta-GGA functionals M06, M06-2X and M08HX, and the more recent range-separated hybrid functional wB97X were also considered. The knowledge of the inhibition mechanism of Mpro by the small protease inhibitors EBS or EBS-OH can enlarge the possibilities for designing more potent and selective inhibitor-based drugs to be used in combination with other antiviral therapies.

Mechanisms of metal-dependent non-redox decarboxylases from quantum chemical calculations

Quantum chemical calculations are today an extremely valuable tool for studying enzymatic reaction mechanisms. In this mini-review, we summarize our recent work on several metal-dependent decarboxylases, where we used the so-called cluster approach to decipher the details of the reaction mechanisms, including elucidation of the identity of the metal cofactors and the origins of substrate specificity.

Decarboxylases are of growing potential for biocatalytic applications, as they can be used in the synthesis of novel compounds of, e.g., pharmaceutical interest. They can also be employed in the reverse direction, providing a strategy to synthesize value‐added chemicals by CO₂ fixation. A number of non-redox metal-dependent decarboxylases from the amidohydrolase superfamily have been demonstrated to have promiscuous carboxylation activities and have attracted great attention in the recent years. The computational mechanistic studies provide insights that are important for the further modification and utilization of these enzymes in industrial processes. The discussed enzymes are: 5‐carboxyvanillate decarboxylase, γ‐resorcylate decarboxylase, 2,3‐dihydroxybenzoic acid decarboxylase, and iso-orotate decarboxylase.

Catalytic mechanism of the colistin resistance protein MCR-1

The mcr-1 gene encodes a membrane-bound Zn2+-metalloenzyme, MCR-1, which catalyses phosphoethanolamine transfer onto bacterial lipid A, making bacteria resistant to colistin, a last-resort antibiotic. Mechanistic understanding of this process remains incomplete. Here, we investigate possible catalytic pathways using DFT and ab initio calculations on cluster models and identify a complete two-step reaction mechanism. The first step, formation of a covalent phosphointermediate via transfer of phosphoethanolamine from a membrane phospholipid donor to the acceptor Thr285, is rate-limiting and proceeds with a single Zn2+ ion. The second step, transfer of the phosphoethanolamine group to lipid A, requires an additional Zn2+. The calculations suggest the involvement of the Zn2+ orbitals directly in the reaction is limited, with the second Zn2+ acting to bind incoming lipid A and direct phosphoethanolamine addition. The new level of mechanistic detail obtained here, which distinguishes these enzymes from other phosphotransferases, will aid in the development of inhibitors specific to MCR-1 and related bacterial phosphoethanolamine transferases.

Mechanistic Studies of Fatty Acid Activation by CYP152 Peroxygenases Reveal Unexpected Desaturase Activity

The majority of cytochrome P450 enzymes (CYPs) predominantly operate as monooxygenases, but recently a class of P450 enzymes was discovered, that can act as peroxygenases (CYP152). These enzymes convert fatty acids through oxidative decarboxylation, yielding terminal alkenes, and through α- and β-hydroxylation to yield hydroxy-fatty acids. Bioderived olefins may serve as biofuels, and hence understanding the mechanism and substrate scope of this class of enzymes is important. In this work, we report on the substrate scope and catalytic promiscuity of CYP OleT_JE and two of its orthologues from the CYP152 family, utilizing α-monosubstituted branched carboxylic acids. We identify α,β-desaturation as an unexpected dominant pathway for CYP OleT_JE with 2-methylbutyric acid. To rationalize product distributions arising from α/β-hydroxylation, oxidative decarboxylation, and desaturation depending on the substrate’s structure and binding pattern, a computational study was performed based on an active site complex of CYP OleT_JE containing the heme cofactor in the substrate binding pocket and 2-methylbutyric acid as substrate. It is shown that substrate positioning determines the accessibility of the oxidizing species (Compound I) to the substrate and hence the regio- and chemoselectivity of the reaction. Furthermore, the results show that, for 2-methylbutyric acid, α,β-desaturation is favorable because of a rate-determining α-hydrogen atom abstraction, which cannot proceed to decarboxylation. Moreover, substrate hydroxylation is energetically impeded due to the tight shape and size of the substrate binding pocket.

Computational Understanding of the Selectivities in Metalloenzymes

Metalloenzymes catalyze many different types of biological reactions with high efficiency and remarkable selectivity. The quantum chemical cluster approach and the combined quantum mechanics/molecular mechanics methods have proven very successful in the elucidation of the reaction mechanism and rationalization of selectivities in enzymes. In this review, recent progress in the computational understanding of various selectivities including chemoselectivity, regioselectivity, and stereoselectivity, in metalloenzymes, is discussed.

How does binuclear zinc amidohydrolase FwdA work in the initial step of methanogenesis: From formate to formyl-methanofuran

The initial step of methanogenesis is the fixation of CO₂ to formyl-methanofuran (formyl-MFR) catalyzed by formyl-MFR dehydrogenase, which can be divided into two half reactions. Herein, the second half reaction catalyzed by FwdA (formyl-methanofuran dehydrogenase subunit A), i.e., from formate to formyl-methanofuran, has been investigated using density functional theory and a chemical model based on the X-ray crystal structure. The calculations indicate that, compared with other well-known di-zinc hydrolases, the FwdA reaction employs a reverse mechanism, including the nucleophilic attack of MFR amine on formate carbon leading to a tetrahedral gem-diolate intermediate, two steps of proton transfer from amine to formate moieties assisted by the Asp385, and the CO bond dissociation to form the formyl-MFR product. The second step of proton transfer from the amine moiety to the Asp385 is rate-limiting with an overall barrier of 21.2 kcal/mol. The two zinc ions play an important role in stabilizing the transition states and intermediates, in particular the negative charge at the formate moiety originated from the nucleophilic attack of the MFR amine. The work here appends a crucial piece in the methanogenic mechanistics and advances the understanding of the global carbon cycle.

Asymmetric abstraction of two chemically-equivalent methylene hydrogens: significant enantioselectivity of endoperoxide presented by fumitremorgin B endoperoxidase

The unique enantioselectivity for an R-chiral endoperoxy ring in verruculogen biosynthesis originates from asymmetric abstraction of two chemically-equivalent methylene hydrogens.

Recent Trends in Quantum Chemical Modeling of Enzymatic Reactions

The quantum chemical cluster approach is a powerful method for investigating enzymatic reactions. Over the past two decades, a large number of highly diverse systems have been studied and a great wealth of mechanistic insight has been developed using this technique. This Perspective reviews the current status of the methodology. The latest technical developments are highlighted, and challenges are discussed. Some recent applications are presented to illustrate the capabilities and progress of this approach, and likely future directions are outlined.

Theoretical Study of the Mechanism of the Nonheme Iron Enzyme EgtB

EgtB is a nonheme iron enzyme catalyzing the C-S bond formation between γ-glutamyl cysteine (γGC) and N-α-trimethyl histidine (TMH) in the ergothioneine biosynthesis. Density functional calculations were performed to elucidate and delineate the reaction mechanism of this enzyme. Two different mechanisms were considered, depending on whether the sulfoxidation or the S-C bond formation takes place first. The calculations suggest that the S-O bond formation occurs first between the thiolate and the ferric superoxide, followed by homolytic O-O bond cleavage, very similar to the case of cysteine dioxygenase. Subsequently, proton transfer from a second-shell residue Tyr377 to the newly generated iron-oxo moiety takes place, which is followed by proton transfer from the TMH imidazole to Tyr377, facilitated by two crystallographically observed water molecules. Next, the S-C bond is formed between γGC and TMH, followed by proton transfer from the imidazole CH moiety to Tyr377, which was calculated to be the rate-limiting step for the whole reaction, with a barrier of 17.9 kcal/mol in the quintet state. The calculated barrier for the rate-limiting step agrees quite well with experimental kinetic data. Finally, this proton is transferred back to the imidazole nitrogen to form the product. The alternative thiyl radical attack mechanism has a very high barrier, being 25.8 kcal/mol, ruling out this possibility.

Exploring the Dependence of QM/MM Calculations of Enzyme Catalysis on the Size of the QM Region

Although QM/MM calculations are the primary current tool for modeling enzymatic reactions, the reliability of such calculations can be limited by the size of the QM region. Thus, we examine in this work the dependence of QM/MM calculations on the size of the QM region, using the reaction of catechol-O-methyl transferase (COMT) as a test case. Our study focuses on the effect of adding residues to the QM region on the activation free energy, obtained with extensive QM/MM sampling. It is found that the sensitivity of the activation barrier to the size of the QM is rather limited, while the dependence of the reaction free energy is somewhat larger. Of course, the results depend on the inclusion of the first solvation shell in the QM regions. For example, the inclusion of the Mg²⁺ ion can change the activation barrier due to charge transfer effects. However, such effects can easily be included in semiempirical approaches by proper parametrization. Overall, we establish that QM/MM calculations of activation barriers of enzymatic reactions are not highly sensitive to the size of the QM region, beyond the immediate region that describes the reacting atoms.

The Use of Multiscale Molecular Simulations in Understanding a Relationship between the Structure and Function of Biological Systems of the Brain: The Application to Monoamine Oxidase Enzymes

HIGHLIGHTS Computational techniques provide accurate descriptions of the structure and dynamics of biological systems, contributing to their understanding at an atomic level.Classical MD simulations are a precious computational tool for the processes where no chemical reactions take place.QM calculations provide valuable information about the enzyme activity, being able to distinguish among several mechanistic pathways, provided a carefully selected cluster model of the enzyme is considered.Multiscale QM/MM simulation is the method of choice for the computational treatment of enzyme reactions offering quantitative agreement with experimentally determined reaction parameters.Molecular simulation provide insight into the mechanism of both the catalytic activity and inhibition of monoamine oxidases, thus aiding in the rational design of their inhibitors that are all employed and antidepressants and antiparkinsonian drugs. Aging society and therewith associated neurodegenerative and neuropsychiatric diseases, including depression, Alzheimer's disease, obsessive disorders, and Parkinson's disease, urgently require novel drug candidates. Targets include monoamine oxidases A and B (MAOs), acetylcholinesterase (AChE), butyrylcholinesterase (BChE), and various receptors and transporters. For rational drug design it is particularly important to combine experimental synthetic, kinetic, toxicological, and pharmacological information with structural and computational work. This paper describes the application of various modern computational biochemistry methods in order to improve the understanding of a relationship between the structure and function of large biological systems including ion channels, transporters, receptors, and metabolic enzymes. The methods covered stem from classical molecular dynamics simulations to understand the physical basis and the time evolution of the structures, to combined QM, and QM/MM approaches to probe the chemical mechanisms of enzymatic activities and their inhibition. As an illustrative example, the later will focus on the monoamine oxidase family of enzymes, which catalyze the degradation of amine neurotransmitters in various parts of the brain, the imbalance of which is associated with the development and progression of a range of neurodegenerative disorders. Inhibitors that act mainly on MAO A are used in the treatment of depression, due to their ability to raise serotonin concentrations, while MAO B inhibitors decrease dopamine degradation and improve motor control in patients with Parkinson disease. Our results give strong support that both MAO isoforms, A and B, operate through the hydride transfer mechanism. Relevance of MAO catalyzed reactions and MAO inhibition in the context of neurodegeneration will be discussed.

Carnosine and Homocarnosine Degradation Mechanisms by the Human Carnosinase Enzyme CN1: Insights from Multiscale Simulations

The endogenous dipeptide l-carnosine, and its derivative homocarnosine, prevent and reduce several pathologies like amytrophic lateral sclerosis (ALS), Alzheimer's disease, and Parkinson's disease. Their beneficial action is severely hampered because of the hydrolysis by carnosinase enzymes, in particular the human carnosinase, hCN1. This belongs to the metallopeptidase M20 family, where a cocatalytic active site is formed by two Zn(2+) ions, bridged by a hydroxide anion. The protein may exist as a monomer and as a dimer in vivo. Here we used hybrid quantum mechanics/molecular mechanics simulations based on the dimeric apoenzyme's structural information to predict the Michaelis complexes with l-carnosine and its derivative homocarnosine. On the basis of our calculations, we suggest that (i) l-carnosine degradation occurs through a nucleophilic attack of a Zn(2+)-coordinated bridging moiety for both monomer and dimer. This mechanistic hypothesis for hCN1 catalysis differs from previous proposals, while it is in agreement with available experimental data. (ii) The experimentally measured higher affinity of homocarnosine for the enzyme relative to l-carnosine might be explained, at least in part, by more extensive interactions inside the monomeric and dimeric hCN1's active site. (iii) Hydrogen bonds at the binding site, present in the dimer but absent in the monomer, might play a role in the experimentally observed higher activity of the dimeric form. Investigations of the enzymatic reaction are required to establish or disprove this hypothesis. Our results may serve as a basis for the design of potent hCN1 inhibitors.

Towards a barrier height benchmark set for biologically relevant systems

We have collected computed barrier heights and reaction energies (and associated model structures) for five enzymes from studies published by Himo and co-workers. Using this data, obtained at the B3LYP/6- 311+G(2d,2p)[LANL2DZ]//B3LYP/6-31G(d,p) level of theory, we then benchmark PM6, PM7, PM7-TS, and DFTB3 and discuss the influence of system size, bulk solvation, and geometry re-optimization on the error. The mean absolute differences (MADs) observed for these five enzyme model systems are similar to those observed for PM6 and PM7 for smaller systems (10-15 kcal/mol), while DFTB results in a MAD that is significantly lower (6 kcal/mol). The MADs for PMx and DFTB3 are each dominated by large errors for a single system and if the system is disregarded the MADs fall to 4-5 kcal/mol. Overall, results for the condensed phase are neither more or less accurate relative to B3LYP than those in the gas phase. With the exception of PM7-TS, the MAD for small and large structural models are very similar, with a maximum deviation of 3 kcal/mol for PM6. Geometry optimization with PM6 shows that for one system this method predicts a different mechanism compared to B3LYP/6-31G(d,p). For the remaining systems, geometry optimization of the large structural model increases the MAD relative to single points, by 2.5 and 1.8 kcal/mol for barriers and reaction energies. For the small structural model, the corresponding MADs decrease by 0.4 and 1.2 kcal/mol, respectively. However, despite these small changes, significant changes in the structures are observed for some systems, such as proton transfer and hydrogen bonding rearrangements. The paper represents the first step in the process of creating a benchmark set of barriers computed for systems that are relatively large and representative of enzymatic reactions, a considerable challenge for any one research group but possible through a concerted effort by the community. We end by outlining steps needed to expand and improve the data set and how other researchers can contribute to the process.

Mechanism and selectivity of the dinuclear iron benzoyl-coenzyme A epoxidase BoxB

DFT calculations are used to elucidate the reaction mechanism and selectivity of BoxB catalyzed benzoyl-CoA epoxidation.

Benzoyl-CoA epoxidase is a dinuclear iron enzyme that catalyzes the epoxidation reaction of the aromatic ring of benzoyl-CoA with chemo-, regio- and stereo-selectivity. It has been suggested that this enzyme may also catalyze the deoxygenation reaction of epoxide, suggesting a unique bifunctionality among the diiron enzymes. We report a density functional theory study of this enzyme aimed at elucidating its mechanism and the various selectivities. The epoxidation is suggested to start with the binding of the O₂ molecule to the diferrous center to generate a diferric peroxide complex, followed by concerted O–O bond cleavage and epoxide formation. Two different pathways have been located, leading to (2S,3R)-epoxy and (2R,3S)-epoxy products, with barriers of 17.6 and 20.4 kcal mol^–1, respectively. The barrier difference is 2.8 kcal mol^–1, corresponding to a diastereomeric excess of about 99 : 1. Further isomerization from epoxide to phenol is found to have quite a high barrier, which cannot compete with the product release step. After product release into solution, fast epoxide–oxepin isomerization and racemization can take place easily, leading to a racemic mixture of (2S,3R) and (2R,3S) products. The deoxygenation of epoxide to regenerate benzoyl-CoA by a diferrous form of the enzyme proceeds via a stepwise mechanism. The C2–O bond cleavage happens first, coupled with one electron transfer from one iron center to the substrate, to form a radical intermediate, which is followed by the second C3–O bond cleavage. The first step is rate-limiting with a barrier of only 10.8 kcal mol^–1. Further experimental studies are encouraged to verify our results.

Theoretical study on the chemical mechanism of enoyl-CoA hydratase and the form of inhibitor binding

Enoyl-CoA hydratase (ECH) catalyzes the second step in the vital β-oxidation pathway of fatty acid metabolism. This enzyme catalyzes the syn-addition of a water molecule across the double bond of 4-(N,N-dimethylamino) cinnamoyl-CoA (DAC-CoA). In this work, the reaction mechanisms of ECH were investigated using the density functional theory (DFT) methods. The different protonation states in which the important residues Glu164 and Glu144 are either neutral or ionized were considered. Four models of the active site were designed based on the X-ray crystal structure of the enzyme. The calculations gave strong support to the proposed mechanism and confirmed that both Glu164 and Glu144 are in a deprotonated state in the reaction mechanism of ECH. In addition, we constructed a model of the active site with the inhibitor acetoacetyl-CoA based on the crystal structure. Caomparison of the calculated energy barriers showed that binding of the keto-enol form of the inhibitor is more reasonable than that of the di-keto form in the inhibition process. Moreover, acetoacetyl-CoA was found to exhibit a keto-enol tautomerism when it acts as an inhibitor in the reaction. The present theoretical results indicated that both residues Glu164 and Glu144 are unprotonated in ECH with the substrate bound, while only Glu164 is unprotonated when the inhibitor binds ECH.

Analyses of cobalt-ligand and potassium-ligand bond lengths in metalloproteins: trends and patterns

Cobalt and potassium are biologically important metal elements that are present in a large array of proteins. Cobalt is mostly found in vivo associated with a corrin ring, which represents the core of the vitamin B12 molecule. Potassium is the most abundant metal in the cytosol, and it plays a crucial role in maintaining membrane potential as well as correct protein function. Here, we report a thorough analysis of the geometric properties of cobalt and potassium coordination spheres that was performed with high resolution on a representative set of structures from the Protein Data Bank and complemented by quantum mechanical calculations realized at the DFT level of theory (B3LYP/ SDD) on mononuclear model systems. The results allowed us to draw interesting conclusions on the structural characteristics of both Co and K centers, and to evaluate the importance of effects such as their association energies and intrinsic thermodynamic stabilities. Overall, the results obtained provide useful data for enhancing the atomic models normally applied in theoretical and computational studies of Co or K proteins performed at the quantum mechanical level, and for developing molecular mechanical parameters for treating Co or K coordination spheres in molecular mechanics or molecular dynamics studies.

Phosphate monoester hydrolysis by trinuclear alkaline phosphatase; DFT study of transition States and reaction mechanism

Alkaline phosphatase (AP) is a trinuclear metalloenzyme that catalyzes the hydrolysis of a broad range of phosphate monoesters to form inorganic phosphate and alcohol (or phenol). In this paper, by using density functional theory with a model based on a crystal structure, the AP-catalyzed hydrolysis of phosphate monoesters is investigated by calculating two substrates, that is, methyl and p-nitrophenyl phosphates, which represent alkyl and aryl phosphates, respectively. The calculations confirm that the AP reaction employs a "ping-pong" mechanism involving two chemical displacement steps, that is, the displacement of the substrate leaving group by a Ser102 alkoxide and the hydrolysis of the phosphoseryl intermediate by a Zn2-bound hydroxide. Both displacement steps proceed via a concerted associative pathway no matter which substrate is used. Other mechanistic aspects are also studied. Comparison of our calculations with linear free energy relationships experiments shows good agreement.

A dominant homolytic O-Cl bond cleavage with low-spin triplet-state Fe(IV)=O formed is revealed in the mechanism of heme-dependent chlorite dismutase

Chlorite dismutase (Cld) is a heme-dependent enzyme that catalyzes the decomposition of toxic chlorite (ClO2(-)) into innocuous chloride and O2. In this paper, using the hybrid B3LYP density functional theory (DFT) method including dispersion interactions, the Cld reaction mechanism has been studied with a chemical model constructed on the X-ray crystal structure. The calculations indicate that the reaction proceeds along a stepwise pathway in the doublet state, i.e. a homolytic O-Cl bond cleavage of the substrate leading to an O-Fe(heme) species and a ClO˙ radical, followed by a rebinding O-O bond formation between them. The O-Fe(heme) species is demonstrated to be a low-spin triplet-state Fe(IV)=O diradicaloid. A low-spin singlet-state Fe(IV)=O is much less stable than the former, with an energy difference of 9.2 kcal mol(-1). The O-Cl bond cleavage is rate-limiting with a barrier of 10.6 kcal mol(-1), in good agreement with the experimental reaction rate of 2.0 × 10(5) s(-1). Furthermore, a heterolytic O-Cl bond dissociation in the initial step is shown to be unreachable, which ensures the high efficiency of the Cld enzyme by avoiding the generation of chlorate byproduct observed in the reactions of synthetic Fe porphyrins. Also, the pathways in the quartet and sextet states are unfavorable for the Cld reaction. The present results reveal a detailed mechanism III (defined in the text) including an interesting di-radical intermediate composed of a low-spin triplet-state Fe(IV)=O and a ClO˙ radical. Compared to a competitive heterolytic Cl-O cleavage in synthetic Fe porphyrins, the revelation of the domination of homolysis in Cld indicates not only the high efficiency of enzyme, but also the sensitivity of a heme and the significance of the enzymatic active-site surroundings (the His170 and Arg183 residues in the present case), which gives more insights into heme chemistry.

Convergence in the QM-only and QM/MM modeling of enzymatic reactions: A case study for acetylene hydratase

We report systematic quantum mechanics-only (QM-only) and QM/molecular mechanics (MM) calculations on an enzyme-catalyzed reaction to assess the convergence behavior of QM-only and QM/MM energies with respect to the size of the chosen QM region. The QM and MM parts are described by density functional theory (typically B3LYP/def2-SVP) and the CHARMM force field, respectively. Extending our previous work on acetylene hydratase with QM regions up to 157 atoms (Liao and Thiel, J. Chem. Theory Comput. 2012, 8, 3793), we performed QM/MM geometry optimizations with a QM region M4 composed of 408 atoms, as well as further QM/MM single-point calculations with even larger QM regions up to 657 atoms. A charge deletion analysis was conducted for the previously used QM/MM model (M3a, with a QM region of 157 atoms) to identify all MM residues with strong electrostatic contributions to the reaction energetics (typically more than 2 kcal/mol), which were then included in M4. QM/MM calculations with this large QM region M4 lead to the same overall mechanism as the previous QM/MM calculations with M3a, but there are some variations in the relative energies of the stationary points, with a mean absolute deviation (MAD) of 2.7 kcal/mol. The energies of the two relevant transition states are close to each other at all levels applied (typically within 2 kcal/mol), with the first (second) one being rate-limiting in the QM/MM calculations with M3a (M4). QM-only gas-phase calculations give a very similar energy profile for QM region M4 (MAD of 1.7 kcal/mol), contrary to the situation for M3a where we had previously found significant discrepancies between the QM-only and QM/MM results (MAD of 7.9 kcal/mol). Extension of the QM region beyond M4 up to M7 (657 atoms) leads to only rather small variations in the relative energies from single-point QM-only and QM/MM calculations (MAD typically about 1-2 kcal/mol). In the case of acetylene hydratase, a model with 408 QM atoms thus seems sufficient to achieve convergence in the computed relative energies to within 1-2 kcal/mol.

Selecting DFT methods for use in optimizations of enzyme active sites: applications to ONIOM treatments of DNA glycosylases

When using a hybrid methodology to treat an enzymatic reaction, many factors contribute to selecting the method for the high-level region, which can be complicated by the presence of dispersion-driven interactions such as π–π stacking. In addition, the proper treatment of the reaction center often requires a large number of heavy atoms to be included in the high-level region, precluding the use of ab initio methods such as MP2 as well as large basis sets, in the optimization step. In the present work, popular DFT methods were tested to identify an appropriate functional for treating the high-level region in ONIOM optimizations of reactions catalyzed by nonmetalloenzymes. Eight different DFT methods (B3LYP, B97-2, MPW1K, MPWB1K, BB1K, B1B95, M06-2X, and ωB97X-D) in combination with four double-ζ quality Pople basis sets were tested for their ability to optimize noncovalent interactions (hydrogen bonding and π–π) and characterize reactions (proton transfer, SN2 hydrolysis, and unimolecular cleavage). Although the primary focus of this study is accurate structure determination, energetics were also examined at both the optimization level of theory, and with triple-ζ quality basis set and select (M06-2X or ωB97X-D) methods. If dispersion-driven interactions exist within the active site, then MPWB1K/6-31G(d,p) or M06-2X/6-31+G(d,p) are recommended for the optimization step with subsequent triple-ζ quality single-point energies. However, since dispersion-corrected functionals (M06-2X and ωB97X-D) generally require diffuse functions to yield appropriate geometries, the possible size of the high-level region is greatly limited with these methods. In contrast, if the model is large enough to recover steric constraints on π–π interactions, then B3LYP with a small basis set performs comparatively well for the optimization step and is significantly less computationally expensive. Interestingly, the functionals that afford the best geometries often do not yield the best energetics, which emphasizes the importance of structural benchmark studies.

Computational modeling of the direct hydride transfer mechanism for the MAO catalyzed oxidation of phenethylamine and benzylamine: ONIOM (QM/QM) calculations

Monoamine oxidases are two isozymic flavoenzymes which are the important targets for drugs used in the treatment of depression, Parkinson and Alzheimer's diseases. The catalytic reaction taking place between the cofactor FAD and amine substrate is still not completely understood. Herein we employed quantum chemical methods on the recently proposed direct hydride transfer mechanism including full active site residues of MAO isoforms in the calculations. Activation free energy barriers of direct hydride transfer mechanism for MAO-A and MAO-B were calculated by ONIOM (our own n-layered integrated molecular orbital + molecular mechanics) method with QM/QM (quantum mechanics:quantum mechanics) approach employing several density functional theory functionals, B3LYP, WB97XD, CAM-B3LYP and M06-2X, for the high layer. The formation of very recently proposed αC-flavin N5 adduct inside the enzyme has been investigated. ONIOM (M06-2X/6-31+G(d,p):PM6) results revealed that such an adduct may form only in MAO-B suggesting slightly different hydride transfer mechanisms for MAO-A and MAO-B.