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 19F 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.

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.

Three Simple Properties Explain Protein Stability Change upon Mutation.

Accurate prediction of protein stability upon mutation enables rational engineering of new proteins and insights into protein evolution and monogenetic diseases caused by single-point amino acid substitutions. Many tools have been developed to this aim, ranging from energy-based models to machine-learning methods that use large amounts of experimental data. However, as the methods become more complex, the interpretation of the chemistry underlying the protein stability effects becomes obscure. It is thus of interest to identify the simplest prediction model that retains complete amino acid specific interpretation; for a given number of input descriptors, we expect such a model to be almost universal. In this study, we identify such a limiting model, SimBa, a simple multilinear regression model trained on a substitution-type-balanced experimental data set. The model accounts only for the solvent accessibility of the site, volume difference, and polarity difference caused by mutation. Our results show that this very simple and directly applicable model performs comparably to other much more complex, widely used protein stability prediction methods. This suggests that a hard limit of ∼1 kcal/mol numerical accuracy and an R ∼ 0.5 trend accuracy exists and that new features, such as account of unfolded states, water colocalization, and amino acid correlations, are required to improve accuracy to, e.g., 1/2 kcal/mol.

A base measure of precision for protein stability predictors: structural sensitivity.

BACKGROUND: Prediction of the change in fold stability (ΔΔG) of a protein upon mutation is of major importance to protein engineering and screening of disease-causing variants. Many prediction methods can use 3D structural information to predict ΔΔG. While the performance of these methods has been extensively studied, a new problem has arisen due to the abundance of crystal structures: How precise are these methods in terms of structure input used, which structure should be used, and how much does it matter? Thus, there is a need to quantify the structural sensitivity of protein stability prediction methods. RESULTS: We computed the structural sensitivity of six widely-used prediction methods by use of saturated computational mutagenesis on a diverse set of 87 structures of 25 proteins. Our results show that structural sensitivity varies massively and surprisingly falls into two very distinct groups, with methods that take detailed account of the local environment showing a sensitivity of ~ 0.6 to 0.8 kcal/mol, whereas machine-learning methods display much lower sensitivity (~ 0.1 kcal/mol). We also observe that the precision correlates with the accuracy for mutation-type-balanced data sets but not generally reported accuracy of the methods, indicating the importance of mutation-type balance in both contexts. CONCLUSIONS: The structural sensitivity of stability prediction methods varies greatly and is caused mainly by the models and less by the actual protein structural differences. As a new recommended standard, we therefore suggest that ΔΔG values are evaluated on three protein structures when available and the associated standard deviation reported, to emphasize not just the accuracy but also the precision of the method in a specific study. Our observation that machine-learning methods deemphasize structure may indicate that folded wild-type structures alone, without the folded mutant and unfolded structures, only add modest value for assessing protein stability effects, and that side-chain-sensitive methods overstate the significance of the folded wild-type structure.

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.

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.

Systematic Investigation of the Data Set Dependency of Protein Stability Predictors.

Prediction of protein stability changes caused by mutation is of major importance to protein engineering and for understanding protein misfolding diseases and protein evolution. The major limitation to these applications is the fact that different prediction methods vary substantially in terms of performance for specific proteins; i.e., performance is not transferable from one type of mutation or protein to another. In this study, we investigated the performance and transferability of eight widely used methods. We first constructed a new data set composed of 2647 mutations using strict selection criteria for the experimental data and then defined a variety of subdata sets that are unbiased with respect to various aspects such as mutation type, stabilization extent, structure type, and solvent exposure. Benchmarking the methods against these subdata sets enabled us to systematically investigate how data set biases affect predictor performance. In particular, we use a reduced amino acid alphabet to quantify the bias toward mutation type, which we identify as the major bias in current approaches. Our results show that all prediction methods exhibit large biases, stemming not from failures of the models applied but mostly from the selection biases of experimental data used for training or parametrization. Our identification of these biases and the construction of new mutation-type-balanced data should lead to the development of more balanced and transferable prediction methods in the future.

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].

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.

Are crystallographic B-factors suitable for calculating protein conformational entropy?

