Epoxide hydrolases are essential enzymes that convert epoxides into 1,2-diols, contributing to detoxification, metabolism, and signaling in a wide range of organisms. In this study, we employed hybrid quantum mechanics/molecular mechanics (QM/MM) calculations to investigate the catalytic mechanism of potato epoxide hydrolase 1 (StEH1), specifically focusing on the hydrolysis of trans-stilbene oxide into chiral diols. Based on the crystal structure of StEH1 (PDB ID: 2CJP), we modeled enzyme-substrate interactions and examined the roles of the catalytic triad (Asp105, His300, and Asp265), two tyrosine residues (Tyr154 and Tyr235) involved in substrate polarization, and a crystallographic water molecule acting as the hydrolytic nucleophile. Our QM/MM calculations revealed a three-step reaction mechanism: alkylation, dealkylation, and proton relay. We also determined the optimal protonation states of several active-site residues, particularly His104 and His300, to ensure an accurate mechanistic picture. These results clarify previously debated aspects of the mechanism, such as the protonation state of His300 and the function of the tyrosine residues, and provide new insights into substrate specificity and offer valuable information for future efforts in inhibitor design and enzyme engineering.

We have compared the performance of 64 different computational methods, based on combined quantum mechanical (QM) and molecular mechanical (QM/MM) or QM-cluster calculations in a continuum solvent, to estimate the acid constant (pKa) of metal-bound ligands in proteins. As a calibration set, we use 12 experimental pKa values from six different proteins that involve Zn2+, Fe3+, or Fe4+. We employ two different density functional theory (DFT) methods (TPSS and B3LYP), two basis sets (def2-SV(P) and def2-TZVPD), QM regions of three different sizes (∼40, ∼100, and ∼350 atoms), relaxed or fixed surroundings, and three different values of the dielectric constant of the continuum-solvation model (ε = 4, 20, or 80). The results clearly show that QM-cluster+continuum-solvation is much better than QM/MM. In general, the most accurate results are obtained with ε = 80 and the minimal QM region. The two DFT methods, the two basis sets, and relaxing or fixing the surroundings give similar results. The best-performing method is TPSS with the minimal QM region, def2-TZVPD, relaxed surroundings, and ε = 80, yielding a mean absolute deviation (after removal of a systematic error of 11.6 pKa units, pu) of 2.0 pu and a maximum deviation of 5.0 pu. The coefficient of determination (R2) and Kendall's τ with respect to the experimental pKa values are both 0.64, while Spearman's rank correlation coefficient is 0.78. This level of accuracy should be sufficient to reliably determine the protonation states of metal-bound ligands in QM-based studies of enzymatic reaction mechanisms.

Quantum refinement (QR) is an approach in which the empirical restraints used in standard structural refinement to ensure that the details of the structure, e.g. bond lengths and angles, make chemical sense are replaced by more accurate quantum mechanical calculations for a small but interesting part of the structure. QR has previously been used for X-ray and neutron crystallography, cryogenic electron microscopy, nuclear magnetic resonance, and extended X-ray absorption fine structure. Here, QR is used for the first time for X-ray free-electron laser (XFEL) crystallography and microcrystal electron diffraction (MicroED). As a test case, we use six structures of the R2a protein of ribonucleotide reductase, concentrating on the binuclear Fe-2 site in either the oxidized (Fe-2(III)) or reduced (Fe-2(II)) state, two each from single-crystal X-ray (SCX) crystallography, XFEL crystallography or MicroED. The results show that QR works well for data from all three radiation sources, even though scattering factors for neutral atoms had to be used for MicroED. QR corrects unrealistically short Fe-O distances in the reduced SCX structure and gives improved real-space Z scores for the reduced MicroED structure. The three methods give similar structures, apart from variation in the weak water ligands and in the binding of carboxylate groups (monodentate, bidentate or a mixture). By performing QR for three protonation states of the bridging solvent molecule, we could show that it is undoubtedly a water molecule in the reduced XFEL and MicroED structures (it is not present in the SCX structure) and that it is not water in the oxidized structures. The XFEL data indicate that it is O2- in the oxidized XFEL structure, in agreement with the spectroscopic results. However, for the SCX structure, O2- and OH- give comparable results, whereas OH- is slightly preferred in the MicroED structure. This indicates that the SCX and MicroED structures may be partly photoreduced during data collection.

