Interplay of halogen bonding and solvation in protein-ligand binding

Halogen bonding is increasingly utilized in efforts to achieve high affinity and selectivity of molecules designed to bind proteins, making it paramount to understand the relationship between structure, dynamics, and thermodynamic driving forces. We present a detailed analysis addressing this problem using a series of protein-ligand complexes involving single halogen substitutions - F, Cl, Br, and I - and nearly identical structures. Isothermal titration calorimetry reveals an increasingly favorable binding enthalpy from F to I that correlates with the halogen size and σ-hole electropositive character, but is partially counteracted by unfavorable entropy, which is constant from F to Cl and Br, but worse for I. Consequently, the binding free energy is roughly equal for Cl, Br, and I. QM and solvation-free-energy calculations reflect an intricate balance between halogen bonding, hydrogen bonds, and solvation. These advances have the potential to aid future drug design initiatives involving halogenated compounds.

CO Oxidation Mechanism of Silver-Substituted Mo/Cu CO-Dehydrogenase - Analogies and Differences to the Native Enzyme

The aerobic oxidation of carbon monoxide to carbon dioxide is catalysed by the Mo/Cu-containing CO-dehydrogenase enzyme in the soil bacterium Oligotropha carboxidovorans, enabling the organism to grow on the small gas molecule as carbon and energy source. It was shown experimentally that silver can be substituted for copper in the active site of Mo/Cu CODH to yield a functional enzyme. In this study, we employed QM/MM calculations to investigate whether the reaction mechanism of the silver-substituted enzyme is similar to that of the native enzyme. Our results suggest that the Ag-substituted enzyme can oxidize CO and release CO2 following the same reaction steps as the native enzyme, with a computed rate-limiting step of 10.4 kcal/mol, consistent with experimental findings. Surprisingly, lower activation energies for C-O bond formation have been found in the presence of silver. Furthermore, comparison of rate constants for reduction of copper- and silver-containing enzymes suggests a discrepancy in the transition state stabilization upon silver substitution. We also evaluated the effects that differences in the water-active site interaction may exert on the overall energy profile of catalysis. Finally, the formation of a thiocarbonate intermediate along the catalytic pathway was found to be energetically unfavorable for the Ag-substituted enzyme. This finding aligns with the hypothesis proposed for the wild-type form, suggesting that the creation of such species may not be necessary for the enzymatic catalysis of CO oxidation.

QM/MM study of the catalytic reaction of aphid myrosinase

Brevicoryne brassicae, an aphid species, exclusively consumes plants from the Brassicaceae family and employs a sophisticated defense mechanism involving a myrosinase enzyme that breaks down glucosinolates obtained from its host plants. In this work, we employed combined quantum mechanical and molecular mechanical (QM/MM) calculations and molecular dynamics (MD) simulations to study the catalytic reaction of aphid myrosinase. A proper QM region to study the myrosinase reaction should contain the whole substrate, models of Gln-19, His-122, Asp-124, Asn-166, Glu-167, Lys-173, Tyr-180, Val-228, Tyr-309, Tyr-346, Ile-347, Glu-374, Glu-423, Trp-424, and a water molecule. The calculations show that Asp-124 and Glu-423 must be charged, His-122 must be protonated on NE2, and Glu-167 must be protonated on OE2. Our model reproduces the anomeric retaining characteristic of myrosinase and indicates that the deglycosylation reaction is the rate-determining step of the reaction. Based on the calculations, we propose a reaction mechanism for aphid myrosinase-mediated hydrolysis of glucosinolates with an overall barrier of 15.2 kcal/mol. According to the results, removing a proton from Arg-312 or altering it to valine or methionine increases glycosylation barriers but decreases the deglycosylation barrier.

Scalar Relativistic All-Electron and Pseudopotential Ab Initio Study of a Minimal Nitrogenase [Fe(SH)4H]- Model Employing Coupled-Cluster and Auxiliary-Field Quantum Monte Carlo Many-Body Methods

Nitrogenase is the only enzyme that can cleave the triple bond in N2, making nitrogen available to organisms. The detailed mechanism of this enzyme is currently not known, and computational studies are complicated by the fact that different density functional theory (DFT) methods give very different energetic results for calculations involving nitrogenase models. Recently, we designed a [Fe(SH)4H]- model with the fifth proton binding either to Fe or S to mimic different possible protonation states of the nitrogenase active site. We showed that the energy difference between these two isomers (ΔE) is hard to estimate with quantum-mechanical methods. Based on nonrelativistic single-reference coupled-cluster (CC) calculations, we estimated that the ΔE is 101 kJ/mol. In this study, we demonstrate that scalar relativistic effects play an important role and significantly affect ΔE. Our best revised single-reference CC estimates for ΔE are 85-91 kJ/mol, including energy corrections to account for contributions beyond triples, core-valence correlation, and basis-set incompleteness error. Among coupled-cluster approaches with approximate triples, the canonical CCSD(T) exhibits the largest error for this problem. Complementary to CC, we also used phaseless auxiliary-field quantum Monte Carlo calculations (ph-AFQMC). We show that with a Hartree-Fock (HF) trial wave function, ph-AFQMC reproduces the CC results within 5 ± 1 kJ/mol. With multi-Slater-determinant (MSD) trials, the results are 82-84 ± 2 kJ/mol, indicating that multireference effects may be rather modest. Among the DFT methods tested, τ-HCTH, r2SCAN with 10-13% HF exchange with and without dispersion, and O3LYP/O3LYP-D4, and B3LYP*/B3LYP*-D4 generally perform the best. The r2SCAN12 (with 12% HF exchange) functional mimics both the best reference MSD ph-AFQMC and CC ΔE results within 2 kJ/mol.

H2 formation from the E2-E4 states of nitrogenase