Conformational entropies are of great interest when studying the binding of small ligands to proteins or the interaction of proteins. Unfortunately, there are no experimental methods available to measure conformational entropies of all groups in a protein. Instead, they are normally estimated from molecular dynamics (MD) simulations, although such methods show problems with convergence and correlation of motions, and depend on the accuracy of the underlying potential-energy function. Crystallographic atomic displacement parameters (also known as B-factors) are available in all crystal structures and contain information about the atomic fluctuations, which can be converted to entropies. We have studied whether B-factors can be employed to extract conformational entropies for proteins by comparing such entropies to those measured by NMR relaxation experiments or obtained from MD simulations in solution or in the crystal. Unfortunately, our results show that B-factor entropies are unreliable, because they include the movement and rotation of the entire protein, they exclude correlation of the movements and they include contributions other than the fluctuations, e.g. static disorder, as well as errors in the model and the scattering factors. We have tried to reduce the first problem by employing translation-libration-screw refinement, the second by employing a description of the correlated movement from MD simulations, and the third by studying only the change in entropy when a pair of ligands binds to the same protein, thoroughly re-refining the structures in exactly the same way and using the same set of alternative conformations. However, the experimental B-factors seem to be incompatible with fluctuations from MD simulations and the precision is too poor to give any reliable entropies.

Geometry and Electronic Structure of the P-Cluster in Nitrogenase Studied by Combined Quantum Mechanical and Molecular Mechanical Calculations and Quantum Refinement.

We have studied the geometry and electronic structure of the P-cluster in nitrogenase in four oxidation states: PN, P1+, P2+, and P3+. We have employed combined quantum mechanical and molecular mechanical (QM/MM) calculations, using two different density-functional theory methods, TPSS and B3LYP. The calculations confirm that the side chain of Ser-188 is most likely deprotonated in the partly oxidized P1+ state, thereby forming a bond to Fe6. Likewise, the backbone amide group of Cys-88 is deprotonated in the doubly oxidized P2+ state, forming a bond to Fe5. The calculations also confirm the two conformations of the P-cluster in the atomic-resolution crystal structure of the enzyme, representing the PN and P2+ states, but show that the finer differences between the two structures are not fully reflected in the crystal structure, because the coordinates of only two atoms differ between the two conformations. However, the recent crystal structure of the P1+ state seems to be of lower quality with many dubious Fe-Fe and Fe-S distances. Quantum refinement of this structure indicates that it is a mixture of the P1+ and P2+ states but confirms that the side chain of Ser-188 is most likely deprotonated in both states. TPSS gives structures that are appreciably closer to the crystal structures than does B3LYP. In addition, we have studied all 16-48 possible broken-symmetry states of the four oxidation states of the P-cluster with DFT in the one or two observed spin states. For the reduced PN state, we can settle the most likely state from the calculated energies and geometries. However, for the more oxidized states there are large differences in the predictions obtained with the two DFT methods.

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

Neutron crystallography is a powerful method to determine the positions of H atoms in macromolecular structures. However, it is sometimes hard to judge what would constitute a chemically reasonable model, and the geometry of H atoms depends more on the surroundings (for example the formation of hydrogen bonds) than heavy atoms, so that the empirical geometry information for the H atoms used to supplement the experimental data is often less accurate. These problems may be reduced by using quantum-mechanical calculations. A method has therefore been developed to combine quantum-mechanical calculations with joint crystallographic refinement against X-ray and neutron data. A first validation of this method is provided by re-refining the structure of the galectin-3 carbohydrate-recognition domain in complex with lactose. The geometry is improved, in particular for water molecules, for which the method leads to better-resolved hydrogen-bonding interactions. The method has also been applied to the active copper site of lytic polysaccharide monooxygenase and shows that the protonation state of the amino-terminal histidine residue can be determined.

Mechanism of hydrogen peroxide formation by lytic polysaccharide monooxygenase.