The use of X-ray structures to determine and interpret the ferryl iron-oxygen bond order in molecular oxygen-activating heme enzymes has, in the past, been controversial. This has mainly stemmed from the susceptibility of ferryl species to X-ray-induced electronic state changes. In this work we establishe using time-resolved serial femtosecond X-ray crystallography (tr-SFX) on a dye-decolourising peroxidase that the ferryl intermediate species (Compounds I and II) captured following in situ mixing of microcrystals with H2O2 have single, rather than the double bond character expected. X-ray emission validated tr-SFX data with quantum refinement, time-dependent-DFT calculations and QM/MM geometry optimizations together support the concept that the single iron-oxygen bond character is not an indication of ferryl reduction or a protonated form (FeIV-OH) but is instead attributed to the existence of accessible excited states possessing ferric-oxyl (FeIII-O•-) character. Such states offer insight into the nature of ferryl heme.

[FeFe]-hydrogenases are highly efficient enzymes in the reversible catalysis of molecular hydrogen production and oxidation. Their active site, the H-cluster, consists of a [4Fe-4S]H subcluster linked to a binuclear [2Fe]H organometallic unit. Many [FeFe]-hydrogenases, such as the one from Desulfovibrio desulfuricans (DdHydAB), possess accessory Fe-S clusters (F and F') that mediate electron transfer. This study employs hybrid quantum mechanics/molecular mechanics (QM/MM) methods to characterize the electronic structure and thermodynamic landscape associated with redox and protonation events in the complete Fe-S cluster network of DdHydAB. Our calculations indicate that the F' cluster plays a key role in the initial reduction of the oxidized resting state, acting as the preferential site for the accumulation of the first electron. Analysis of protonated states, upon reduction events, reveals a strong correlation between protonation and electron transfer (PCET), with protonation at the H-cluster inducing electron transfer from the F' cluster to the H-cluster. Calculations indicate that the formation of a terminal hydride is energetically favored over ADT protonation, and subsequent isomerization to a bridging hydride (μ-H) is further stabilizing, albeit potentially kinetically limiting. The study highlights how accessory clusters influence the electronic distribution and redox properties of the H-cluster, underscoring the importance of considering the entire Fe-S cluster system for a complete understanding of the catalytic mechanism of [FeFe]-hydrogenases.

We have compared the performance of 64 different computational methods, based on combined quantum mechanical (QM) and molecular mechanical (QM/MM) or QM-cluster calculations in a continuum solvent, to estimate the acid constant (pK a) of metal-bound ligands in proteins. As a calibration set, we use 12 experimental pK a values from six different proteins that involve Zn2+, Fe3+, or Fe4+. We employ two different density functional theory (DFT) methods (TPSS and B3LYP), two basis sets (def2-SV(P) and def2-TZVPD), QM regions of three different sizes (similar to 40, similar to 100, and similar to 350 atoms), relaxed or fixed surroundings, and three different values of the dielectric constant of the continuum-solvation model (epsilon = 4, 20, or 80). The results clearly show that QM-cluster+continuum-solvation is much better than QM/MM. In general, the most accurate results are obtained with epsilon = 80 and the minimal QM region. The two DFT methods, the two basis sets, and relaxing or fixing the surroundings give similar results. The best-performing method is TPSS with the minimal QM region, def2-TZVPD, relaxed surroundings, and epsilon = 80, yielding a mean absolute deviation (after removal of a systematic error of 11.6 pK a units, pu) of 2.0 pu and a maximum deviation of 5.0 pu. The coefficient of determination (R 2) and Kendall's tau with respect to the experimental pK a values are both 0.64, while Spearman's rank correlation coefficient is 0.78. This level of accuracy should be sufficient to reliably determine the protonation states of metal-bound ligands in QM-based studies of enzymatic reaction mechanisms.

In 2021, neutron structures of oxidised and reduced human manganese superoxide dismutase were published, suggesting several unexpected features, including deprotonated glutamate, tyrosine and histidine residues and a reduced Mn ion with two hydroxide ligands. We have used quantum refinement to evaluate whether alternative interpretations of the structures are possible, comparing many structural models with different possible protonation states or other structural interpretations. Quantum refinement is standard crystallographic refinement in which the empirical restraints, which are employed to ensure that the structure makes chemical sense and gives reasonable bond lengths and angles, are replaced by quantum mechanical calculations for a small but interesting part of the structure. We show that in all cases, there are more chemically reasonable interpretations of the structures, not involving any deprotonated residues, which give slightly improved structures in terms of real-space Z scores based on the difference maps and strain energies. Weak nuclear densities do not necessarily indicate that an atom is not present; instead it may indicate that the atom may have several conformations or shows extensive dynamics. For the Tyr-34 residue, our results indicate that the phenolic hydroxide H atom should preferably be within the aromatic ring plane (by 20 kJ/mol), but with two possible conformations. Moreover, an Mn-bound water molecule can readily accept a hydrogen bond from the nearby Gln-143 residue.

Hirshfeld atom refinement (HAR) provides a more realistic interpretation of crystallographic data than the standard independent atom model (IAM) by using aspherical atomic form factors derived from quantum mechanical (QM) calculations. With this aspherical description, it is possible to obtain improved atomic positions, atomic displacement parameters and correct bond lengths even for hydrogen atoms. Unfortunately, HAR is computationally very demanding for larger molecules. Recently, we suggested how this can be solved by calculating aspherical atomic form factors for small overlapping fragments of the system, the fragHAR approach. Here, we have created a new implementation of fragHAR in Olex2 within the NoSpherA2 interface. We have also solved previous issues with hydrogen bonds by automatically extending the fragments with all hydrogen-bond acceptors. This implementation was successfully tested on three oligopeptides, demonstrating that fragHAR yields indistinguishable results in terms of atomic charge, residual density or R values compared with full HAR. Subsequently, fragHAR was applied to the proteins crambin and rubredoxin, with 843 and 1014 atoms, respectively, showing improved results in terms of egross, which decreases from 0.350 with the IAM to 0.318 with fragHAR for crambin, and from 0.195 to 0.176 for rubredoxin, although it turned out to be necessary to keep all bond lengths involving hydrogen atoms constrained for the latter protein. FragHAR shows near-linear scaling and 46-fold speedup for rubredoxin compared with HAR. It also provides a convenient solution to alternative conformations and positional disorder, which cause an exponential increase in the time consumption of the conventional HAR approach. The successful refinement of rubredoxin marks a significant milestone, presenting the first HAR application of a metalloprotein, and further underlines the relevance of fragHAR in protein crystallography.

Particulate methane monooxygenase (pMMO) is the most efficient of the two groups of enzymes that can hydroxylate methane. The enzyme is membrane bound and therefore hard to study experimentally. For that reason, there is still no consensus regarding the location and nature of the active site. We have used combined quantum mechanical and molecular mechanical (QM/MM) methods to study the reactivity of the CuB site with a histidine brace and two additional histidine ligands. We compare it with the similar active site of lytic polysaccharide monooxygenases. We show that the CuB site can form a reactive [CuO]+ state by the addition of three electrons and two protons, starting from a resting Cu(ii) state, with a maximum barrier of 72 kJ mol-1. The [CuO]+ state can abstract a proton from methane, forming a Cu-bound OH- group, which may then recombine with the CH3 group, forming the methanol product. The two steps have barriers of 59 and 52 kJ mol-1, respectively. However, in many of the steps, formation and dissociation of H2O2 or HO2- compete with the formation of the [CuO]+ state and the former steps are typically more favourable. Thus, our calculations indicate that the CuB site is not employed for methane oxidation, but may rather be used for the formation of hydrogen peroxide. This conclusion concurs with recent experimental investigations that excludes the CuB site as the site for methane oxidation.

Density functional theory (DFT) thermochemistry of 3d transition-metal complexes is well-known to be sensitive to the amount of exact Hartree–Fock exchange incorporated into the exchange–correlation functional. For example, relative energies of different protonation states of iron–sulfur complexes may vary by hundreds of kJ/mol among different DFT methods. In the present study, we examine the relative energies of four protonation isomers of the [CH3S4Fe2IIIS2H]− [2Fe–2S] ferredoxin model. Compared to many-body ab initio phaseless auxiliary-field quantum Monte Carlo with multi-Slater determinant trial wavefunctions and fully connected singles and doubles coupled-cluster with perturbative triples methods, the r2SCAN12-D4, B3LYP-D4, and B97-1-D3(OP) approaches perform the best. We also demonstrate that density-corrected DFT on top of KS-CCSD electronic densities provides reliable results with the r2SCAN functional. Moreover, the direct random phase approximation on top of the TPSSh, O3LYP, and r2SCAN12 hybrid functionals performs well.

Lytic polysaccharide monooxygenases (LPMOs) are unique mono-copper enzymes that boost the degradation of different polysaccharides and play important roles in the sustainable production of biofuels, in human and plant pathogens and potentially also in plastic degradation. Their activity depends on a co-substrate, where recent results show that hydrogen peroxide is the preferred co-substrate. Under typical experimental conditions, no hydrogen peroxide is added and it is instead produced in situ by LPMOs themselves, which could possibly be the rate-limiting step. Previous theoretical investigations of the oxidase reaction have been highly inhomogeneous, and focused on different aspects of LPMO reactivity. In this paper, we systematically investigate how LPMOs generate hydrogen peroxide using accurate quantum mechanics/molecular mechanics (QM/MM) hybrid methods with extended QM regions. We find that protonation of a generated superoxide intermediate at the active-site is most likely, but that this requires further reduction of the superoxide.

Alchemical free-energy perturbation (FEP) is an accurate and thermodynamically stringent way to estimate relative energies for the binding of small ligands to biological macromolecules. It has repeatedly been pointed out that a single simulation normally stays near the starting point in phase space and therefore underestimates the uncertainty of the results. Therefore, it is better to run an ensemble of independent simulations. Traditionally, such an ensemble has been generated by using different starting velocities. We argue that it is better to use also other random choices made during the setup of the simulations, in particular the solvation of the solute. We show here that such solvent-induced independent simulations (SIS) sometimes give a larger standard deviation and slightly different results for the binding of 42 ligands to five different proteins, viz. human N-terminal bromodomain 4, the Leu99Ala mutant of T4 lysozyme, dihydrofolate reductase, blood-clotting factor Xa, and ferritin. SIS does not involve any increase in the time consumption. Therefore, we strongly recommend the use of SIS (in addition to different velocities) to start independent simulations. Other random or uncertain choices in the setup of the simulated systems, e.g., the selection of residues with alternative conformations or positions of added protons, may also be used to enhance the variation in independent simulations.

Particulate methane monooxygenase (pMMO) is an enzyme that converts methane into methanol at ambient temperature and pressure. Over the past three decades, the metal content and location of the active site have been highly controversial. Recent single-particle cryogenic electron-microscopy (cryo-EM) structures have furthered this debate. In this study, three cryo-EM structures (PDB entries 7s4h, 7s4j and 7ev9) are analysed by quantum refinement (QR). This approach augments traditional structural refinement with quantum-mechanical (QM) calculations for a small but interesting part of the protein (in this case, the copper sites). Our results indicate that the bis-His (CuA) site is correctly modelled as a mononuclear copper site in all three structures. The His-brace (CuB) site is also best modelled as mononuclear in all structures, although it was suggested to be a binuclear site in PDB entry 7ev9. The CuC site, which is observed only in PDB entry 7s4j, is correctly modelled and is probably reduced in the structure. The CuD putative active site, observed only in PDB entry 7s4h, is also mononuclear, but a water molecule might at least intermittently coordinate to the copper ion. On the other hand, our study does not find any support for the five additional copper ions suggested to be present in PDB entry 7ev9, including the suggested trinuclear active site and two sites in the so-called copper sponge. Instead, more chemically reasonable structures and better fit to both the cryo-EM and QM data are obtained if these copper ions are replaced with water molecules. This study illustrates the potential of QR as a standard component of cryo-EM studies for metal sites, for which reliable empirical restraints are missing.

Combined quantum mechanics and molecular mechanics (QM/MM) calculations are a popular approach to study reaction mechanisms of enzymes. However, recently, the reproducibility of such calculations has been questioned, comparing the results of two software: NWChem and Q-Chem. Here, we continue and extend this study by including three additional software─ComQum, ORCA, and AMBER─using the same test case, the covalent attachment of the carmofur inhibitor to the catalytic Cys-145 residue of the SARS-CoV-2 main protease, using a quantum region of 83 atoms. We confirm that the various software programs give varying results for the reaction (ΔE) and activation (ΔE‡) energies. The main reason for the variation is how charges around the cleaved bonds between the QM and MM regions are treated, i.e., the charge-redistribution scheme. However, there are still differences of ∼10 kJ/mol between different implementations of the same method in ComQum and ORCA. Some of these problems can be solved by calculating the final energies with larger QM systems. We show that energies calculated with the big-QM approach are reasonably converged if atoms within 8 Å of the minimal QM region are included (∼1400 atoms), solvent-exposed charged residues are neutralized, and the calculation is performed in a continuum solvent with a dielectric constant of 80. On the other hand, we show that different setups of the protein lead to even larger differences in the calculated energies, by up to 114 kJ/mol. Even if the same approach is used and the only difference is how water molecules are added (by random) to the crystal structure, energies differ by 18-57 kJ/mol. The results also strongly depend on how much of the surrounding protein and solvent are relaxed in the calculations. Therefore, it seems that for a solvent-exposed active site, QM/MM calculations with minimized structures cannot be recommended. Instead, methods that incorporate dynamic effects and calculate free energies seem preferable.

β-Alanine synthase (βAS), which is a dizinc metalloenzyme, catalyzes the irreversible hydrolysis of N-carbamyl-β-alanine (NCβA) to β-alanine. This enzyme has potential applications for β-amino acid production. Understanding the reaction mechanism and selectivity of βAS at atomic details can help design and engineer the enzyme for cascade biocatalysis. Here, the protonation states of two conserved active-site histidine residues (His262 and His397 in Saccharomyces kluyveri) of βAS were investigated by means of combined quantum mechanical and molecular mechanical (QM/MM) molecular dynamics (MD) simulation, as well as the ONIOM QM/QM' approach. The calculations predicted that both His262 and His397 should be neutral for efficient catalysis. Furthermore, the βAS reaction mechanism and its stereospecificity toward a series of NCβA substrates containing different β2 and β3-β-alanine substitutions were studied, which suggested factors governing the origin of stereoselectivity of this enzyme. The mechanism for the conversion of NCβA into β-alanine, carbon dioxide, and ammonia by βAS involved four reaction steps: nucleophilic attack by a hydroxide ion, substrate protonation and formation of a zwitterionic intermediate, and C-N bond cleavage to produce β-alanine and carbamate, which is finally decomposed into carbon dioxide and ammonia. The rate-limiting step is the protonation of the amide nitrogen of the substrate by Glu159, with the overall reaction barrier (16.5 kcal/mol) consistent with the experimental data. In silico alanine scanning analysis of the reaction mechanism for four variants (His262Ala, His397Ala, Asn309Ala, and Arg322Ala) is performed, showing increased activation energies compared to the wild-type enzyme, which confirms the roles of these residues in catalysis. The results explain the enzyme's preference for linear N-carbamyl substrates, as large and branched substrates cannot fit in the active site, restricted by the residue of the loop/region of the enzyme. Overall, we have demonstrated that a combined use of QM/MM MD and ONIOM models can be a promising strategy to elucidate possible protonation states of the ionizable residues in the enzyme active site prior to catalysis.

Solvent cage-escape dynamics of bimolecular photoredox products in solution has been investigated computationally through a combination of molecular dynamics simulations and quantum chemical calculations. The present work focuses on the photoinduced oxidation of the organic electron donor dimethylaniline (DMA) by a Fe(III) N-heterocyclic carbene photosensitizer (Fe(III)NHC+) in two different solvents, serving as an example of current interest due to their relevance for the development of earth-abundant photocatalytic systems. Calculated solvent cage-escape yields of radical-cation and neutral photoproducts (DMA•+ and Fe(II)NHC, respectively) by molecular dynamics simulations reveal more favorable solvation in acetonitrile than in dichloromethane following the initial photoinduced charge-separation. These results agree with basic expectations from solvent polarity considerations but give an opposite trend compared to experimentally reported cage-escape yields. Alternative cage-escape mechanisms were therefore considered computationally to account for the anomalous experimental cage-escape yields. Both quantum chemical calculations and molecular dynamics simulations support the formation of radical-cation dimers ([(DMA)2]•+), allowing for more efficient charge migration involving the radical-cations as the polarity of the solvent is decreased. The results further demonstrate the ability of the counterion (PF6-) to stabilize the photoproducts through radical-cation-anion pairing, suggesting that these bimolecular interactions can also play an important role to preferentially promote photoproduct formation in less polar solvents. Both radical-cation dimer formation and radical-cation-counterion interactions are therefore proposed to provide additional pathways that help to explain the experimental observations of anomalous solvation dependence of the cage-escape dynamics in the investigated system. The broader implications of the bimolecular cage-escape processes on photocatalytic reaction dynamics are also considered based on our findings about light-induced intermolecular interactions.

The catalytic activity of the binuclear glyoxalase II (GlxII) enzyme is closely linked to the type and charge of metal ions in its active site. Using hybrid quantum mechanics/molecular mechanics (QM/MM) calculations, we investigated the reaction mechanism of human GlxII, which features two Zn(II) ions in its active site. By systematically replacing these Zn(II) ions with Fe(II), Fe(III), or Co(II), we evaluated the impact of metal substitutions on reaction energetics and active-site geometry. Our results reveal that the type and position of the metal ions are critical to the catalytic activity of GlxII. Substitution of the Zn(II) ion in the three-histidine site with Fe(II), Fe(III), or Co(II) significantly increased the activation barrier, indicating that these configurations are less favorable. In contrast, substituting Zn(II) in the two-histidine site with either Fe(II) or Co(II) resulted in a reduced activation barrier and produced geometries closely resembling those observed when both metal sites are occupied by Zn(II). Additionally, moving the metal ions from the QM to the MM region inhibited the reaction, highlighting their direct chemical involvement in catalysis beyond electrostatic stabilization. These results underscore that the metal ions chemically participate in the catalytic process beyond their electrostatic contributions. Collectively, our results provide insights into the structural and electronic factors governing GlxII catalysis, offering a theoretical framework to complement and refine experimental studies.

Uronate isomerase (EC 5.3.1.12; URI) catalyzes uronate sugar interconversion, a key step in bacterial metabolism, yet its reaction mechanism remains poorly understood. This study delineates the detailed mechanism for the isomerization of d-glucuronate to d-fructuronate catalyzed by a Zn2+-dependent URI enzyme (from Bacillus halodurans). Using quantum mechanical (QM) cluster calculations, three mechanistic pathways were evaluated, all involving Asp355 as the catalytic base but differing in the proton-shuttle mechanism. The most favorable mechanism features C5 deprotonation of the substrate, followed by a water-mediated 1,2-proton transfer via a stabilized cis-enediol intermediate, with the C2-C5 intramolecular transfer as the rate-determining step. The calculated activation barrier (17.3 kcal mol- 1) aligns well with experimental data. Alternative pathways involving Tyr50 or a Tyr48-water relay were found to be less favorable due to higher net barriers (26-36 kcal mol- 1). A comparative analysis of d-glucuronate and d-galacturonate revealed that although both proceed through similar steps, d-glucuronate has ∼2 kcal mol- 1 lower overall barrier due to enhanced stabilization of late-stage intermediates. These findings clarify the roles of solvent and active-site residues in URI catalysis and contribute to a broader understanding of proton-transfer mechanisms in Zn2+-dependent enzymes across the amidohydrolase superfamily.

Lytic polysaccharide monooxygenases (LPMOs) are copper-dependent enzymes that have fueled the hope for sustainable biofuel production since they enhance the breakdown of recalcitrant polysaccharides like cellulose. In the consensus mechanism, their catalytic activity relies on forming an 'oxyl', [CuO center dot-]+, species at the active site, followed by subsequent hydrogen atom abstraction (HAA) from the substrate. Some studies report rather high barriers for this reaction, identifying it as the rate-limiting step in the oxidation process, whereas other investigations have reported significantly lower barriers. In this study, we have constructed a force field for the active site and show through extensive sampling from molecular dynamics simulations that the QM/MM reaction barrier depends critically on the underlying structural conformations of the enzyme-substrate complex. The results support low-energy barriers for the HAA step and help to explain previous discrepancies in the literature, which may be attributed to insufficient conformational sampling.