Nitrogenase is the only enzyme that can cleave the strong triple bond in N2, making nitrogen available for biological lifeforms. The active site is a MoFe7S9C cluster (the FeMo cluster) that binds eight electrons and protons during one catalytic cycle, giving rise to eight intermediate states E0-E7. It is experimentally known that N2 binds to the E4 state and that H2 is a compulsory byproduct of the reaction. However, formation of H2 is also an unproductive side reaction that should be avoided, especially in the early steps of the reaction mechanism (E2 and E3). Here, we study the formation of H2 for various structural interpretations of the E2-E4 states using combined quantum mechanical and molecular mechanical (QM/MM) calculations and four different density-functional theory methods. We find large differences in the predictions of the different methods. B3LYP strongly favours protonation of the central carbide ion and H2 cannot form from such structures. On the other hand, with TPSS, r2SCAN and TPSSh, H2 formation is strongly exothermic for all structures and En and therefore need strict kinetic control to be avoided. For the E2 state, the kinetic barriers for the low-energy structures are high enough to avoid H2 formation. However, for both the E3 and E4 states, all three methods predict that the best structure has two hydride ions bridging the same pair of Fe ions (Fe2 and Fe6) and these two ions can combine to form H2 with an activation barrier of only 29-57 kJ mol-1, corresponding to rates of 7 × 102 to 5 × 107 s-1, i.e. much faster than the turnover rate of the enzyme (1-5 s-1). We have also studied H-atom movements within the FeMo cluster, showing that the various protonation states can quite freely be interconverted (activation barriers of 12-69 kJ mol-1).

Protonation of Homocitrate and the E1 State of Fe-Nitrogenase Studied by QM/MM Calculations

Nitrogenase is the only enzyme that can cleave the strong triple bond in N2, making nitrogen available for biological life. There are three isozymes of nitrogenase, differing in the composition of the active site, viz., Mo, V, and Fe-nitrogenase. Recently, the first crystal structure of Fe-nitrogenase was presented. We have performed the first combined quantum mechanical and molecular mechanical (QM/MM) study of Fe-nitrogenase. We show with QM/MM and quantum-refinement calculations that the homocitrate ligand is most likely protonated on the alcohol oxygen in the resting E0 state. The most stable broken-symmetry (BS) states are the same as for Mo-nitrogenase, i.e., the three Noodleman BS7-type states (with a surplus of β spin on the eighth Fe ion), which maximize the number of nearby antiferromagnetically coupled Fe-Fe pairs. For the E1 state, we find that protonation of the S2B μ2 belt sulfide ion is most favorable, 14-117 kJ/mol more stable than structures with a Fe-bound hydride ion (the best has a hydride ion on the Fe2 ion) calculated with four different density-functional theory methods. This is similar to what was found for Mo-nitrogenase, but it does not explain the recent EPR observation that the E1 state of Fe-nitrogenase should contain a photolyzable hydride ion. For the E1 state, many BS states are close in energy, and the preferred BS state differs depending on the position of the extra proton and which density functional is used.

Multireference Protonation Energetics of a Dimeric Model of Nitrogenase Iron-Sulfur Clusters

Characterizing the electronic structure of the iron-sulfur clusters in nitrogenase is necessary to understand their role in the nitrogen fixation process. One challenging task is to determine the protonation state of the intermediates in the nitrogen fixing cycle. Here, we use a dimeric iron-sulfur model to study relative energies of protonation at C, S, or Fe. Using a composite method based on coupled cluster and density matrix renormalization group energetics, we converge the relative energies of four protonated configurations with respect to basis set and correlation level. We find that accurate relative energies require large basis sets as well as a proper treatment of multireference and relativistic effects. We have also tested ten density functional approximations for these systems. Most of them give large errors in their relative energies. The best performing functional in this system is B3LYP, which gives mean absolute and maximum deviations of only 10 and 13 kJ/mol with respect to our correlated wave function estimates, respectively, comparable to the uncertainty in our correlated estimates. Our work provides benchmark results for the calibration of new approximate electronic structure methods and density functionals for these problems.

Assessment of DFT functionals for a minimal nitrogenase [Fe(SH)4H]- model employing state-of-the-art ab initio methods

We have designed a [Fe(SH)4H]- model with the fifth proton binding either to Fe or S. We show that the energy difference between these two isomers (∆E) is hard to estimate with quantum-mechanical (QM) methods. For example, different density functional theory (DFT) methods give ∆E estimates that vary by almost 140 kJ/mol, mainly depending on the amount of exact Hartree-Fock included (0%-54%). The model is so small that it can be treated by many high-level QM methods, including coupled-cluster (CC) and multiconfigurational perturbation theory approaches. With extrapolated CC series (up to fully connected coupled-cluster calculations with singles, doubles, and triples) and semistochastic heat-bath configuration interaction methods, we obtain results that seem to be converged to full configuration interaction results within 5 kJ/mol. Our best result for ∆E is 101 kJ/mol. With this reference, we show that M06 and B3LYP-D3 give the best results among 35 DFT methods tested for this system. Brueckner doubles coupled cluster with perturbaitve triples seems to be the most accurate coupled-cluster approach with approximate triples. CCSD(T) with Kohn-Sham orbitals gives results within 4-11 kJ/mol of the extrapolated CC results, depending on the DFT method. Single-reference CC calculations seem to be reasonably accurate (giving an error of ∼5 kJ/mol compared to multireference methods), even if the D1 diagnostic is quite high (0.25) for one of the two isomers.

Two local minima for structures of [4Fe-4S] clusters obtained with density functional theory methods

[4Fe-4S] clusters are essential cofactors in many proteins involved in biological redox-active processes. Density functional theory (DFT) methods are widely used to study these clusters. Previous investigations have indicated that there exist two local minima for these clusters in proteins. We perform a detailed study of these minima in five proteins and two oxidation states, using combined quantum mechanical and molecular mechanical (QM/MM) methods. We show that one local minimum (L state) has longer Fe-Fe distances than the other (S state), and that the L state is more stable for all cases studied. We also show that some DFT methods may only obtain the L state, while others may obtain both states. Our work provides new insights into the structural diversity and stability of [4Fe-4S] clusters in proteins, and highlights the importance of reliable DFT methods and geometry optimization. We recommend r2SCAN for optimizing [4Fe-4S] clusters in proteins, which gives the most accurate structures for the five proteins studied.

Comparison of force fields to study the zinc-finger containing protein NPL4, a target for disulfiram in cancer therapy