Lytic polysaccharide monooxygenases (LPMOs) are copper-containing metalloenzymes that can cleave the glycosidic link in polysaccharides. This could become crucial for production of energy-efficient biofuels from recalcitrant polysaccharides. Although LPMOs are considered oxygenases, recent investigations have shown that H2O2 can also act as a co-substrate for LPMOs. Intriguingly, LPMOs generate H2O2 in the absence of a polysaccharide substrate. Here, we elucidate a new mechanism for H2O2 generation starting from an AA10-LPMO crystal structure with an oxygen species bound, using QM/MM calculations. The reduction level and protonation state of this oxygen-bound intermediate has been unclear. However, this information is crucial to the mechanism. We therefore investigate the oxygen-bound intermediate with quantum refinement (crystallographic refinement enhanced with QM calculations), against both X-ray and neutron data. Quantum refinement calculations suggest a Cu(ii)-O-2 system in the active site of the AA10-LPMO and a neutral protonated -NH2 state for the terminal nitrogen atom, the latter in contrast to the original interpretation. Our QM/MM calculations show that H2O2 generation is possible only from a Cu(i) center and that the most favourable reaction pathway is to involve a nearby glutamate residue, adding two electrons and two protons to the Cu(ii)-O-2 system, followed by dissociation of H2O2.

Interplay between Conformational Entropy and Solvation Entropy in Protein-Ligand Binding.

Understanding the driving forces underlying molecular recognition is of fundamental importance in chemistry and biology. The challenge is to unravel the binding thermodynamics into separate contributions and to interpret these in molecular terms. Entropic contributions to the free energy of binding are particularly difficult to assess in this regard. Here we pinpoint the molecular determinants underlying differences in ligand affinity to the carbohydrate recognition domain of galectin-3, using a combination of isothermal titration calorimetry, X-ray crystallography, NMR relaxation, and molecular dynamics simulations followed by conformational entropy and grid inhomogeneous solvation theory (GIST) analyses. Using a pair of diastereomeric ligands that have essentially identical chemical potential in the unbound state, we reduced the problem of dissecting the thermodynamics to a comparison of the two protein-ligand complexes. While the free energies of binding are nearly equal for the R and S diastereomers, greater differences are observed for the enthalpy and entropy, which consequently exhibit compensatory behavior, ΔΔ H°(R - S) = -5 ± 1 kJ/mol and - TΔΔ S°(R - S) = 3 ± 1 kJ/mol. NMR relaxation experiments and molecular dynamics simulations indicate that the protein in complex with the S-stereoisomer has greater conformational entropy than in the R-complex. GIST calculations reveal additional, but smaller, contributions from solvation entropy, again in favor of the S-complex. Thus, conformational entropy apparently dominates over solvation entropy in dictating the difference in the overall entropy of binding. This case highlights an interplay between conformational entropy and solvation entropy, pointing to both opportunities and challenges in drug design.

Protonation and Reduction of the FeMo Cluster in Nitrogenase Studied by Quantum Mechanics/Molecular Mechanics (QM/MM) Calculations.

We have performed a systematic computational study of the relative energies of possible protonation states of the FeMo cluster in nitrogenase in the E0-E4 states, i.e., the resting state and states with 1-4 electrons and protons added but before N2 binds. We use the combined quantum mechanics and molecular mechanics (QM/MM) approach, including the complete solvated heterotetrameric enzyme in the calculations. The QM system consisted of 112 atoms, i.e., the full FeMo cluster, as well all groups forming hydrogen bonds to it within 3.5 Å. It was treated with either the TPSS-D3 or B3LYP-D3 methods with the def2-SV(P) or def2-TZVPD basis sets. For each redox state, we calculated relative energies of at least 50 different possible positions for the proton, added to the most stable protonation state of the level with one electron less. We show quite conclusively that the resting E0 state is not protonated using quantum refinement and by comparing geometries to the crystal structure. The E1 state is protonated on S2B, in agreement with most previous computational studies. However, for the E2-E4 states, the two QM methods give diverging results, with relative energies that differ by over 300 kJ/mol for the most stable E4 states. TPSS favors hydride ions binding to the Fe ions. The first bridges Fe2 and Fe6, whereas the next two bind terminally to either Fe4, Fe5, or Fe6 with nearly equal energies. On the other hand, B3LYP disfavors hydride ions and instead suggests that 1-3 protons bind to the central carbide ion.