Molecular dynamics (MD) simulations are a powerful approach to studying the structure and dynamics of proteins related to health and disease. Advances in the MD field allow modeling proteins with high accuracy. However, modeling metal ions and their interactions with proteins is still challenging. NPL4 is a zinc-binding protein and works as a cofactor for p97 to regulate protein homeostasis. NPL4 is of biomedical importance and has been proposed as the target of disulfiram, a drug recently repurposed for cancer treatment. Experimental studies proposed that the disulfiram metabolites, bis-(diethyldithiocarbamate)‑copper and cupric ions, induce NPL4 misfolding and aggregation. However, the molecular details of their interactions with NPL4 and consequent structural effects are still elusive. Here, biomolecular simulations can help to shed light on the related structural details. To apply MD simulations to NPL4 and its interaction with copper the first important step is identifying a suitable force field to describe the protein in its zinc-bound states. We examined different sets of non-bonded parameters because we want to study the misfolding mechanism and cannot rule out that the zinc may detach from the protein during the process and copper replaces it. We investigated the force-field ability to model the coordination geometry of the metal ions by comparing the results from MD simulations with optimized geometries from quantum mechanics (QM) calculations using model systems of NPL4. Furthermore, we investigated the performance of a force field including bonded parameters to treat copper ions in NPL4 that we obtained based on QM calculations.

N2 binding to the E0-E4 states of nitrogenase

Nitrogenase is the only enzyme that can convert N₂ into NH₃. The reaction requires the addition of eight electrons and protons to the enzyme and the mechanism is normally described by nine states, E₀-E₈, differing in the number of added electrons. Experimentally, it is known that three or four electrons need to be added before the enzyme can bind N₂. We have used combined quantum mechanical and molecular mechanics methods to study the binding of N₂ to the E₀-E₄ states of nitrogenase, using four different density functional theory (DFT) methods. We test many different structures for the E₂-E₄ states and study binding both to the Fe2 and Fe6 ions of the active-site FeMo cluster. Unfortunately, the results depend quite strongly on the DFT methods. The TPSS method gives the strongest bonding and prefers N₂ binding to Fe6. It is the only method that reproduces the experimental observation of unfavourable binding to the E₀-E₂ states and favourable binding to E₃ and E₄. The other three methods give weaker binding, preferably to Fe2. B3LYP strongly favours structures with the central carbide ion triply protonated. The other three methods suggest that states with the S2B ligand dissociated from either Fe2 or Fe6 are competitive for the E₂-E₄ states. Moreover, such structures with two hydride ions both bridging Fe2 and Fe6 are the best models for E₄ and also for the N₂-bound E₃ and E₄ states. However, for E₄, other structures are often close in energy, e.g. structures with one of the hydride ions bridging instead Fe3 and Fe7. Finally, we find no support for the suggestion that reductive elimination of H₂ from the two bridging hydride ions in the E₄ state would enhance the binding of N₂.

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.

Histidine oxidation in lytic polysaccharide monooxygenase

The lytic polysaccharide monooxygenases (LPMOs) comprise a super-family of copper enzymes that boost the depolymerisation of polysaccharides by oxidatively disrupting the glycosidic bonds connecting the sugar units. Industrial use of LPMOs for cellulose depolymerisation has already begun but is still far from reaching its full potential. One issue is that the LPMOs self-oxidise and thereby deactivate. The mechanism of this self-oxidation is unknown, but histidine residues coordinating to the copper atom are the most susceptible. An unusual methyl modification of the NE2 atom in one of the coordinating histidine residues has been proposed to have a protective role. Furthermore, substrate binding is also known to reduce oxidative damage. We here for the first time investigate the mechanism of histidine oxidation with combined quantum and molecular mechanical (QM/MM) calculations, with outset in intermediates previously shown to form from a reaction with peroxide and a reduced LPMO. We show that an intermediate with a [Cu-O]⁺ moiety is sufficiently potent to oxidise the nearest C-H bond on both histidine residues, but methylation of the NE2 atom of His-1 increases the reaction barrier of this reaction. The substrate further increases the activation barrier. We also investigate a [Cu-OH]²⁺ intermediate with a deprotonated tyrosine radical. This intermediate was previously proposed to have a protective role, and we also find it to have higher barriers than the corresponding a [Cu-O]⁺ intermediate.

A computational study of the reaction mechanism and stereospecificity of dihydropyrimidinase

Dihydropyrimidinase (DHPase) is a key enzyme in the pyrimidine pathway, the catabolic route for synthesis of β-amino acids. It catalyses the reversible conversion of 5,6-dihydrouracil (DHU) or 5,6-dihydrothymine (DHT) to the corresponding N-carbamoyl-β-amino acids. This enzyme has the potential to be used as a tool in the production of β-amino acids. Here, the reaction mechanism and origin of stereospecificity of DHPases from Saccharomyces kluyveri and Sinorhizobium meliloti CECT4114 were investigated and compared using a quantum mechanical cluster approach based on density functional theory. Two models of the enzyme active site were designed from the X-ray crystal structure of the native enzyme: a small cluster to characterize the mechanism and the stationary points and a large model to probe the stereospecificity and the role of stereo-gate-loop (SGL) residues. It is shown that a hydroxide ion first performs a nucleophilic attack on the substrate, followed by the abstraction of a proton by Asp358, which occurs concertedly with protonation of the ring nitrogen by the same residue. For the DHT substrate, the enzyme displays a preference for the L-configuration, in good agreement with experimental observation. Comparison of the reaction energetics of the two models reveals the importance of SGL residues in the stereospecificity of catalysis. The role of the conserved Tyr172 residue in transition-state stabilization is confirmed as the Tyr172Phe mutation increases the activation barrier of the reaction by ∼8 kcal mol^-1. A detailed understanding of the catalytic mechanism of the enzyme could offer insight for engineering in order to enhance its activity and substrate scope.

How general is the effect of the bulkiness of organic ligands on the basicity of metal-organic catalysts? H2-evolving Mo oxides/sulphides as case studies

Tailoring the activity of an organometallic catalyst usually requires a targeted ligand design. Tuning the ligand bulkiness and tuning the electronic properties are popular approaches, which are somehow interdependent because substituents of different sizes within ligands can determine inter alia the occurrence of different degrees of inductive effects. Ligand basicity, in particular, turned out to be a key property for the modulation of protonation reactions occurring in vacuo at the metals in complexes bearing organophosphorus ligands; however, when the same reactions take place in a polar organic solvent, their energetics becomes dependent on the trade-off between ligand basicity and bulkiness, with the polarity of the solvent playing a key role in this regard [Bancroft et al., Inorg. Chem., 1986, 25, 3675; Rovaletti et al., J. Phys. Org. Chem., 2018, 31, e3748]. In the present contribution, we carried out molecular dynamics and density functional theory calculations on water-soluble Mo-based catalysts for proton reduction, in order to study the energetics of protonation reactions in complexes where the incipient proton binds a catalytically active ligand (i.e., an oxide or a disulphide). We considered complexes either soaked in water or in a vacuum, and featuring N-based ancillary ligands of different bulkiness (i.e. cages constituted either by pyridine or isoquinoline moieties). Our results show that the energetics of protonation events can be affected by ancillary ligand bulkiness even when the metal center does not play the role of the H⁺ acceptor. In vacuo, protonation at the O or S atom in the α position relative to the metal in complexes featuring the bulky isoquinoline-based ligand is more favored by around 10 kcal mol^-1 when compared to the case of the pyridine-based counterparts, a difference that is almost zero when the same reactions occur in water. Such an outcome is rationalized in light of the different electrostatic properties of complexes bearing ancillary ligands of different sizes. The overall picture from theory indicates that such effects of ligand bulkiness can be relevant for the design of green chemistry catalysts that undergo protonation steps in water solutions.

Proton Transfer Pathways in Nitrogenase with and without Dissociated S2B

Nitrogenase is the only enzyme that can convert N₂ to NH₃ . Crystallographic structures have indicated that one of the sulfide ligands of the active-site FeMo cluster, S2B, can be replaced by an inhibitor, like CO and OH^- , and it has been suggested that it may be displaced also during the normal reaction. We have investigated possible proton transfer pathways within the FeMo cluster during the conversion of N₂ H₂ to two molecules of NH₃ , assuming that the protons enter the cluster at the S3B, S4B or S5A sulfide ions and are then transferred to the substrate. We use combined quantum mechanical and molecular mechanical (QM/MM) calculations with the TPSS and B3LYP functionals. The calculations indicate that the barriers for these reactions are reasonable if the S2B ligand remains bound to the cluster, but they become prohibitively high if S2B has dissociated. This suggests that it is unlikely that S2B reversibly dissociates during the normal reaction cycle.

Computationally enhanced X-ray diffraction analysis of a gold(III) complex interacting with the human telomeric DNA G-quadruplex. Unravelling non-unique ligand positioning

The crystal structure of the human telomeric DNA Tel24 G-quadruplex (Tel24 = TAG₃(T₂AG₃)₃T) in complex with the novel [AuL] species (with L = 2,4,6-tris(2-pyrimidyl)-1,3,5-triazine - TPymT-α) was solved by a novel joint molecular mechanical (MM)/quantum mechanical (QM) innovative approach. The quantum-refinement crystallographic method (crystallographic refinement enhanced with quantum mechanical calculation) was adapted to treat the [AuL]/G-quadruplex structure, where each gold complex in the binding site was found spread over four equally occupied positions. The four positions were first determined by docking restrained to the crystallographically determined metal ions' coordinates. Then, the quantum refinement method was used to resolve the poorly defined density around the ligands and improve the crystallographic determination, revealing that the binding preferences of this metallodrug toward Tel24 G-quadruplex arise from a combined effect of pyrimidine stacking, metal-guanine interactions and charge-charge neutralizing action of the π-acid triazine. The occurrence of interaction in solution with the Tel24 G-quadruplex DNA was further proved through DNA melting experiments, which showed a slight destabilisation of the quadruplex upon adduct formation.

Comparison of Grand Canonical and Conventional Molecular Dynamics Simulation Methods for Protein-Bound Water Networks

Water molecules play important roles in all biochemical processes. Therefore, it is of key importance to obtain information of the structure, dynamics, and thermodynamics of water molecules around proteins. Numerous computational methods have been suggested with this aim. In this study, we compare the performance of conventional and grand-canonical Monte Carlo (GCMC) molecular dynamics (MD) simulations to sample the water structure, as well GCMC and grid-based inhomogeneous solvation theory (GIST) to describe the energetics of the water network. They are evaluated on two proteins: the buried ligand-binding site of a ferritin dimer and the solvent-exposed binding site of galectin-3. We show that GCMC/MD simulations significantly speed up the sampling and equilibration of water molecules in the buried binding site, thereby making the results more similar for simulations started from different states. Both GCMC/MD and conventional MD reproduce crystal-water molecules reasonably for the buried binding site. GIST analyses are normally based on restrained MD simulations. This improves the precision of the calculated energies, but the restraints also significantly affect both absolute and relative energies. Solvation free energies for individual water molecules calculated with and without restraints show a good correlation, but with large quantitative differences. Finally, we note that the solvation free energies calculated with GIST are ∼5 times larger than those estimated by GCMC owing to differences in the reference state.

Benchmark Study of Redox Potential Calculations for Iron-Sulfur Clusters in Proteins

Redox potentials have been calculated for 12 different iron–sulfur sites of 6 different types with 1–4 iron ions. Structures were optimized with combined quantum mechanical and molecular mechanical (QM/MM) methods, and the redox potentials were calculated using the QM/MM energies, single-point QM methods in a continuum solvent or by QM/MM thermodynamic cycle perturbations. We show that the best results are obtained with a large QM system (∼300 atoms, but a smaller QM system, ∼150 atoms, can be used for the QM/MM geometry optimization) and a large value of the dielectric constant (80). For absolute redox potentials, the B3LYP density functional method gives better results than TPSS, and the results are improved with a larger basis set. However, for relative redox potentials, the opposite is true. The results are insensitive to the force field (charges of the surroundings) used for the QM/MM calculations or whether the protein and solvent outside the QM system are relaxed or kept fixed at the crystal structure. With the best approach for relative potentials, mean absolute and maximum deviations of 0.17 and 0.44 V, respectively, are obtained after removing a systematic error of −0.55 V. Such an approach can be used to identify the correct oxidation states involved in a certain redox reaction.