Binding free energies in the SAMPL6 octa-acid host-guest challenge calculated with MM and QM methods.

We have estimated free energies for the binding of eight carboxylate ligands to two variants of the octa-acid deep-cavity host in the SAMPL6 blind-test challenge (with or without endo methyl groups on the four upper-rim benzoate groups, OAM and OAH, respectively). We employed free-energy perturbation (FEP) for relative binding energies at the molecular mechanics (MM) and the combined quantum mechanical (QM) and MM (QM/MM) levels, the latter obtained with the reference-potential approach with QM/MM sampling for the MM → QM/MM FEP. The semiempirical QM method PM6-DH+ was employed for the ligand in the latter calculations. Moreover, binding free energies were also estimated from QM/MM optimised structures, combined with COSMO-RS estimates of the solvation energy and thermostatistical corrections from MM frequencies. They were performed at the PM6-DH+ level of theory with the full host and guest molecule in the QM system (and also four water molecules in the geometry optimisations) for 10-20 snapshots from molecular dynamics simulations of the complex. Finally, the structure with the lowest free energy was recalculated using the dispersion-corrected density-functional theory method TPSS-D3, for both the structure and the energy. The two FEP approaches gave similar results (PM6-DH+/MM slightly better for OAM), which were among the five submissions with the best performance in the challenge and gave the best results without any fit to data from the SAMPL5 challenge, with mean absolute deviations (MAD) of 2.4-5.2 kJ/mol and a correlation coefficient (R2) of 0.77-0.93. This is the first time QM/MM approaches give binding free energies that are competitive to those obtained with MM for the octa-acid host. The QM/MM-optimised structures gave somewhat worse performance (MAD = 3-8 kJ/mol and R2 = 0.1-0.9), but the results were improved compared to previous studies of this system with similar methods.

QM/MM study of the reaction mechanism of sulfite oxidase.

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

Quantum Refinement Does Not Support Dinuclear Copper Sites in Crystal Structures of Particulate Methane Monooxygenase.

Particulate methane monooxygenase (pMMO) is one of the few enzymes that can activate methane. The metal content of this enzyme has been highly controversial, with suggestions of a dinuclear Fe site or mono-, di-, or trinuclear Cu sites. Crystal structures have shown a mono- or dinuclear Cu site, but the resolution was low and the geometry of the dinuclear site unusual. We have employed quantum refinement (crystallographic refinement enhanced with quantum-mechanical calculations) to improve the structure of the active site. We compared a number of different mono- and dinuclear geometries, in some cases enhanced with more protein ligands or one or two water molecules, to determine which structure fits two sets of crystallographic raw data best. In all cases, the best results were obtained with mononuclear Cu sites, occasionally with an extra water molecule. Thus, we conclude that there is no crystallographic support for a dinuclear Cu site in pMMO.

Protonation States of Homocitrate and Nearby Residues in Nitrogenase Studied by Computational Methods and Quantum Refinement.

Nitrogenase is the only enzyme that can break the triple bond in N2 to form two molecules of ammonia. The enzyme has been thoroughly studied with both experimental and computational methods, but there is still no consensus regarding the atomic details of the reaction mechanism. In the most common form, the active site is a MoFe7S9C(homocitrate) cluster. The homocitrate ligand contains one alcohol and three carboxylate groups. In water solution, the triply deprotonated form dominates, but because the alcohol (and one of the carboxylate groups) coordinate to the Mo ion, this may change in the enzyme. We have performed a series of computational calculations with molecular dynamics (MD), quantum mechanical (QM) cluster, combined QM and molecular mechanics (QM/MM), QM/MM with Poisson-Boltzmann and surface area solvation, QM/MM thermodynamic cycle perturbations, and quantum refinement methods to settle the most probable protonation state of the homocitrate ligand in nitrogenase. The results quite conclusively point out a triply deprotonated form (net charge -3) with a proton shared between the alcohol and one of the carboxylate groups as the most stable at pH 7. Moreover, we have studied eight ionizable protein residues close to the active site with MD simulations and determined the most likely protonation states.