We have studied redox potentials of 12 iron−sulfur sites of 6 types with 1−4 iron ions. Structures were optimized with combined quantum mechanical and molecular mechanical (QM/MM) methods, and the redox potentials were calculated with QM/MM, QM calculations in a continuum solvent or by QM/MM thermodynamic cycle perturbations. The best results are obtained with the second approach using ∼300 atoms in the QM model and a large dielectric constant.

Thermodynamically Favourable States in the Reaction of Nitrogenase without Dissociation of any Sulfide Ligand

We have used combined quantum mechanical and molecular mechanical (QM/MM) calculations to study the reaction mechanism of nitrogenase, assuming that none of the sulfide ligands dissociates. To avoid the problem that there is no consensus regarding the structure and protonation of the E₄ state, we start from a state where N₂ is bound to the cluster and is protonated to N₂H₂, after dissociation of H₂. We show that the reaction follows an alternating mechanism with HNNH (possibly protonated to HNNH₂) and H₂NNH₂ as intermediates and the two NH₃ products dissociate at the E₇ and E₈ levels. For all intermediates, coordination to Fe6 is preferred, but for the E₄ and E₈ intermediates, binding to Fe2 is competitive. For the E₄, E₅ and E₇ intermediates we find that the substrate may abstract a proton from the hydroxy group of the homocitrate ligand of the FeMo cluster, thereby forming HNNH₂, H₂NNH₂ and NH₃ intermediates. This may explain why homocitrate is a mandatory component of nitrogenase. All steps in the suggested reaction mechanism are thermodynamically favourable compared to protonation of the nearby His‐195 group and in all cases, protonation of the NE2 atom of the latter group is preferred.

Can water act as a nucleophile in CO oxidation catalysed by Mo/Cu CO-dehydrogenase? Answers from theory

The aerobic CO dehydrogenase from Oligotropha carboxidovorans is an environmentally crucial bacterial enzyme for maintenance of subtoxic concentration of CO in the lower atmosphere, as it allows for the oxidation of CO to CO₂ which takes place at its Mo−Cu heterobimetallic active site. Despite extensive experimental and theoretical efforts, significant uncertainties still concern the reaction mechanism for the CO oxidation. In this work, we used the hybrid quantum mechanical/molecular mechanical approach to evaluate whether a water molecule present in the active site might act as a nucleophile upon formation of the new C−O bond, a hypothesis recently suggested in the literature. Our study shows that activation of H₂O can be favoured by the presence of the Mo=O_eq group. However, overall our results suggest that mechanisms other than the nucleophilic attack by Mo=O_eq to the activated carbon of the CO substrate are not likely to constitute reactive channels for the oxidation of CO by the enzyme.

Combining crystallography with quantum mechanics

In standard crystallographic refinement of biomacromolecules, the crystallographic raw data are supplemented by empirical restraints that ensure that the structure makes chemical sense. These restraints are typically accurate for amino acids and nucleic acids, but less so for cofactors, substrates, inhibitors, ligands and metal sites. In quantum refinement, this potential is replaced by more accurate quantum mechanical (QM) calculations. Several implementations have been presented, differing in the level of QM and whether it is used for the entire structure or only for a site of particular interest. It has been shown that the method can improve and correct errors in crystal structures and that it can be used to determine protonation and tautomeric states of various ligands and to decide what is really seen in the structure by refining different interpretations and using standard crystallographic and QM quality measures to decide which fits the structure best.

On the Use of Interaction Entropy and Related Methods to Estimate Binding Entropies

Molecular mechanics combined with Poisson-Boltzmann or generalized Born and solvent-accessible area solvation energies (MM/PBSA and MM/GBSA) are popular methods to estimate the free energy for the binding of small molecules to biomacromolecules. However, the estimation of the entropy has been problematic and time-consuming. Traditionally, normal-mode analysis has been used to estimate the entropy, but more recently, alternative approaches have been suggested. In particular, it has been suggested that exponential averaging of the electrostatic and Lennard-Jones interaction energies may provide much faster and more accurate entropies, the interaction entropy (IE) approach. In this study, we show that this exponential averaging is extremely poorly conditioned. Using stochastic simulations, assuming that the interaction energies follow a Gaussian distribution, we show that if the standard deviation of the interaction energies (σIE) is larger than 15 kJ/mol, it becomes practically impossible to converge the interaction entropies (more than 10 million energies are needed, and the number increases exponentially). A cumulant approximation to the second order of the exponential average shows a better convergence, but for σIE > 25 kJ/mol, it gives entropies that are unrealistically large. Moreover, in practical applications, both methods show a steady increase in the entropy with the number of energies considered.

Exploring ligand dynamics in protein crystal structures with ensemble refinement

Understanding the dynamics of ligands bound to proteins is an important task in medicinal chemistry and drug design. However, the dominant technique for determining protein-ligand structures, X-ray crystallography, does not fully account for dynamics and cannot accurately describe the movements of ligands in protein binding sites. In this article, an alternative method, ensemble refinement, is used on six protein-ligand complexes with the aim of understanding the conformational diversity of ligands in protein crystal structures. The results show that ensemble refinement sometimes indicates that the flexibility of parts of the ligand and some protein side chains is larger than that which can be described by a single conformation and atomic displacement parameters. However, since the electron-density maps are comparable and Rfree values are slightly increased, the original crystal structure is still a better model from a statistical point of view. On the other hand, it is shown that molecular-dynamics simulations and automatic generation of alternative conformations in crystallographic refinement confirm that the flexibility of these groups is larger than is observed in standard refinement. Moreover, the flexible groups in ensemble refinement coincide with groups that give high atomic displacement parameters or non-unity occupancy if optimized in standard refinement. Therefore, the conformational diversity indicated by ensemble refinement seems to be qualitatively correct, indicating that ensemble refinement can be an important complement to standard crystallographic refinement as a tool to discover which parts of crystal structures may show extensive flexibility and therefore are poorly described by a single conformation. However, the diversity of the ensembles is often exaggerated (probably partly owing to the rather poor force field employed) and the ensembles should not be trusted in detail.

Neutron structures of Leishmania mexicana triosephosphate isomerase in complex with reaction-intermediate mimics shed light on the proton-shuttling steps

Triosephosphate isomerase (TIM) is a key enzyme in glycolysis that catalyses the interconversion of glyceraldehyde 3-phosphate and dihydroxy-acetone phosphate. This simple reaction involves the shuttling of protons mediated by protolysable side chains. The catalytic power of TIM is thought to stem from its ability to facilitate the deprotonation of a carbon next to a carbonyl group to generate an enediolate intermediate. The enediolate intermediate is believed to be mimicked by the inhibitor 2-phosphoglycolate (PGA) and the subsequent enediol intermediate by phosphoglycolohydroxamate (PGH). Here, neutron structures of Leishmania mexicana TIM have been determined with both inhibitors, and joint neutron/X-ray refinement followed by quantum refinement has been performed. The structures show that in the PGA complex the postulated general base Glu167 is protonated, while in the PGH complex it remains deprotonated. The deuteron is clearly localized on Glu167 in the PGA-TIM structure, suggesting an asymmetric hydrogen bond instead of a low-barrier hydrogen bond. The full picture of the active-site protonation states allowed an investigation of the reaction mechanism using density-functional theory calculations.

Entropy-Entropy Compensation between the Protein, Ligand, and Solvent Degrees of Freedom Fine-Tunes Affinity in Ligand Binding to Galectin-3C

Molecular recognition is fundamental to biological signaling. A central question is how individual interactions between molecular moieties affect the thermodynamics of ligand binding to proteins and how these effects might propagate beyond the immediate neighborhood of the binding site. Here, we investigate this question by introducing minor changes in ligand structure and characterizing the effects of these on ligand affinity to the carbohydrate recognition domain of galectin-3, using a combination of isothermal titration calorimetry, X-ray crystallography, NMR relaxation, and computational approaches including molecular dynamics (MD) simulations and grid inhomogeneous solvation theory (GIST). We studied a congeneric series of ligands with a fluorophenyl-triazole moiety, where the fluorine substituent varies between the ortho, meta, and para positions (denoted O, M, and P). The M and P ligands have similar affinities, whereas the O ligand has 3-fold lower affinity, reflecting differences in binding enthalpy and entropy. The results reveal surprising differences in conformational and solvation entropy among the three complexes. NMR backbone order parameters show that the O-bound protein has reduced conformational entropy compared to the M and P complexes. By contrast, the bound ligand is more flexible in the O complex, as determined by ¹⁹F NMR relaxation, ensemble-refined X-ray diffraction data, and MD simulations. Furthermore, GIST calculations indicate that the O-bound complex has less unfavorable solvation entropy compared to the other two complexes. Thus, the results indicate compensatory effects from ligand conformational entropy and water entropy, on the one hand, and protein conformational entropy, on the other hand. Taken together, these different contributions amount to entropy-entropy compensation among the system components involved in ligand binding to a target protein.

Critical evaluation of a crystal structure of nitrogenase with bound N2 ligands

Recently, a 1.83 Å crystallographic structure of nitrogenase was suggested to show N₂-derived ligands at three sites in the catalytic FeMo cluster, replacing the three [Formula: see text] bridging sulfide ligands (two in one subunit and the third in the other subunit) (Kang et al. in Science 368: 1381-1385, 2020). Naturally, such a structure is sensational, having strong bearings on the reaction mechanism of the enzyme. Therefore, it is highly important to ensure that the interpretation of the structure is correct. Here, we use standard crystallographic refinement and quantum refinement to evaluate the structure. We show that the original crystallographic raw data are strongly anisotropic, with a much lower resolution in certain directions than others. This, together with the questionable use of anisotropic B factors, give atoms an elongated shape, which may look like diatomic atoms. In terms of standard electron-density maps and real-space Z scores, a resting-state structure with no dissociated sulfide ligands fits the raw data better than the interpretation suggested by the crystallographers. The anomalous electron density at 7100 eV is weaker for the putative N₂ ligands, but not lower than for several of the [Formula: see text] bridging sulfide ions and not lower than what can be expected from a statistical analysis of the densities. Therefore, we find no convincing evidence for any N₂ binding to the FeMo cluster. Instead, a standard resting state without any dissociated ligands seems to be the most likely interpretation of the structure. Likewise, we find no support that the homocitrate ligand should show monodentate binding.

Quantum-refinement studies of the bidentate ligand of V‑nitrogenase and the protonation state of CO-inhibited Mo‑nitrogenase

Nitrogenase is the only enzyme that can cleave the triple bond in N₂, making nitrogen available to plants (although the enzyme itself is strictly microbial). It has been studied extensively with both experimental and computational methods, but many details of the reaction mechanism are still unclear. X-ray crystallography is the main source of structural information for biomacromolecules, but it has problems to discern hydrogen atoms or to distinguish between elements with the same number of electrons. These problems can sometimes be alleviated by introducing quantum chemical calculations in the refinement, providing information about the ideal structure (in the same way as the empirical restraints used in standard crystallographic refinement) and comparing different interpretations of the structure with normal crystallographic and quantum mechanical quality measures. We have performed such quantum-refinement calculations to address two important issues for nitrogenase. First, we show that the bidentate ligand of the active-site FeV cluster in V‑nitrogenase is carbonate, rather than bicarbonate or nitrate. Second, we study the CO-inhibited structure of Mo‑nitrogenase. CO binds to a reduced and protonated state of the enzyme by replacing one of the sulfide ions (S2B) in the active-site FeMo cluster. We examined if it is possible to deduce from the crystal structure the location of the protons. Our results indicates that the crystal structure is best modelled as fully deprotonated.

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.

Quantum refinement with multiple conformations: application to the P-cluster in nitrogenase

X-ray crystallography is the main source of atomistic information on the structure of proteins. Normal crystal structures are obtained as a compromise between the X-ray scattering data and a set of empirical restraints that ensure chemically reasonable bond lengths and angles. However, such restraints are not always available or accurate for nonstandard parts of the structure, for example substrates, inhibitors and metal sites. The method of quantum refinement, in which these empirical restraints are replaced by quantum-mechanical (QM) calculations, has previously been suggested for small but interesting parts of the protein. Here, this approach is extended to allow for multiple conformations in the QM region by performing separate QM calculations for each conformation. This approach is shown to work properly and leads to improved structures in terms of electron-density maps and real-space difference density Z-scores. It is also shown that the quality of the structures can be gauged using QM strain energies. The approach, called ComQumX-2QM, is applied to the P-cluster in two different crystal structures of the enzyme nitrogenase, i.e. an Fe8S7Cys6 cluster, used for electron transfer. One structure is at a very high resolution (1.0 Å) and shows a mixture of two different oxidation states, the fully reduced PN state (Fe82+, 20%) and the doubly oxidized P2+ state (80%). In the original crystal structure the coordinates differed for only two iron ions, but here it is shown that the two states also show differences in other atoms of up to 0.7 Å. The second structure is at a more modest resolution, 2.1 Å, and was originally suggested to show only the one-electron oxidized state, P1+. Here, it is shown that it is rather a 50/50% mixture of the P1+ and P2+ states and that many of the Fe-Fe and Fe-S distances in the original structure were quite inaccurate (by up to 0.8 Å). This shows that the new ComQumX-2QM approach can be used to sort out what is actually seen in crystal structures with dual conformations and to give locally improved coordinates.

Automated orientation of water molecules in neutron crystallographic structures of proteins

The structure and function of proteins are strongly affected by the surrounding solvent water, for example through hydrogen bonds and the hydrophobic effect. These interactions depend not only on the position, but also on the orientation, of the water molecules around the protein. Therefore, it is often vital to know the detailed orientations of the surrounding ordered water molecules. Such information can be obtained by neutron crystallography. However, it is tedious and time-consuming to determine the correct orientation of every water molecule in a structure (there are typically several hundred of them), which is presently performed by manual evaluation. Here, a method has been developed that reliably automates the orientation of a water molecules in a simple and relatively fast way. Firstly, a quantitative quality measure, the real-space correlation coefficient, was selected, together with a threshold that allows the identification of water molecules that are oriented. Secondly, the refinement procedure was optimized by varying the refinement method and parameters, thus finding settings that yielded the best results in terms of time and performance. It turned out to be favourable to employ only the neutron data and a fixed protein structure when reorienting the water molecules. Thirdly, a method has been developed that identifies and reorients inadequately oriented water molecules systematically and automatically. The method has been tested on three proteins, galectin-3C, rubredoxin and inorganic pyrophosphatase, and it is shown that it yields improved orientations of the water molecules for all three proteins in a shorter time than manual model building. It also led to an increased number of hydrogen bonds involving water molecules for all proteins.

Is density functional theory accurate for lytic polysaccharide monooxygenase enzymes?

The lytic polysaccharide monooxygenase (LPMO) enzymes boost polysaccharide depolymerization through oxidative chemistry, which has fueled the hope for more energy-efficient production of biofuel. We have recently proposed a mechanism for the oxidation of the polysaccharide substrate (E. D. Hedegård and U. Ryde, Chem. Sci., 2018, 9, 3866-3880). In this mechanism, intermediates with superoxide, oxyl, as well as hydroxyl (i.e. [CuO2]+, [CuO]+ and [CuOH]2+) cores were involved. These complexes can have both singlet and triplet spin states, and both spin-states may be important for how LPMOs function during catalytic turnover. Previous calculations on LPMOs have exclusively been based on density functional theory (DFT). However, different DFT functionals are known to display large differences for spin-state splittings in transition-metal complexes, and this has also been an issue for LPMOs. In this paper, we study the accuracy of DFT for spin-state splittings in superoxide, oxyl, and hydroxyl intermediates involved in LPMO turnover. As reference we employ multiconfigurational perturbation theory (CASPT2).

Does the crystal structure of vanadium nitrogenase contain a reaction intermediate? Evidence from quantum refinement

Recently, a crystal structure of V-nitrogenase was presented, showing that one of the µ2 sulphide ions in the active site (S2B) is replaced by a lighter atom, suggested to be NH or NH2, i.e. representing a reaction intermediate. Moreover, a sulphur atom is found 7 Å from the S2B site, suggested to represent a storage site for this ion when it is displaced. We have re-evaluated this structure with quantum refinement, i.e. standard crystallographic refinement in which the empirical restraints (employed to ensure that the final structure makes chemical sense) are replaced by more accurate quantum-mechanical calculations. This allows us to test various interpretations of the structure, employing quantum-mechanical calculations to predict the ideal structure and to use crystallographic measures like the real-space Z-score and electron-density difference maps to decide which structure fits the crystallographic raw data best. We show that the structure contains an OH--bound state, rather than an N2-derived reaction intermediate. Moreover, the structure shows dual conformations in the active site with ~ 14% undissociated S2B ligand, but the storage site seems to be fully occupied, weakening the suggestion that it represents a storage site for the dissociated ligand.

N2H2 binding to the nitrogenase FeMo cluster studied by QM/MM methods

We have made a systematic combined quantum mechanical and molecular mechanical (QM/MM) investigation of possible structures of the N₂ bound state of nitrogenase. We assume that N₂ is immediately protonated to a N₂H₂ state, thereby avoiding the problem of determining the position of the protons in the cluster. We have systematically studied both end-on and side-on structures, as well as both HNNH and NNH₂ states. Our results indicate that the binding of N₂H₂ is determined more by interactions and steric clashes with the surrounding protein than by the intrinsic preferences of the ligand and the cluster. The best binding mode with both the TPSS and B3LYP density-functional theory methods has trans-HNNH terminally bound to Fe2. It is stabilised by stacking of the substrate with His-195 and Ser-278. However, several other structures come rather close in energy (within 3-35 kJ/mol) at least in some calculations: The corresponding cis-HNNH structure terminally bound to Fe2 is second best with B3LYP. A structure with HNNH₂ terminally bound to Fe6 is second most stable with TPSS (where the third proton is transferred to the substrate from the homocitrate ligand). Structures with trans-HNNH, bound to Fe4 or Fe6, or cis-HNNH bound to Fe6 are also rather stable. Finally, with the TPSS functional, a structure with cis-HNNH side-on binding to the Fe3-Fe4-Fe5-Fe7 face of the cluster is also rather low in energy, but all side-on structures are strongly disfavoured by the B3LYP method.

fragHAR: towards ab initio quantum-crystallographic X-ray structure refinement for polypeptides and proteins

The first ab initio aspherical structure refinement against experimental X-ray structure factors for polypeptides and proteins using a fragmentation approach to break up the protein into residues and solvent, thereby speeding up quantum-crystallographic Hirshfeld atom refinement (HAR) calculations, is described. It it found that the geometric and atomic displacement parameters from the new fragHAR method are essentially unchanged from a HAR on the complete unfragmented system when tested on dipeptides, tripeptides and hexapeptides. The largest changes are for the parameters describing H atoms involved in hydrogen-bond interactions, but it is shown that these discrepancies can be removed by including the interacting fragments as a single larger fragment in the fragmentation scheme. Significant speed-ups are observed for the larger systems. Using this approach, it is possible to perform a highly parallelized HAR in reasonable times for large systems. The method has been implemented in the TONTO software.

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.

What Is the Structure of the E4 Intermediate in Nitrogenase?

Nitrogenase is the only enzyme that can cleave the strong triple bond in N₂. The active site contains a complicated MoFe₇S₉C cluster. It is believed that it needs to accept four protons and electrons, forming the E₄ state, before it can bind N₂. However, there is no consensus on the atomic structure of the E₄ state. Experimental studies indicate that it should contain two hydride ions bridging two pairs of Fe ions, and it has been suggested that both hydride ions as well as the two protons bind on the same face of the cluster. On the other hand, density functional theory (DFT) studies have indicated that it is energetically more favorable with either three hydride ions or with a triply protonated carbide ion, depending on the DFT functional. We have performed a systematic combined quantum mechanical and molecular mechanical (QM/MM) study of possible E₄ states with two bridging hydride ions. Our calculations suggest that the most favorable structure has hydride ions bridging the Fe2/6 and Fe3/7 ion pairs. In fact, such structures are 14 kJ/mol more stable than structures with three hydride ions, showing that pure DFT functionals give energetically most favorable structures in agreement with experiments. An important reason for this finding is that we have identified a new type of broken-symmetry state that involves only two Fe ions with minority spin, in contrast to the previously studied states with three Fe ions with minority spin. The energetically best structures have the two hydride ions on different faces of the FeMo cluster, whereas better agreement with ENDOR data is obtained if they are on the same face; such structures are only 6–22 kJ/mol less stable.

Water structure in solution and crystal molecular dynamics simulations compared to protein crystal structures

The function of proteins is influenced not only by the atomic structure but also by the detailed structure of the solvent surrounding it. Computational studies of protein structure also critically depend on the water structure around the protein. Herein we compare the water structure obtained from molecular dynamics (MD) simulations of galectin-3 in complex with two ligands to crystallographic water molecules observed in the corresponding crystal structures. We computed MD trajectories both in a water box, which mimics a protein in solution, and in a crystallographic unit cell, which mimics a protein in a crystal. The calculations were compared to crystal structures obtained at both cryogenic and room temperature. Two types of analyses of the MD simulations were performed. First, the positions of the crystallographic water molecules were compared to peaks in the MD density after alignment of the protein in each snapshot. The results of this analysis indicate that all simulations reproduce the crystallographic water structure rather poorly. However, if we define the crystallographic water sites based on their distances to nearby protein atoms and follow these sites throughout the simulations, the MD simulations reproduce the crystallographic water sites much better. This shows that the failure of MD simulations to reproduce the water structure around proteins in crystal structures observed both in this and previous studies is caused by the problem of identifying water sites for a flexible and dynamic protein (traditionally done by overlaying the structures). Our local clustering approach solves the problem and shows that the MD simulations reasonably reproduce the water structure observed in crystals. Furthermore, analysis of the crystal MD simulations indicates a few water molecules that are close to unmodeled electron density peaks in the crystal structures, suggesting that crystal MD could be used as a complementary tool for identifying and modelling water in protein crystallography.

Molecular dynamics simulations can reproduce the water structure around proteins in crystal structure only if a local clustering is performed.

Refinement of protein structures using a combination of quantum-mechanical calculations with neutron and X-ray crystallographic data. Corrigendum

Corrections are published for the article by Caldararu et al. [(2019), Acta Cryst. D75, 368–380].

Why does sulfite reductase employ siroheme?

Sulfite reductase (SiR) contains in the active site a unique assembly of siroheme and a [4Fe4S] cluster, linked by a cysteine residue. Siroheme is a doubly reduced variant of heme that is not used for a catalytic function in any other enzyme. We have used non-equilibrium Green's function methods coupled with density functional theory computations to explain why SiR employs siroheme rather than heme. The results show that direct, through vacuum, charge-transfer routes are inhibited when heme is replaced by siroheme. This ensures more efficient channelling of the electrons to the catalytic iron during the six-electron reduction of sulfite to sulfide, limiting potential side-reactions that could occur if the incoming electrons were delocalized onto the macrocyclic ring.

Structure and Energetics of Ligand-Fluorine Interactions with Galectin-3 Backbone and Side-Chain Amides: Insight into Solvation Effects and Multipolar Interactions

Multipolar fluorine–amide interactions with backbone and side‐chain amides have been described as important for protein–ligand interactions and have been used to improve the potency of synthetic inhibitors. In this study, fluorine interactions within a well‐defined binding pocket on galectin‐3 were investigated systematically using phenyltriazolyl‐thiogalactosides fluorinated singly or multiply at various positions on the phenyl ring. X‐ray structures of the C‐terminal domain of galectin‐3 in complex with eight of these ligands revealed potential orthogonal fluorine–amide interactions with backbone amides and one with a side‐chain amide. The two interactions involving main‐chain amides seem to have a strong influence on affinity as determined by fluorescence anisotropy. In contrast, the interaction with the side‐chain amide did not influence affinity. Quantum mechanics calculations were used to analyze the relative contributions of these interactions to the binding energies. No clear correlation could be found between the relative energies of the fluorine–main‐chain amide interactions and the overall binding energy. Instead, dispersion and desolvation effects play a larger role. The results confirm that the contribution of fluorine–amide interactions to protein–ligand interactions cannot simply be predicted, on geometrical considerations alone, but require careful consideration of the energetic components.