Protonation States of Homocitrate and Nearby Residues in

Subscriber access provided by UNIVERSITY OF ADELAIDE LIBRARIES

Article

Protonation States of Homocitrate and Nearby Residues in Nitrogenase Studied by Computational Methods and Quantum Refinement Lili Cao, Octav Caldararu, and Ulf Ryde J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b02714 • Publication Date (Web): 07 Aug 2017 Downloaded from http://pubs.acs.org on August 18, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

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

Lili Cao, Octav Caldararu and Ulf Ryde *

Department of Theoretical Chemistry, Lund University, Chemical Centre, P. O. Box 124, SE-221 00 Lund, Sweden

Correspondence to Ulf Ryde, E-mail: [email protected], Tel: +46 – 46 2224502, Fax: +46 – 46 2228648

2017-08-07

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract 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 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 most stable at pH 7. Moreover, we have studied eight ionisable protein residues close to the active site with MD simulations and determined the most likely protonation states.


Page 2 of 46

Page 3 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction Nitrogen is an essential component of many biomolecules, including all amino acids and the nucleotide components of DNA and RNA. Molecular nitrogen gas (N2) is the major component of the atmosphere (78%), but still nitrogen is often the limiting element for plant growth and a main component of synthetic fertilisers.1 The reason for this is that the triple bond in N2 is so strong that it cannot be cleaved by the great majority of the organisms. Instead, it needs to be taken up in the form of ammonia or nitrate. Industrially, ammonia is formed from N2 and H2 by the Haber–Bosch process, which requires high temperature, high pressure and an iron catalyst.2 The process was invented in 1909–1910 and it is one of the most important industrial processes today, consuming ~1% of the world’s energy supplies and providing half of the total biologically available nitrogen on earth.3 Thus, it is a major factor in the recent agricultural revolution and the human population explosion.1 The enzyme nitrogenase (EC 1.18/19.6.1) is present in a few groups of bacteria and archaea.1,4,5 Several of these live in symbiosis with higher plants, such as legumes, rice, sugarcane and alder. It can perform the conversion of N2 to ammonia at ambient pressure and temperature. However, the process is very energy demanding, consuming 16 molecules of ATP for each molecule of N2 fixed according to the reaction: N2 + 8 e– + 8 H+ + 16 ATP → 2 NH3 + H2 + 16 ADP + 16 Pi

(1)

from which it can be seen that H2 is a byproduct of the reaction. Several crystal structures of nitrogenases have been determined6–10 showing that it is a α2β2 heterotetramer. It is now agreed that the active site consists of the unusual MoFe7CS9(homocitrate) MoFe cluster,8,11 which is connected to the enzyme by one cysteine (Cys) and one histidine (His) ligand, as can be seen in Figure 1a. In some enzymes, the Mo ion is replaced by vanadium or iron.12 The protein also contains an [8Fe7S] electron-transfer cluster, called the P cluster, which is shown in Figure 1b. The electrons are provided by another protein, called the Fe protein, which also binds two ATP molecules, triggering the docking to the MoFe protein and electron transfer to the MoFe cluster, via the P cluster. Hydrolysis of the ATP molecules triggers the dissociation of the Fe protein, opening up for additional electron transfers. The nitrogenases have been extensively studied by a great wealth of biochemical, spectroscopic and kinetic methods.1,4,5,11,13 The consensus is that the enzyme follows the Lowe–Thorneley scheme, involving nine intermediates E0–E8 differing in the number of electrons and protons delivered to the MoFe cluster.14 Several of these intermediates have been trapped and spectroscopically characterised.1 Of particular interest is E4, which is believed to be the species that binds N2. The experimental studies have also been supplemented by numerous quantum-mechanical (QM) studies, which may give an atomic picture of the reaction mechanism.1,13,15–27 Unfortunately, they have involved strongly diverging mechanisms. For example, Dance has suggested that N2 binds side-on to one Fe ion and is protonated alternatively on the two N atoms, forming intermediates that bridge two Fe ions after the first proton has been added.28 On the other hand, Nørskov and coworkers have proposed a mechanism in which N2 binds with one N atom bridging two Fe ions, after the dissociation of one of the sulfide ions, and the non-coordinating N atom is fully protonated and dissociates before the other N atom is protonated.24 Furthermore, Siegbahn has propounded that Fe binds side-on, bridging between two Fe ions, after a triple protonation of the central carbide ion, which thereby moves to the periphery of the cluster.25 Moreover, McKee has suggested that N2 instead binds end-on to the central C atom and is protonated first on the distal N atom.29 Likewise, Rao et al. has proposed that the central C is singly protonated and covalently interacts with N2, which is 3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

protonated in a sequential manner.27 Thus, there is no consensus regarding the detailed mechanism of nitrogenase, besides that four electrons and protons need to be added to the cluster before N2 can bind. A great majority of the QM investigations have ignored the surrounding protein or modelled it as a featureless continuum solvent (a few studies have included a small set of amino acids close to the MoFe cluster13,19,25,27). This may give a too large flexibility of the cluster, allowing for extensive conformational changes of the cluster that are not possible inside the protein. Moreover, it ignores the influence of the detailed electrostatics on the protonation of the various sulfide ions in the cluster, which are crucial for the detailed mechanism. The only exception is the combined QM and molecular mechanical (MM; QM/MM) by Xie et al.20 They noted that the “protein environment is important for computational characterisation of the structure of the MoFe cofactor”, but reached the incorrect conclusion that the central atom was O, rather than C. Rao et al. used an intermediate approach in which all atoms within 8 Å of the MoFe cluster (780 atoms) were considered at the PM6 semiempirical level of theory in QM/QM calculations, but keeping most atoms fixed.27 There are several problems when setting up QM/MM calculations of the whole nitrogenase. First, the protein is large, with around 2000 amino acids. Second, the protonation state of all residues in the protein needs to be determined. This is particularly important for the homocitrate ligand of the Mo ion in the MoFe cofactor, which can be expected to strongly influence the reactivity of the active site. It involves one alcohol and three carboxylate groups, of which one carboxylates and the alcohol coordinate to Mo. In water solution, the alcohol is protonated and the three carboxylate groups are deprotonated, although the third pKa value is close to 7, at least for the closely related citric-acid molecule.30 However, this may change when the molecule coordinates to Mo (the pKa of water changes from 15.6 to 2.2 when bound to Fe3+ 31). In previous QM studies, the homocitrate ligand has been modelled in many different ways, e.g. by two OH– ions,24 by a singly32 or doubly17,28 deprotonated glycolate molecule, or by the whole homocitrate molecule with a charge of –2,25 –313 or –4.13,19,20,27 Likewise, there are several buried histidine and aspartate or glutamate groups close to the MoFe cluster that may have unusual protonation states. Third, nitrogenase contains three metal sites, the MoFe cluster, the P cluster and a Ca2+ site (Figure 1), which are problematic to model with MM methods. Forth, both the MoFe and P clusters have complicated electronic structures, involving antiferromagnetically coupled metal ions. These are normally treated by the broken-symmetry (BS) approach in QM density-functional theory (DFT) calculations.19 However, these complicated clusters may have many BS solutions. For an isolated MoFe cluster, there are ten BS states, owing to the approximate three-fold symmetry of the cluster,19 but in the protein, this symmetry disappears and there are instead 35 possible BS states. These may change when the cluster is reduced or protonated. In this paper, we examine these problems with a series of molecular dynamics (MD) simulations, QM-cluster, QM/MM, big-QM,33 QTCP (QM/MM thermodynamic cycle perturbation),34 as well as quantum-refinement calculations.35 In particular, we decide the most likely protonation state of the homocitrate ligand, as well as all the protein residues around the MoFe and P clusters.


Page 4 of 46

Page 5 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Methods The protein All calculations were based on the 1.0-Å crystal structure of nitrogenase from Azotobacter vinelandii (PDB code 3U7Q).8 The entire α2β2 heterotetramer was included in the calculations, because the various subunits are entangled without any natural way to separate them. No attempt was made to model the missing starting and ending residues of the subunits A and C, whereas the B and D residues were complete and were modelled with charged amino and carboxy terminals. Eight imidazole molecules from the buffer were deleted, as well as two Mg2+ ions on the surface of the protein. On the other hand, the two Ca2+ ions were kept, because they are deeply buried in the protein and stabilise the interface between the subunits. All crystal-water molecules were kept in the calculations, except 26 that overlapped with each other or with protein atoms: HOH-1157A, 1158A, 1729A, 1743A, 2215A, 2228A, 1264B, 1650B, 1661B, 1861B, 2138B, 2214B, 2231B, 2302B, 2325B, 2340B, 2383B, 1976, 2037, 2262, 612D, 1513D, 1591D, 1770D, 2175D and 2235D. When discussing residues, a letter “A, “B”, or “D” after the residue number indicates that it belongs to that subunit. If no letter is given, it belongs to subunit C. In general, we will not explicitly discuss subunits A and B (e.g. for the protonation states), because they are identical to the C and D subunits. The QM calculations were concentrated on the MoFe cluster in the C subunit, because there is a buried imidazole molecule from the solvent rather close to the MoFe cluster (~11 Å) in the A subunit. MoFe and P clusters not treated by QM were modelled by MM in the resting and fully reduced state, respectively. The protonation states of all the residues were determined from a detailed study of the hydrogen-bond pattern and the solvent accessibility. It was checked by the PROPKA36 and Maestro37 software. All Arg, Lys, Asp, and Glu residues were assumed to be charged, except Glu-153, 440, and 231D. Cys residues coordinating to Fe ions were assumed to be deprotonated. A thorough manual investigation of all the other His residues gave the following protonation assignment: His-271, 451, 297D, 359D and 519D were assumed to be protonated on the ND1 atom (we denote atoms by their PDB names), His-31, 196, 285, 383, 90D, 185D, 363D and 457D were presumed to be protonated on both the ND1 and NE2 atoms (and therefore positively charged), whereas the remaining 14 His residues were modelled with a proton on the NE2 atom. Residues His-31, 363D and 429D were flipped (i.e. the C and N atoms in the imidazole ring were exchanged). Protons were added by the Maestro software,37 optimising the hydrogen-bond network. Maestro also flipped the Asn and Gln sidechains of residues 49, 119, 199, 252, 271, 384, 432, 468, 469, 58D, 71D, 104D, 130D, 136D, 167D, 225D, 268D, 278D, 286D, 294D, 317D, 445D, 499D and 513D.

QM calculations All QM calculations were performed with the Turbomole software (versions 7.1 and 7.2).38 We employed two DFT methods, TPSS39 and B3LYP,40–42 and two different basis sets of increasing size, def2-SV(P)43 and def2-TZVPD.44 The calculations were sped up by expanding the Coulomb interactions in an auxiliary basis set, the resolution-of-identity (RI) approximation.45,46 Empirical dispersion corrections were included with the DFT-D3 approach47 and Becke–Johnson damping,48 as implemented in Turbomole. The MoFe cluster was modelled by MoFe7S9C(homocitrate)(CH3S)(imidazole) (53 atoms, cf. Figure 1a), where the two last groups are models of Cys-275 and His-442. Following previous QM investigations,13,25 we used the oxidation state-assignment MoIII FeII FeIII . Four different protonation states of homocitrate were tested: fully deprotonated (0H), with one 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 46

proton either on the alcohol (1Ha) or on the O2 carboxylate atom (1Hc) or with protons on both these groups (2H). These protonation states are shown in Figure 2. They give net charges for the QM cluster of –6 to –4. The cluster was modelled in the S = 3/2 spin state with three unpaired electrons, which is the experimental ground state.1,13 The electronic structure was obtained with the BS approach.19 In practice, this means that each of the seven Fe ions are in the high-spin state, with either a surplus of α (four of them) or β (three of them) spin. This can be selected in 35 ! different ways (! !), which were all tested. A starting wavefunction was obtained by first optimising the all-high-spin state with 35 unpaired electrons and then changing the occupation numbers to the quartet state. Alternatively, we also used the fragment approach suggested by Szilagyi and Winslow.49 The other spin states were then obtained by simply swapping the coordinates of the Fe ions.50 A detailed account on the influence of the protein, QM method and basis set on the various BS states will be published elsewhere.51 All calculations presented in this paper are based on the lowest-energy BS7 state,52 in the nomenclature of Noodleman,19 with a surplus of β spin on the Fe2, Fe3 and Fe5 atoms (atom names are shown in Figure 1a and spin densities are given in Table S1 in the supplementary material). For the 1Hc state, BS7 is 26 kJ/mol more stable than the second most stable BS6 state at the TPSSD3/def2-TZVPD level of theory.51 The results in Table S2 show that the same applies for the 2H state (41 and 33 kJ/mol more stable than the next BS state with the TPSS and B3LYP functionals, respectively). The P cluster was modelled by [Fe8S7(CH3S)6]4–, i.e. the all-reduced state. It was also treated by the BS approach with all eight Fe ions in the high-spin state, antiferromagnetically coupled to an open-shell singlet state. All possible 6*6 = 36 broken-symmetry states were considered, assuming that each of the two Fe4 sub-clusters couple internally to a singlet state. The lowest energy was found for the state with a surplus of β spin on the Fe1, Fe3, Fe5 and Fe8 atoms (cf. Figure 1b). In some calculations, the QM system was immersed into a continuum solvent, employing the conductor-like screening model (COSMO),53,54 implemented in Turbomole. The default optimised COSMO radii were employed and a water solvent radius of 1.3 Å,55 whereas a radius of 2.0 Å was used for the metals.56 All geometry optimisations and QTCP calculations were performed at the TPSS-D3/def2SV(P) level of theory. The basis-set and DFT-functional dependence was examined by performing single-point energy calculations with the TPSS-D3/def2-TZVPD and B3LYPD3/def2-TZVPD methods. The acidity constant (pKa value) of an AH/A– couple was calculated according to p =

∆A– ∆AH ∆H

(2)

where ∆G(H+) = −1131.0 kJ/mol, which includes the hydration free energy of a proton, the translational Gibbs free energy of a proton and the change in reference state from 1 atm to 1 M, all calculated at 300 K and 1 atm pressure.57 The QM energies of the protonated (AH) and deprotonated molecule (A–) were calculated in the COSMO continuum solvent and thermal corrections to the Gibbs free energy (including the zero-point energy) were obtained from the vibrational frequencies.58

MD simulations All MM and MD simulations were performed with the Amber software.59 For the protein, we used the Amber ff14SB force field60 and water molecules were described by the TIP3P 6 ACS Paragon Plus Environment

Page 7 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


model.61 For the metal sites, restrained electrostatic potential charges were employed,62 using electrostatic potentials calculated at the TPSS/def2-SV(P) level of theory39,43 and sampled with the Mertz–Kollman scheme,63 although at a higher-than-default point density (~2000/atom56). The calculations were performed on the [Ca(OOC(CH2)3NHCOCH3)(CH3COO)2(H2O)3]–, [Fe8S7(CH3S)6]4– and MoFe7S9C(homocitrate)(CH3S)(imidazole) cluster models, shown in Figure 1. The charges are listed in Tables S3–S5 and the van der Waals parameters in Table S6.64 The metal sites were treated by a non-bonded model with restraints between the metal and the ligands,65 as well as between the nearby metals. This is necessary to keep the structure intact for these complicated metal clusters. For the restraints, we used the metal–ligand distances observed in the crystal structure (average over the two subunits) and the force constants listed in Table S7. The latter were obtained from TPSS/def2-SV(P) frequency calculations of the models in Figure 1, calculated by the method of Seminario66,67 and averaged over interactions of the same type. For all metal–metal distances, the same force constant of 125.5 kJ/mol/Å2 was used with four different distances: 2.59 Å for the three shortest Fe–Fe distances in the MoFe cluster, 2.65 Å for all the other Fe–Fe interactions in the MoFe cluster, 2.66 for the three Mo–Fe interactions in the MoFe cluster and also for all Fe–Fe interactions in the P cluster, except the two longest, for which 2.94 Å was used. Two types of simulations were set up. For the MD simulations, the protein was solvated in a periodic truncated octahedral box of explicit water molecules, extending at least 10 Å from the solute using the leap program in the Amber suite.59 185 Cl– ions and around 207 Na+ ions (depending on the protonation state of homocitrate and the examined His, Asp and Glu residues) were added to neutralise the protein and give an ionic strength of 0.2 M.68 The ions were added by replacing random water molecules with the leap program in the Amber software. However, some of them then end up inside the protein. This was avoided by using local software, ensuring that all counter ions were in the solvent. The final system contained ~154 850 atoms. After the solvation, we performed 1000 cycles of minimisation, with the heavy atoms of the protein restrained towards the starting structure with a force constant of 418 kJ/mol/Å2. This was followed by 20 ps constant-volume and 20 ps constant-pressure equilibration with the same restraints, but a force constant of 209 kJ/mol/Å2. Finally, the system was equilibrated for 1 ns without any restraints, followed by a 10-ns production simulation, during which coordinates were sampled every 10 ps. The temperature was kept constant at 300 K using Langevin dynamics, with a collision frequency of 2 ps–1.69 The pressure was kept constant at 1 atm using Berendsen’s weakcoupling isotropic algorithm with a relaxation time of 1 ps.70 Long-range electrostatics were handled by particle-mesh Ewald summation71 with a fourth-order B spline interpolation and a tolerance of 10–5. The cut-off radius for Lennard‒Jones interactions was set to 8 Å. All bonds involving hydrogen atoms were constrained to their equilibrium values using the SHAKE algorithm,72 allowing for a time step of 2 fs during the simulations. For the QM/MM calculations, the protein was instead solvated in a sphere with a radius of 70 Å around the geometrical centre of the protein. Moreover, all solvent-exposed Asp, Glu, Lys, Arg, and His residues were neutralised (a list of the buried charged groups are given in Table S8). 160 Cl– and ~180 Na+ ions were added in the same way as described above. The final system contained ~133 900 atoms. The added protons and water molecules were optimised by a 240 ps simulated annealing calculation (up to 370 K), followed by a minimisation, keeping the other atoms fixed at the crystal-structure positions. No SHAKE restraints were employed and the time step was 0.5 fs. Water molecules were kept inside the spherical system by a half-harmonic potential with a force constant of 6.3 kJ/mol/Å2, but the radius was reduced to 65 Å to avoid voids in the system. 7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

In the analysis of the MD simulations, the cut-off for hydrogen bonds was a H–O/N distance of less than 2.5 Å and a H–O/N–X angle larger than 90º.

QM/MM-2QM calculations The QM/MM-2QM calculations were performed with the ComQum software.73,74 In this approach, the protein and solvent are split into three subsystems: the two QM systems (systems 1a and 1b) and the MM system (system 2). In the QM calculations, systems 1a and 1b were represented by a wavefunction, whereas all the other atoms were represented by an array of partial point charges, one for each atom, taken from MM libraries. Thereby, the polarisation of the QM system by the surroundings is included in a self-consistent manner (electrostatic embedding, EE). In the QM/MM-2QM calculations, the MM system was kept fixed at the original (crystallographic) coordinates (but it was relaxed in the subsequent QTCP calculations. When there is a bond between systems 2 and 1a or 1b (a junction), the hydrogen linkatom approach was employed: The QM system was capped with hydrogen atoms (hydrogen link atoms, HL), the positions of which are linearly related to the corresponding carbon atoms (carbon link atoms, CL) in the full system.73,75 All atoms were included in the point-charge model, except the CL atoms.76 The total QM/MM-2QM energy in ComQum was calculated as73,74,77 CL HL HL HL HL !QM/MM = !QM1a+ptch1b2 − !MM1a,q + !QM1b+ptch1a2 − !MM1b,q + !MM1ab2,q 1a =0 1a =0,q 1b =0 1b =0 (3) HL where !QM1a+ptch1b2 is the QM energy of the QM system 1a truncated by HL atoms and embedded in the set of point charges modelling systems 1b and 2 (but excluding the selfHL energy of the point charges). !QM1a+ptch1b2 is the corresponding QM energy of system 1b. HL 67 !223,435 and !MM1b,q1b =0 are the MM energies of the QM systems, still truncated by HL :7 atoms, but without any electrostatic interactions. Finally, !22389,435 ,485 is the classical energy of all atoms with CL atoms and with the charges of the two QM regions set to zero (to avoid double counting of the electrostatic interactions). However, the electrostatic interactions HL HL between the two QM systems is still accounted for in both !QM1a+ptch1b2 and !QM1b+ptch1a2 . This problem was solved by scaling down the point charges of the other QM system by a factor of 2 in both terms.77 The point-charge model of each QM system was obtained by a fit to the electrostatic potential (ESP) calculated for a wavefunction polarised by a point-charge model of the surroundings, but without the point charges when the ESP was calculated.77,78 The ESP points were sampled with the Merz–Kollman approach63 as implemented in Turbomole.38 The charges were updated in each step of the geometry optimisation. The geometry optimisations were continued until the energy change between two iterations was less than 2.6 J/mol (10–6 a.u.) and the maximum norm of the Cartesian gradients was below 10–3 a.u. The geometry optimisations were performed using the TPSS-D3 method39,47 and the def2-SV(P)43 basis set.

QM/MM-PBSA calculations The QM/MM-PBSA method77,78 can be viewed as a post-processing of QM/MM calculations to obtain more stable energies, including a solvation energy obtained by continuum methods. It is an adaptation of the widely used MM/PBSA approach79,80 for QM/MM calculations. In this approach, an approximation to the total free energy for each 8 ACS Paragon Plus Environment

Page 8 of 46

Page 9 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


state is obtained from: ME ; = !QM/MM-2QM + ;PB + ;np − AB

(4)

ME CL CL CL HL HL !QM/MM-2QM = !QM1a − !MM1a + !QM1b − !MM1b + !MM1a1b2

(5)

Here

HL HL where !QM1a and !QM1b are the QM energy of the two QM systems, truncated with HL atoms, CL CL but without any point-charge model. !MM1a and !MM1b are the corresponding MM energies of CL the two QM systems, but now with CL atoms and with full charges. Finally, !MM1a1b2 is the MM energy of the total system with CL atoms and full charges. This energy is similar to the QM/MM-2QM energy in Eqn. (1), but it employs mechanical embedding (ME), i.e. the electrostatic interactions between the QM and MM systems are calculated at the MM level. GSolv is the polar solvation energy, estimated by a continuum approach, obtained either by solving the Poisson–Boltzmann (PB) equation using the PB solver in Amber.59 Gnp is the nonpolar solvation energy, estimated from the solvent-accessible surface area (SASA), through the linear relation Gnp = 0.0227 SASA (Å2) + 3.85 kJ/mol.80,81 For the entropy term, we included only the entropy of the two QM systems, obtained by QM-cluster calculations at the TPSS-D3/def2-SV(P) level of the isolated QM systems, whereas the contribution from the surroundings was ignored, because such entropy effects are expected to be quite small in proton-transfer reactions.77 Two variants of QM/MM-PBSA-2QM were used. They differ in how the QM energy and the charges on the QM atoms are calculated.78,82,83 In the first approach, the QM energy and the charges were calculated from a vacuum wavefunction. In the second approach, the QM wavefunction is obtained with a point-charge model of the other systems. However, the final energy and the QM charges are calculated without this point-charge model and without reoptimising the wavefunction (obtained by setting the number of SCF iterations to 1 and using the final DFT integration grid in this single step). Thereby, the QM system is polarised by the MM surroundings and the cost of the polarisation is included in the calculations. The QM/MM-PBSA-2QM calculations were based on the spherical QM/MM-2QM structures. The calculations were automatized and performed by Linux shell scripts, which are available from the authors on request. Further details of the calculations can be found in http://www.teokem.lu.se/~ulf/Methods/qmmm_pbsa.html.

QTCP calculations The QTCP approach (QM/MM thermodynamic cycle perturbation) is a method to calculate free energies between two states, A and B, with a high-level QM/MM method, using free-energy perturbation (FEP) with sampling at only the MM level.34,84,85 It employs the thermodynamic cycle shown in Figure 3a, showing that the QM/MM free-energy difference between A and B can be obtained from three calculations: a FEP from A to B at the MM level, a FEP in methods space from MM to QM/MM for the A state and a similar calculation for the B state: ∆;CΜ/ΜΜA → B = Δ;22 A → B − Δ;22→C2/22 A + Δ;22→C2/22 B

(6)

In this study, we used QTCP to calculate the energy difference between two protonation states of the homocitrate ligand of the MoFe cluster calculated. To make the energies more 9 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

stable, we did not simply delete the proton. Instead, the proton was assumed to move to the Glu-212A residue on the opposite side of the protein (the considered MoFe cluster is in the C subunit). Our previous studies have shown that the results (in particular differences in pKa values between different state of the active site) do not depend much on which residue is used, as long as it is solvent exposed and does not form any H-bond with any other protein residue.86 Experimental studies have shown that solvent-exposed carboxylate groups have essentially the same pKa values of 4.0–4.4.87 The calculations were further stabilised by neutralising all solvent-exposed groups, as is often done in pKa calculations, because they are experimentally known not to affect the properties of buried active sites.88–90 The QTCP calculations were performed as has been described before:82,84 First, each state of interest was optimised by QM/MM-2QM, keeping system 2 fixed at the crystal structure. The first QM system (1a) contained the MoFe cluster, as shown in Figure 1a, and the other (1b) contained Glu-212A, modelled by an acetate group. For both QM systems, the TPSSD3/def2-SV(P) level of theory was used. The protein was then further solvated in an octahedral box of TIP3P water molecule, extending at least 9 Å from the QM/MM system (~186 100 atoms in total). For the A state (homocitrate protonated and Glu-212A deprotonated), the system was subjected to a 1000step minimisation, keeping atoms in the QM region fixed and restraining all heavy atoms towards the crystal structure with a force constant of 418 kJ/mol/Å2. Next, two 20 ps MD simulations were run with the heavy atoms still restrained. The first was run with a constant volume, the second with a constant pressure. Next, the size of the periodic box was equilibrated by a 100-ps MD simulation with a constant pressure and only the heavy atoms in QM system restrained to the QM/MM structure. The final structure of this simulation was used as the starting structure also for the other state (B, i.e. homocitrate deprotonated and Glu212A protonated). Finally, an equilibration of 200 ps and a production simulation of 400 ps were run with a constant volume for each state and with the QM systems kept fixed at the QM/MM-2QM structure. During the production run, 200 snapshots were collected every 2 ps. Based on these snapshots, a number of FEPs were performed as is shown in Figure 3b. First, FEPs were performed at the MM level between the two states A and B by changing the charges and coordinates of the QM system to those of the QM/MM calculations.82 We started from a protonated homocitrate and a deprotonated Glu-212A residue. Then, we introduced a dummy atom on Glu-212A, calculating its contributions to the free energy analytically from local configurational integrals.77,91,92 Next, the van der Waals parameters of this appearing proton were perturbed from those of a dummy atom with a zeroed a theoretical depth to polar hydrogen by a single-step perturbation. After that, the charges of the two QM sites were perturbed from those of the starting state to those of the final state (and zero charges for the disappearing proton). Such a reaction does not change the net charge of the protein, thereby making the calculated energies more stable. We performed nine simulations with the charges of the QM systems linearly transformed from the reactant state to the product state by a coupling parameter, λ, assuming values of 0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, and 1. This was accomplished by simply changing the charges in the parameter file, whereas the coordinates of the QM system were kept to those of the reactant state. Subsequently, the van der Waals parameters of the disappearing proton of the TNC were perturbed from those of a polar hydrogen to a dummy atom with a zeroed well depth in a single step. Next, the resulting non-interacting dummy atom was deleted, again calculating its contributions to the free energy analytically. In fact, the two analytical contributions will cancel in the net result. Finally, the coordinates of the two QM systems were transformed from those of the starting state (without the dissociating proton) to those of the final state. This was accomplished in five FEP simulations using another coupling parameter attaining the values of 0, 0.25, 0.5, 0.75, and 1. 10 ACS Paragon Plus Environment

Page 10 of 46

Page 11 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Two methods were used to estimate long-range electrostatic effects outside the simulated system. The results presented in the text and the tables were obtained with periodic systems and Ewald summation (as in the MD simulations). However, non-periodic calculations with an infinite cut-off and solvation energies outside the simulated system calculated by the Onsager models77,93 gave results that differed by only a few pKa units, but with a slightly larger hysteresis. For the two end states (A and B), we also performed a perturbation from MM to QM/MM-2QM energies, as has been described in detail before.34,84 These calculations were performed on the same octahedral systems as all the other QTCP-2QM calculations, but without any periodicity or Ewald summation, but instead with an infinite cutoff. The uncertainty of each step in Figure 3b was estimated as the standard error over the 200 MD snapshots. The net precision of the QTCP calculations was then obtained by error propagation over all the steps. All FEP calculations were performed with the local software calcqtcp. Further details of the QTCP calculations can be found in http://www.teokem.lu.se/~ulf/Methods/qtcp.html.

Quantum refinement Quantum refinement is standard crystallographic refinement supplemented by quantum chemical calculations for a small, but interesting part of the protein.35,94 Crystallographic refinement programs change the protein model (coordinates, occupancies, B factors, etc.) to improve the fit of the observed and calculated structure-factor amplitudes (usually estimated as the residual disagreement, the R factor). Owing to the limited resolution normally obtained for biomolecules, the experimental data are supplemented by some chemical information, usually in the form of a MM force field.95 Then, the refinement takes the form of a minimization or simulated annealing calculation by molecular dynamics using an energy function of the form Ecryst = wA EXray + EMM

(7)

where EXray is a penalty function that describes how well the model agrees with the experimental data (we have used a maximum-likelihood refinement target using amplitudes, MLF96,97). EMM is an MM energy function with bond, angle, dihedral and non-bonded terms, and wA is a weight factor, which is necessary because EMM is in energy units whereas EXray is unit-less. It determines the relative importance of the crystallographic raw data and the MM force field for the final structure. Quantum chemistry can be introduced in this function by replacing the MM potential for a small region of the protein (system 1) by a quantum mechanics (QM) calculation (in analogy to the QM/MM calculations), yielding a QM energy for system 1, EQM1. To avoid double counting we must then subtract the MM energy of system 1, EMM1, Ecqx = wA EXray + EMM12 + wQM EQM1 – EMM1

(8)

Thereby, we introduce an accurate energy function for the system of interest. Such an energy function is implemented in the software ComQum-X,35 which is a combination of the software Turbomole38 and the crystallography and NMR system (CNS), 98,99 version 1.3. The factor wQM in Eq. 5 is another weight, which is needed because the CNS MM force field is based on a statistical analysis of crystal structures.100 Therefore, the force constants are not energy-derived, as is the QM term, but they are in arbitrary statistical units. Experience has shown that the CNS force constants are typically three times larger than energy-based force 11 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

constants,100 and wQM = 3 has therefore been used throughout this work.35 Crystallographic refinement is traditionally performed without any electrostatic interactions, because hydrogen atoms are not discerned in the structure. We followed this custom and excluded electrostatics and hydrogen atoms from all crystallography and MM calculations (but they are of course included in the QM calculations). In analogy with the QM/MM calculations, the QM system was truncated by H atoms. The quantum-refinement calculations were based on the same crystal structure as all the other calculations in this paper.8 Coordinates, occupancies, B factors, and structure factors were obtained from the Protein Data Bank files 3U7Q. From these files, we also obtained the space group, unit-cell parameters, resolution limits, R factors and the test set used for the evaluation of the Rfree factor. The full protein was used in all calculations, including all crystal water molecules. In each cycle of the geometry optimization, the surrounding protein was allowed to relax by one cycle of crystallographic minimization and one cycle of individual B-factor refinement. However, the new coordinates and B factors were accepted only if the R factor was reduced. For the protein, we used the standard CNS force field (protein_rep.param, water_rep.param, and ion.param). However, CNS does not support anisotropic B factors, so these were ignored. The MM force field for non-standard residues were downloaded from the hetero-compound information centre Uppsala.101 The wA factor was determined by CNS to 0.0793. After quantum-refinement, anisotropic B factor and occupancy refinement was performed using phenix.refine102. Electron density maps were generated using phenix.maps. The QM calculations were performed at the TPSS/def2-SV(P) level of theory. Test calculations were also performed at the B3LYP/SV(P) level to check that they are consistent. The quality of the models was compared using the real-space difference density Zscore103 (RSZD), calculated by EDSTATS, which measures the local accuracy of the model. The maximum of the absolute values of the positive and negative RSZD (combined RSZD) for the homocitrate residue was taken as the quality metric. In all calculations, the maximum value was positive. RSZD is typically < 3.0 for a good model.

Result and Discussion The properties of nitrogenase (and most other enzymes) strongly depend on electrostatic interactions, both within the MoFe cluster and between the cluster and the surrounding protein. Such interactions crucially depend on the protonation state of the various residues, i.e. whether they are charged or not. Around the MoFe cluster, there are several residues that have an unclear protonation state. Of particular interest is the homocitrate molecule, which is a ligand of the Mo ion and therefore directly affects the electronic structure of the cluster. We have therefore employed six methods of varying sophistication to try to determine the protonation state of this ligand: QM-cluster, QM/MM, QM/MM-PBSA and QTCP calculations, as well as MD simulations and quantum-refinement calculations. There are eight additional protein residues close to the MoFe cluster, for which our protonation-state assignment was not fully conclusive, four His residues (His-195, 274, 362 and 451), three Glu residues (Glu-153, 380 and 440) and one Asp residue (Asp-445). We studied these by performing MD simulations of various protonation state and then comparing the H-bond pattern and root-mean-squared deviation (RMSD) from the starting crystal structure, as has been done before for His residues in three proteins.104 We performed one simulation with the preferred protonation state for all residues and one simulation for the alternative protonation state of each of the eight residues and three alternative protonation states of homocitrate. This was performed in an iterative manner until we reached a consensus result (i.e. until the reference calculation involved the preferred protonation state of all 12 ACS Paragon Plus Environment

Page 12 of 46

Page 13 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


residues). We describe only the results of the final calculations. We will describe the results for the eight residues, as well as for homocitrate, in separate sections. Then, we will describe the QM calculations for homocitrate.

His-195 In the crystal structure, His-195 forms a H-bond to the S2B atom of the MoFe cluster (the S and Fe atoms of the MoFe cluster are named as in the crystal structure; the names are shown in Figure 1a) by the NE2 atom (3.2 Å) and to a water molecule by the ND1 atom (2.8 Å; cf. Figure 4a). Thus, the NE2 atom must be protonated, but the protonation state of ND1 is unclear. This residue is often assumed to deliver protons to the MoFe cluster.1,25 We run MD simulations with the ND1 atom either deprotonated (called HIE) or protonated (called HIP). First, we studied H-bonds involving the sidechain of His-195 in the MD simulation: 1000 snapshots were collected from the 10 ns MD simulation and for each snapshot, all H-bonds were recorded (with the criteria that the H–O/N distance is less than 2.5 Å and the H–O/N–X angle is larger than 90º). In Table 1, these H-bonds are reported as the percentage of the snapshots in which the H-bond was observed (%). The philosophy is that if the protonation state is correct, the MD simulations will reproduce the H-bonds observed in the crystal structure. Nitrogenase is a dimer of two identical (heterodimeric) structures (consisting of subunits A and B, as well as C and D). Both were included in the MD simulations. Therefore, there are two copies of each residue (from subunits A and C) and we report the results from each of them individually. They give an indication of the reproducibility of the results and the flexibility of the protein around the studied group. In the HIE simulation, His-195 keeps the H-bonds observed in the crystal structure, i.e. to S2B and a water molecule with occurrences of 78–82% and 84–94%, respectively. In the HIP simulation, the same H-bonds are also observed (in 94–97% and 92–94% of the snapshots, respectively). However, there is also a H-bond between HD1 and the backbone O atom of Arg-277 in 35–47% of the snapshots, which is not observed in the crystal structure. Further information about which state is more realistic can be obtained by studying the root-mean-squared deviation (RMSD) of various atoms from the positions observed in the crystal structure. In this case, the philosophy is that if an incorrect protonation state is employed, the atoms in that residue or in nearby residues will move to release steric or electrostatic clashes or to form new favourable interactions. Therefore, the RMSD will be higher for incorrect protonation states. For each residue, we calculated the RMSD for three sets of atoms: For all atoms in the protein, for the studied residue alone, as well as for the studied residue and all residues with at least one atom within 3.5 Å of the studied residue. The first RMSD measures the global structure of the whole enzyme. The second RMSD shows whether the studied residue stays at the crystal position or not, whereas the third indicates whether any of the neighbouring residues moves to release any unfavourable interaction. For the latter two, one measure for each of the two subunits of the protein was obtained. The RMSD for the whole protein is similar in the two simulations (1.45 and 1.46 Å for HIE and HIP), which is expected for such a small change in a single residue. The same applies in all MD simulations (cf. Table S9), so the RMSD of the whole enzyme will not be further discussed (but the reference simulation has the lowest RMSD). However, the RMSD of the His-195 residue alone is much lower in the HIE simulation (0.59 and 0.69 Å in the two subunits) than in the HIP simulation (0.65 and 0.83 Å; cf. Figure 5). Likewise, if all residues within 3.5 Å of His-195 are considered, the RMSD calculation, it is lower in the HIE simulation than in the HIP simulation (0.92–0.99 Å compared to 1.09–1.13 Å). Thus, the results quite clearly point out that His-195 should be protonated only on the NE2 atom and not on ND1, at least in the resting state of the MoFe cluster. 13 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

His-274 His-274 is located ~7 Å from the MoFe cluster with mainly water molecules filling the gap. The ND1 atom forms a H-bond to a water molecule in the crystal structure (2.9 Å). The NE2 atom forms a H-bond to the backbone N atom of Phe-299 (3.0 Å), although the geometry is not fully ideal (Figure 4b). The residue is buried in the protein without any interaction with negatively charged groups, so it is unlikely that it is doubly protonated. Therefore, we run simulations for the HID and HIE states of this residue. In the HID simulations, the H-bonds observed in the crystal structure are seen in most of the snapshots (96–98% for HD1 and 69–83% for NE2). In the HIE simulation, ND1 still interacts with water (92–94%). However, the HE2 atom interacts instead with a water molecule in 96–97% of the snapshots. This is also reflected in the RMSD results. It is much smaller in the HID simulation (0.58–0.69 Å for His-274 alone and 0.75–0.81 Å with the 3.5-Å surrounding) than in the HIE simulations (2.43 Å for His-274 and 0.98–0.99 Å with the surroundings). Therefore, we can conclude that His-274 should undoubtedly be protonated only on the ND1 atom and not NE2.

His-362 His-362 is also located ~7 Å from the MoFe cluster, close to the surface of the protein. It is 3.8 Å from His-274, but without any direct H-bonds. In the crystal structure, the ND1 atom forms a H-bond with a water molecule (2.8 Å), whereas the NE2 atom interacts with the backbone O atom of Tyr-297 (2.9 Å; cf. Figure 4c). Consequently, the NE2 atom must be protonated, but the protonation status of the ND1 atom is unclear. Since the residue is rather close to the MoFe cluster, we performed MD simulations with the ND1 atom either deprotonated (HIE) or protonated (HIP). In the HIE simulation, the ND1 atom forms a H-bond with a water molecule in 74–75% of the snapshots. The HE2 atom forms the H-bond with the backbone O atom of Tyr-297 in 68–78% of the snapshots. In addition, it forms occasional H-bonds with the backbone O and the sidechain NE2 atoms of His-274 and with a water molecule (4–11% of the snapshots). In the HIP simulation, HD1 still interacts with one or occasionally two water molecules in essentially all snapshots. However, the HE2 atom forms a H-bond with Tyr-297 only in few snapshots (0–5%). Instead, it interacts with the NE2 atoms of His-274 and 451 (in 20–34% and 66–76% of the snapshots, respectively). Sometimes, it also forms a H-bond to a water molecule, (2–13% of the snapshots). The RMSD of His-362 is 0.74–0.76 Å in the HIE simulation, but very varying in the HIP simulation, 0.64 and 2.24 Å. However, if residues within 3.5 Å are included in the RMSD, the HIE simulation gave a smaller RMSD (0.76–0.79 Å) than the HIP simulation (1.42–1.45 Å). Thus, the HIE simulation gave the better and more stable structures.

His-451 His-451 is also located ~7 Å from the MoFe cluster with four water molecules filling the gap, and close to His-274 (~3.6 Å), although they do not interact directly. The NE2 atom forms a H-bond with a water molecule (2.8 Å). The ND1 atom is 2.9 Å from the carboxylate group of Asp-228 in the crystal structure, but the geometry is rather poor for a H-bond (Figure 4d). The carboxylate group of Asp-234 essentially stacks with the imidazole ring (~3.5 Å distance), without forming any H-bond. The close interaction with two carboxylate groups could indicate that the residue is positively charged, but the H-bonds are not as short and with an ideal geometry as you expect from ionic pairs. Therefore, we started MD simulations with His-451 either singly (HID) or doubly (HIP) protonated. 14 ACS Paragon Plus Environment

Page 14 of 46

Page 15 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In the HID simulation, the H1D atom forms a H-bond with OD1 or OD2 of Asp-228 in 75–84% of the snapshots. It also forms more occasional H-bonds with the backbone O atom of Asp-234 (7–22%) and the OG atom of Ser-238 (0–40%). NE2 forms a H-bond with water in 95–97% of the snapshots. In the HIP simulation, the two subunits of nitrogenase gave somewhat differing results: In the A subunit, H1D forms a H-bond with OD2 of Asp-228 in only 24% of the snapshots. Instead, it forms H-bonds with the backbone O atom of Asp-234 (45%) or the OG atom of Ser-238 (89%). However, in the C subunit, H1D forms a H-bond with OD2 of Asp-228 in 78% of the snapshots and H-bonds with the backbone O atom of Asp-234 or the OG atom of Ser-238 in 64 and 21% of the snapshots. The HE2 atom forms H-bonds with water in only 40% of the snapshots and instead interacts with NE2 of His-274 in 48% of the snapshots in the C subunit. In the A subunit, HE2 forms a H-bond to water molecules in 90% of the snapshots. The RMSD of His-451 is appreciably lower in the HID simulations (0.66 Å) than in the HIP simulations (0.82–0.93 Å). The same applies if the closest residues are included also in the comparison, although the results are more varying 0.73–0.78 Å for HID, compared to 0.68–0.94 Å for the HIP simulations. Thus, we conclude that the HID simulation gave the better and more stable results.

Glu-153 Glu-153 is buried inside the protein, with the carboxylate group ~5 Å from the P cluster (cf. Figure 4e). The OE atoms form H-bonds with a water molecule (2.7 Å) and the backbone O atom of Pro-85 (2.7 Å). The latter interaction indicates that the sidechain may be protonated. Therefore, we performed two simulations, one with the residue protonated (GLH) and one with it deprotonated (GLU). The two H-bonds observed in the crystal structure are preserved in the GLH simulation: 87–99% of the snapshots had the H-bond to water and 95–100% of the snapshots had the Hbond to Pro-85. Likewise, the H-bond to water is kept in essentially all the snapshots of the GLU simulation. However, the interaction with Pro-85 is completely broken. Instead, OE2 forms H-bonds with the backbone H atoms of Gly-84 and Glu-153 (97–100% and 95–96% of the snapshots, respectively). In addition, in subunit A, the OE1 also interacts with the HH atoms of Tyr-64 (91%) and HG of Ser-82 (60%). Consequently, there are also large differences in the RMSD results: The RMSD for Glu153 is appreciably lower in the GLH simulation (1.23–1.24 Å) than in the GLU simulation (1.54–1.67 Å). The same applies if all residues within 3.5 Å of Glu-153 are included in the comparison: The RMSD is 0.95–0.96 Å in the GLH simulation but 1.16–1.22 Å in the GLU simulation. Therefore, we can conclude that Glu-153 is protonated in nitrogenase.

Glu-380 Glu-380 is located ~6 Å from the Mo ion and it interacts with the homocitrate ligand via a water molecule. Therefore, the protonation state of this residue is potentially crucial for the active site. Glu-380 is located inside the protein, but there are several water molecules around it. One of the OE atoms interacts with two water molecules in the crystal structure (2.6 and 3.0 Å), one of which bridges to His-442 and the other to homocitrate (Figure 4f). The other OE atom forms a weak H-bond to the NE2 atom of Gln-191 (3.0 Å). Again, we have performed two MD simulations with Glu-380 either deprotonated (GLU) or protonated (GLH). In both MD simulations, there are considerable differences between the two subunits. In 15 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the GLU simulation of subunit C, the crystal interactions are retained: OE2 interacts with HE22 of Gln-191 in 63% of the snapshots and OE1 atom forms H-bonds with two water molecules. In addition, OE2 also interacts with two water molecules and occasionally with the backbone H atom of Gly-356 (20% of the snapshots). However, in the A subunit, the H-bonds to water are retained, but the interaction with Gln-191 is essentially broken (8%; in most snapshots one of the water molecules bridges to Gln-191). Instead, OE2 interacts with H of Gly356 and HE2 of His-442 (58 and 18% of the snapshots, respectively). Likewise, in the GLH simulation, the OE atoms interact with HE22 of Gln-191 in 88% of the snapshots in subunit A, but only in 18 % of the snapshots in subunit C. The OE1 atom interacts with a water molecule (10 and 71% in the two subunits). The HE2 atom forms a Hbond to a single water molecule in both subunits (98–99%), whereas the OE2 atom forms only occasional H-bonds (4–6% of the snapshots). This analysis shows that Glu-380 is quite flexible in both simulations. The interaction with Gln-191 is quite weak, whereas the interaction with two water molecules is best reproduced in the GLU simulation. Somewhat unexpectedly, the RMSD analysis shows that the GLU simulation of the A subunit (with the broken H-bond to Gln-191) gives the lowest RMSD (0.74 Å), whereas the two GLH simulations give intermediate values 1.01–1.04 Å, and subunit of C of the GLU simulation (with the intact H-bond to Gln-191) gives the highest RMSD of 1.18 Å. This may indicate that the H-bond to Gln-191 in the crystal structure is not ideal. If all atoms within 3.5 Å of Glu-380 is included in the RMSD, both GLU simulations give a lower RMSD (0.87– 0.88 Å) than the GLH simulations (1.03–1.09 Å). Therefore, we tend to prefer to keep Glu380 deprotonated, although this conclusion is not very strong.

Glu-440 and Asp-445 Glu-440 is located 9 Å from the Mo ion. It is connected to the homocitrate ligand via a water molecule. Like Glu-380, it is mainly buried, but with a number of water molecules around the carboxylate group. In the crystal structure, the OE1 atom forms a H-bond to a water molecule (2.8 Å, bridging to homocitrate) and to the OG atom of Ser-443 (2.8 Å; cf. Figure 4g). OE2 forms a rather weak H-bond to OH of Tyr-446 (3.0 Å) and a stronger H-bond to OD1 of Asp-445 (2.9 Å). The latter indicates that either Asp-445 or Glu-440 should be protonated. Asp-445 receives H-bonds from the NZ group of Lys-428 and the backbone N group of Glu-440 (both 2.8 Å), as well as from two water molecules (2.8 and 2.9 Å). Therefore, we performed two simulations: one with Glu-440 protonated and Asp-445 deprotonated (GLH/ASP) and another with Glu-440 deprotonated and Asp-445 protonated (GLU/ASH). In the GLH/ASP simulation, OE1, forms primarily H-bonds with water (71–137%), but also occasional interactions with HZ of Lys-428 (0–17%), HG of Ser-443 (0–40%) and HH of Tyr-446 (4–21%). The HE2 proton forms a H-bond with OD1 of Asp-445 in 26–99% of the snapshot, but also with a water molecule in 2–50% of the snapshots and with OH of Tyr-446 in subunit A. The OE2 atom forms more occasional interactions with water (25–29%), HZ of Lys-428 (12–20%), HG of Ser-443 (9–14%) and HH of Tyr-446 (0–5%). In the GLU/ASH simulation, both OE atoms show similar and varying interactions, indicating that the carboxylate group rotates slowly during the MD simulations. The OE atoms form many H-bonds, with several water molecules (170–265%, sum of both OE atoms), HZ of Lys-428 (4–13%), H of His-442 (54–71%), H of Ser-443 (80–97%), H of Asp445 (1–28%) and HH of Tyr-446 (4–21%). In addition, there is a H-bond between OE2 in Glu-440 and HE2 in Asp-445 in subunit C in 73% of the snapshots. For Asp-445 in the GLH/ASP simulation, the two OD atoms seem to interchange rather freely, except that only the OD2 atoms forms H-bonds with the HE and HH2 atoms of Arg16 ACS Paragon Plus Environment

Page 16 of 46

Page 17 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


439 (19–26%). The two OD atoms forms H-bonds with 2–3 water molecules (196–323%, sum of both atoms), with the HZ atoms of Lys-428 (2–91%) and with the backbone H atom of Glu-440 (123–160%). In the GLU/ASH simulation, OD1 of Asp-445 still interacts with water, but only in 61– 98% of the snapshots and with the backbone H atom of Glu-440 (20–99%). It also interacts with HZ of Lys-428, but only in 58% of the snapshots of subunit C. The HD2 atoms interacts with a water molecule in 10–102% of the snapshots and also with its own backbone O atom (0–22%), and with OE2 of Glu-440 in 73% of the snapshots of subunit C, as was already mentioned. The OD2 atom sometimes interacts with a water molecule (17–18%), with HZ of Lys-428 (6–66%), and with the backbone H atom of Glu-440 in 50% of the snapshots of subunit A. In summary, the H-bond analysis show that these two residues and the surroundings are very flexible. None of the simulations reproduce all H-bonds in the crystal structure. However, the GLH/ASP simulation, especially of subunit C, seems to come closest to the crystal structure (for four of the eight H-bonds), whereas the GLU/ASH simulation is better for only one (subunit A) or two (subunit C) of the H-bonds. The RMSD values are quite varying between the two subunits in the two simulations. However, the RMSD of both Glu-440 and Asp-445 are smaller in the GLH/ASP simulation (0.73–1.34 and 0.83–1.14 Å) than in the GLU/ASH simulation (1.72–2.00 and 1.31–1.39 Å). The same applies if all residues within 3.5 Å of each residue are included in the calculation, although the difference is smaller (0.79–1.05 and 0.81–1.08 Å, compared to 1.10 and 1.11– 1.13 Å). The RMSD is always lowest for subunit C in the GLH/ASP simulation, in agreement that this simulation also reproduces the H-bonds in the crystal structure best. Therefore, we prefer to let Glu-440 be protonated and Asp-445 deprotonated.

Protonation of homocitrate The homocitrate molecule is a bidentate ligand of the Mo ion in the MoFe cluster. Therefore, it is of key importance to decide its protonation status. It has three carboxylic groups which probably are deprotonated at pH 7 (the corresponding pKa values for citric acid are 3.1, 4.8, and 6.4).30 However, the situation is complicated by the fact that the molecule also contains a hydroxyl group, which coordinates directly to Mo (in addition to one of the carboxylate groups). This opens the possibility that the hydroxyl group is also deprotonated (giving a net –4 charge of the molecule) or that the hydroxylate group is deprotonated but one of the two non-coordinating carboxylate groups is protonated (giving a net –3 charge). Moreover, the homocitrate ligand is buried inside the protein, which could lead to further protonation of the carboxylate groups. In the protein, the hydroxyl O7 and carboxyl O5 atoms coordinate to Mo (the numbering of the atoms is shown in Figure 4h), both with a Mo–O distance of 2.2 Å. The O1 atom receives a H-bond from the sidechain NE2 atom of Gln-191 (2.9 Å) and the O3 atom receives a H-bond from the backbone N atom of Ile-425 (2.8 Å). In addition, there are a number of crystal-water molecules around the homocitrate ligand, which form H-bonds with O2 (two at 2.7–2.8 Å), O3 (two at 2.9 Å), O4 (three at 2.7–2.8), O5 (one at 2.8 Å) and O6 (two at 2.7–2.8 Å). There are two positively charged residues close to the homocitrate molecule, forming water-mediated H-bonds: NZ of Lys-426 interacts with O4 (2.8+2.8 Å), whereas NH2 of Arg105B interacts with both O3 (3.0+2.9 Å) and O6 (3.0+2.8 Å). Moreover, three negatively charged amino acids form water-mediated H-bonds with homocitrate: Glu-440 interacts with O3 (2.8+2.9 Å), Glu-427 with O4 (2.7+2.7 Å) and Glu-380 interacts also with O4 (2.6+2.8 Å). Thus, the surroundings do not compensate for the negative charge on the homocitrate ligand, but rather enhance it, which may indicate that homocitrate does not have a full –4 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

charge. In particular, it shows that there is a complicated network of both positively and negatively charged residues around the homocitrate ligand. After a thorough study of the H-bond network, we decided to test three different charge states of the homocitrate ligand. The O5/O6 carboxylate group binds to Mo and interacts indirectly with Arg-105B. Therefore, it is unlikely to be protonated. Likewise, the O3/O4 carboxylate group interacts indirectly with both Arg-105B and Lys-426 and forms three Hbonds each, as can be expected for a charged group. Therefore, this group was also assumed to be deprotonated. However, the O1/O2 carboxylate group forms only one or two H-bonds and it does not interact with any positively charged groups. Hence, it could very well be protonated. O1 receives a H-bond from Gln-191 in the syn-position and the anti-position is directed towards the alcoholic O7 atom with a O–O distance of 2.5 Å. The O2 atom forms Hbonds with two water molecules. Protonation in the anti-position would lead to rather poor Hbonding geometries, whereas protonation in the syn-position would give better geometries. Consequently, we decided to test four different protonation states, combining protonations in the syn-position of O2, and on O7 with the proton directed towards O1. Hence, we started four MD simulations with either both O2 and O7 protonated (2H, net charge –2), both deprotonated (0H; net charge –4), or either O2 (1Hc) or O7 (1Ha) protonated (“c” and “a” indicating protonation of the carboxylate or alcohol; net charge –3). These states are shown in Figure 2. The H-bonds observed in the four MD simulations are listed in Table 2. It can be seen that the H-bond between the O1 atom of homocitrate and HE21 of Gln191, present in the crystal structure, is observed only in 83% of the snapshots of the C subunit in the 2H simulation. In all the other simulations, O1 instead forms H-bonds to 1–3 water molecules. In most of the simulations, it also forms a H-bond to the HH atoms of Arg-105B, but in only 1– 72% of the snapshots. The O2 atom of homocitrate forms similar interactions, especially in the 1Ha and H0 simulations, in which the O1 and O2 atoms are identical carboxylate atoms that seem to interchange rapidly on the time scale of the MD simulations. In particular, O2 interacts with 1–3 water molecules, in accordance with the crystal structure, in which two Hbonds to water are seen. In the 2H and 1Hc simulations, O2 is protonated and the HO2 atom forms a H-bond to a water molecule in 59–91% of the snapshots, except in subunit C of the 2H simulation. The O2 atom also forms a H-bond to a water molecule in 16–101% of the snapshots. The O3 and O4 atoms reside on the same carboxylate group and seem to interchange rather freely during the MD simulations. Therefore, we will in this paragraph discuss the sum of the H-bond occurrences for these two atoms. In the crystal structure O3 makes a H-bond to the backbone H atom of Ile-425 and to two water molecules, whereas O4 forms H-bonds to three water molecules. In the MD simulations, the H-bond to Ile-425 is seen in 3–155% of the snapshots, less in the 2H simulation than in the other three simulations. The two subunits in the 0H simulation give strongly varying results, with only 3% occurrence of this H-bond in subunit A, but 155% in the other (indicating that it interacts with both O atoms in half of the snapshots). The two O atoms form H-bonds to 2–4 water molecules, somewhat more in the 2H simulation and less in the 0H simulation. In addition, they form H-bonds with HZ of Lys426 (1–97%, less in 1Hc), HG of Ser-101B (0–109%, especially in the A subunit of the 0H simulation), and HH of Arg-105B (6–258%, especially in the C subunit of the 2H simulation). In the 1Hc simulation, H-bonds with HG of the backbone H atom of Lys-426 are also observed in 4–32% of the snapshots. The O5 atom binds to the Mo ion both in the crystal structure and the MD simulations. In the crystal, it also forms a H-bond with a water molecule. This H-bond is retained in the MD simulations (90–119%). Likewise, the other O atom of that carboxylate group, O6, forms Hbonds to two water molecules both in the crystal structure and in the MD simulations, with 18 ACS Paragon Plus Environment

Page 18 of 46

Page 19 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


only small differences between the various simulations (177–211%). O7 is the alcohol oxygen in homocitrate and also a Mo ligand. In the crystal structure, it does not form any H-bonds with the surroundings. However, in the MD simulations, it forms a H-bond with 1–2 water molecules. In the H1c and H0 simulations (in which O7 is deprotonated), O7 interacts directly with water in 118–135% of the snapshots. In the other two simulation, in which O7 is protonated, it is instead mainly the HO7 atom that interacts with water (in 72–97 % of the snapshots). In addition, HO7 interacts with the O1 atom of homocitrate in some of the simulations (in up to 59% of the snapshots). In conclusion, the H-bond analysis shows that the homocitrate ligand is quite flexible in the MD simulations. In particular, we see only rather small differences between the various protonation states of homocitrate. Only a single subunit of the 2H simulation shows a stable H-bond between O1 and HE2 of Gln-191, whereas the H-bond between O3 and H of Ile-425 is somewhat better reproduced, except in the 2H simulation. H-bonds to water molecules seem to be well described in all simulations. There are indications that the MM force field may give somewhat too strong interactions between the carboxylate groups of homocitrate and the positively charged groups of Lys-426 and Arg-105B. Fortunately, the RMSD analysis gives more conclusive results. The RMSD of all atoms in the protein is lowest for the 1Ha simulation (1.45 Å), followed by 1Hc, 2H, and 0H, but the differences are small (1.47–1.49 Å; Table S8). However, the differences are larger for the RMSD of the homocitrate ligand itself, as can be seen in Figure 5: It is lowest for the 1Ha simulation (1.73–1.85 Å in the two subunits), followed by the 2H simulation (1.78–1.97 Å). The other two simulations give appreciably higher RMSDs, 2.24 Å for the 1Hc simulation, and 2.19–2.33 Å for the 0H simulation. If the RMSD is calculated for all residues within 3.5 Å of homocitrate, the same ordering is observed: The 1Ha simulation gives the lowest RMSD (0.88–0.93 Å), followed by the 2H (0.93–0.95 Å) and 1Hc simulation (0.98–1.11 Å), whereas the 0H simulation gives the largest RMSD (1.15–1.22 Å). Thus, the 1Ha state gives the lowest RMSD compared to the crystal structure.

QM-cluster calculations of homocitrate To obtain further support for the protonation state of the homocitrate residue, we performed also QM calculations at several levels of approximation. First, we performed QM cluster calculations in the COSMO continuum solvent with the standard QM model of the MoFe cluster in Figure 1a. We tested the four protonation states in Figure 2. The pKa was calculated according to Eqn. 2 and the results are collected in Table 3. It can be seen that the calculated pKa values are rather insensitive to the theoretical treatment: Optimising the structures in a vacuum or in continuum solvent changes the pKa values by less than 3 pKa units. Likewise, a change from the pure TPSS functional to the hybrid B3LYP functional changes the pKa values by less than 3 pKa units. Increasing the basis set from def2-SV(P) to the much larger def2-TZVPD (including diffuse functions on all atoms) has a slightly larger effect, up to 5 pKa units. On the other hand, it is important to include the thermal effects (including the entropy and zero-point energy), which reduce the pKa values by 5–7 pKa units. Changing the dielectric constant of the continuum solvent in the energy calculation has an even larger effect, 47 pKa units for the two pKa values involving 2H, and 65 pKa units for the two pKa values involving 0H. Considering that the homocitrate ligand is in a very polar surrounding with several both positively and negatively charged groups and many water molecules (e.g. 3 Asp/Glu, 8 Lys/Arg and ~38 waters within 9 Å), it is most likely that the results with a high dielectric constant (80) are more reliable. These results (with both DFT methods and basis sets) show that at pH 7, 1Ha is the most stable protonation state of homocitrate. However, it is only 19 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

slightly more stable than the 2H state, with a pKa(2H/1Ha) value of 6.1–6.4 (with the large basis set). The pKa values involving the fully deprotonated 0H state (17–23) show that it is highly unlikely. The two singly protonated states (1Ha and 1Hc) contain the same atoms and therefore give directly comparable energies. It can be seen that the 1Ha state is 30–32 kJ/mol more stable than the 1Hc state with both DFT methods and the large basis set. This energy difference is more reliable than the pKa calculations, because it does involve calculations with differing net charges and therefore does not depend on the dielectric constant (less than 0.1 kJ/mol difference between calculations with dielectric constants of 4 and 80). Interestingly, when the 1Ha state is optimised in the continuum solvent, the proton moves from the alcohol group (O7) atom, to the carboxylate O1 atom: In the vacuum-optimised structure, the O7–H and O1–H distances are 1.09 and 1.44 Å, whereas in the solventoptimised (ε = 80) structure, they are essentially inverted, 1.40 and 1.09 Å, respectively. The same happens when optimised with a dielectric constant of 4, but not for the 2H state, because the other O atom (O2) in the carboxylate group is then also protonated (the O7–H and O1–H distances are 1.67 and 1.01 Å when optimised with ε = 80). However, it can be seen from Table 3 that this movement of the H atom from O7 to O2 does not affect the pKa values by more than 1 pKa unit. Still, it shows that the binding of O2 to Mo changes its pKa values to below that of a free carboxylate group.

QM/MM calculations on homocitrate Next, we performed QM/MM-2QM calculations on the MoFe cluster in the full nitrogenase tetramer. The aim was to estimate the relative stability of the four protonation states of the homocitrate ligand (i.e. the pKa values). To make the calculations more stable, they were performed on an enzyme with all solvent-exposed charges neutralised.81–83 Moreover, we studied the transfer of a proton from homocitrate to a solvent-exposed Glu group on the opposite side of the protein (Glu-212A).86 Thereby, we avoid the change in the net charge of the protein and excessive electrostatic terms. All four states of homocitrate in Figure 2 were optimised with QM/MM-2QM, including the full MoFe cluster in the first QM system (Figure 1a) and an acetate model of Glu-212A in a second QM system. Then, singlepoint energies were calculated with the larger def2-TZVPD basis set and also with the B3LYP method. In addition, a better solvation was obtained with the QM/MM-PBSA approach, employing charges either from a vacuum or a polarised wavefunction. All results are collected in Table 4. It can be seen that the QM/MM calculations still predict that 1Ha state is the most stable state at pH 7. The effect of the basis set and DFT functional is still small (up to 4 pKa units). The 1Ha state is still more stable than the 1Hc state, but the energy difference is somewhat smaller (9–28 kJ/mol) than in the QM-cluster calculations and more dependent on the functional and basis set. There are two estimates of this energy difference, because separate QM/MM-2QM structures were obtained for the 2H/1H and 1H/0H calculations (to keep the net charge of simulated system neutral). The proton is still transferred from O7 to O1 during the geometry optimisation of the 1Ha state (O7–H and O1–H distances of 1.55 and 1.05 Å, respectively). The QM/MM-2QM calculations provide an atomistic account of the surrounding protein and solvent, but they are not fully relaxed for each change in the protonation state. A somewhat better picture of the solvation can be obtained with the QM/MM-PBSA approach, in which explicit water molecules are replaced by a Poisson–Boltzmann continuum solvent. Such calculations are also included in Table 4, obtained with charges calculated either for a vacuum (vac) or a polarised (ptch) wavefunction. It can be seen that such calculations 20 ACS Paragon Plus Environment

Page 20 of 46

Page 21 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


strongly favour the protonated states. With the vacuum charges, the 1Ha state is still most stable, but only by 1–2 pKa units. However, with the polarised charges, the 2H state is actually more stable by 2–5 pKa units. Unfortunately, the results strongly depend on the continuum solvation model – generalised Born solvation instead indicates that the 0H state is most stable. Finally, it can also be seen from Table 4 that the 1Ha state is strongly stabilised by the continuum solvation, so that it is 70–88 kJ/mol more stable than the 1Hc state. The most accurate results are obtained with the QTCP approach, in which we perform FEP calculations both between the two protonation states of interest at the MM level and from the MM to the QM/MM description for the two end states, as is illustrated in Figure 3. These calculations have the advantage that they include a full flexibility of the surrounding protein and solvent, allowing them to fully relax with respect to the change in the charge distribution when the protonation state changes. Moreover, they do not involve continuum-solvation calculations with an ill-defined dielectric constant. The results of these calculations are listed in Table 5. It can be seen that QTCP still predicts that 1Ha is the most stable state. It is 10–27±12 kJ/mol more stable than the 1Hc state (after corrections for the thermal effects and for a change of the theoretical method to B3LYP-D3/def2-TZVPD). Moreover, the pKa(2H/1Ha) value is 3±1, indicating that 1Ha is the preferred protonation state at pH 7. Likewise, the pKa(1Ha/0H) value is 19±2, showing that the 1Ha state is preferred before 0H at pH 7. The uncertainty of the calculations is 1.3–1.7 pKa units, coming primarily from the charge and MM→QM/MM perturbations with approximately equal contributions. The hysteresis of the calculations is 2±3 pKa units, confirming the estimated uncertainties.

Quantum refinement Finally, we investigated the protonation state of the homocitrate ligand with an additional independent method, quantum refinement.35 In this approach, we employ the raw data (the structure factors) of the starting crystal structure8 and re-refine it with the four investigated protonation states of homocitrate, replacing the empirical MM potential used in the original refinement with QM/MM calculations of the MoFe cluster, which provide an ideal structure of the active site for a given protonation state (or at least a considerably better structure than the MM potential). We have previously shown that quantum refinement allows us to locally improve medium-resolution crystal structures and to determine the protonation state of metalbound water molecules in protein crystal structures.94,105,106 The results of this approach are shown in Table 6. It can be seen that the four structures give similar R and Rfree values, 0.15 and 0.16–0.17. This is expected, because they are global quality factors. Therefore, we instead used the local RSZD score, which measures the accuracy of the model through an analysis of difference density. As can be seen in Table 6, the homocitrate RSZD score shows extensive differences between the four states and it is lowest (and therefore best) for the 1Ha state (3.0). It is slightly higher for the 2H state (3.2) and much higher for the 1Hc and 0H states (8 and 10). This is also reflected in the electrondensity maps shown in Figure 6. In particular, it can be seen that the O1 atom is too far from the O7 atom in the 2H state (too little density between the two atoms, shown as a green volume in the Fo – Fc difference map in Figure 6b and too much density on the opposite side of the O1 atom, shown as a red volume; such volumes are not present for the 1Ha state in Figure 6a). This reflects a change in the O1–O7 distance from 2.44 Å in the 1Ha structure to 2.56 Å in the 2H structure, as an effect of the weakened H-bond (1.32 and 1.64 Å, respectively). This is the prime difference between the two protonation states and it is very satisfying that the electron-density maps so clearly reflect this difference. Calculations at the B3LYP/def2-SV(P) level of theory gave very similar results (Table S12 and Figure S1). 21 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Conclusions In this article, we have studied the protonation state of several groups in nitrogenase. In particular, we have studied the homocitrate group, which is a bidentate ligand of the Mo ion in the MoFe cluster. We have employed several different computational methods. For eight residues close to the MoFe and P clusters, we have used MD simulations to determine the most probably protonation state (besides the standard procedure of studying the surroundings, the solvent accessibility and the H-bond network around the residue). We have run one MD simulation for each alternative protonation state and then examined which of the states gives the lowest RMSD from the starting crystal structure, for the whole protein, for the studied residue and for all residues with at least one atom within 3.5 Å of the studied residue. This analysis gave clear and consistent results, as can be seen in Figure 5 and Table S9: His-195 and 362 are protonated on the NE2 atom and not on the N1D atom, His-274 and 451 are protonated on the N1D atom and not on NE2, Glu-153 is protonated, Glu-380 deprotonated, whereas Glu-440 is protonated with the HE2 atom directed towards the deprotonated Asp-445 residue. We have also studied the H-bond pattern during the MD simulations. For the four His residues and Glu-153, this confirmed the above assignment of the protonation state. However, for the other three carboxylate residues, the MD simulations showed that the sites are flexible, often with large differences between the two halves of the protein. Therefore, it was hard to decide which of the two simulations gave the better results. In particular, there was a tendency of the MD simulations to give additional H-bonds with water molecules (which may not be seen in the crystal structures owing to the dynamics) and to charged protein residues. The situation was similar for homocitrate. Based on the H-bond pattern in the crystal structure, we tested the four protonation states in Figure 2. In the MD simulations, they showed large flexibility with no clear conclusions from the H-bond analysis. However, the RMSD results indicated that the 1Ha state is most favourable. Owing to the proximity of the homocitrate group to the MoFe cluster, we also performed a set of QM calculations of the four protonation states. QM-cluster calculations (Table 3) again pointed out 1Ha as the most stable state at pH 7, but it is less than one pKa unit more stable than the 2H state. Moreover, the results depended strongly on the dielectric constant assumed for the COSMO continuum solvent. QM/MM calculations also supported 1Ha as the most stable state, but QM/MM-PBSA gave more varying results, owing to the strong solvation effects. However, the more accurate QTCP calculations confirmed that the 1Ha state is most stable, four pKa units more stable than the 2H state. These calculations clearly show how hard it is to estimate pKa values of groups in proteins, especially inside the protein and close to metal ions.107–109 In particular, continuum methods are problematic because the results strongly depend on the dielectric constant of the protein, which is unknown and may vary in different parts of the protein. Moreover, dynamic effects may be significant. Therefore, we prefer to trust the results obtained with atomistic simulations, like QTCP, which includes explicit dynamics and avoids the use of a protein dielectric constant. Finally, we have also performed quantum refinement of the original crystal structure with the four protonation states of homocitrate. These calculations again showed that the 1Ha protonation state fits the crystallographic raw data best, slightly better than the 2H state. In particular, the difference-density maps in Figure 6b convincingly show that the O1–O7 distance is too large in the 2H structure. Consequently, we can with quite some confidence conclude that homocitrate most likely is singly protonated (i.e. with a net charge of –3) in nitrogenase at pH 7. This state has been used only in one single previous study.13 The proton is shared between the alcohol O7 and the carboxylate O1 atoms and both the QM-cluster and 22 ACS Paragon Plus Environment

Page 22 of 46

Page 23 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the QM/MM calculations indicate that it actually resides on the O1 atom with a short H-bond to O7 (1.39–1.55 Å). This provides invaluable information for future QM/MM studies of the reaction mechanism of nitrogenase.

Supporting Information MM charges, van der Waals parameters and restraints used for the metal clusters, list of buried charges in the protein, RMSD of the entire protein in the various simulations, more detailed account of the hydrogen bonds around the studied residues in the MD simulations

Acknowledgements This investigation has been supported by grants from the Swedish research council (project 2014-5540), from COST through Action CM1305 (ECOSTBio) and a scholarship to LC from the China Scholarship Council. The computations were performed on computer resources provided by the Swedish National Infrastructure for Computing (SNIC) at Lunarc at Lund University.

References (1) Hoffman, B. M.; Lukoyanov, D.; Yang, Z.-Y.; Dean, D. R.; Seefeldt, L. C. Mechanism of Nitrogen Fixation by Nitrogenase: The Next Stage. Chem. Rev. 2014, 114, 4041– 4062. (2) Smil, V. Enriching the Earth: Fritz Haber, Carl Bosch, and the Transformation of World Food Production; MIT Press: Cambridge, MA, 2004. (3) Smith, B. E. Nitrogenase Reveals Its Inner Secrets. Science 2002, 297, 1654–1655. (4) Burgess, B. K.; Lowe, D. J. Mechanism of Molybdenum Nitrogenase. Chem. Rev. 1996, 96, 2983–3012. (5) Schmid, B.; Chiu, H.-J.; Ramakrishnan, V.; Howard, J. B.; Rees, D. C. Nitrogenase. In Handbook of Metalloproteins; John Wiley & Sons, Ltd, 2006; pp 1025–1036. (6) Kim, J.; Rees, D. C. Structural Models for the Metal Centers in the Nitrogenase Molybdenum-Iron Protein. Science 1992, 257, 1677–1682. (7) Einsle, O.; Tezcan, F. A.; Andrade, S. L. A.; Schmid, B.; Yoshida, M.; Howard, J. B.; Rees, D. C. Nitrogenase MoFe-Protein at 1.16 Å Resolution: A Central Ligand in the FeMo-Cofactor. Science 2002, 297, 1696. (8) Spatzal, T.; Aksoyoglu, M.; Zhang, L.; Andrade, S. L. A.; Schleicher, E.; Weber, S.; Rees, D. C.; Einsle, O. Evidence for Interstitial Carbon in Nitrogenase FeMo Cofactor. Science 2011, 334, 940–940. (9) Spatzal, T.; Perez, K. A.; Einsle, O.; Howard, J. B.; Rees, D. C. Ligand Binding to the FeMo-Cofactor: Structures of CO-Bound and Reactivated Nitrogenase. Science 2014, 345, 1620–1623. (10) Einsle, O. Nitrogenase FeMo Cofactor: An Atomic Structure in Three Simple Steps. J. Biol. Inorg. Chem. 2014, 19, 737–745. (11) Lancaster, K. M.; Roemelt, M.; Ettenhuber, P.; Hu, Y.; Ribbe, M. W.; Neese, F.; Bergmann, U.; DeBeer, S. X-Ray Emission Spectroscopy Evidences a Central Carbon in the Nitrogenase Iron-Molybdenum Cofactor. Science 2011, 334, 974–977. (12) Eady, R. R. Structure−Function Relationships of Alternative Nitrogenases. Chem. Rev. 1996, 96, 3013–3030. (13) Bjornsson, R.; Lima, F. A.; Spatzal, T.; Weyhermüller, T.; Glatzel, P.; Bill, E.; Einsle, O.; Neese, F.; DeBeer, S. Identification of a Spin-Coupled Mo(III) in the Nitrogenase 23 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(14) (15)

(16) (17)

(18) (19)

(20)

(21) (22) (23) (24)

(25)

(26)

(27)

(28) (29)

(30) (31) (32) (33)

Iron–molybdenum Cofactor. Chem. Sci. 2014, 5, 3096–3103. Thorneley, R. N. F.; Lowe, D. J. No Title. In Molybdenum Enzymes; Spiro, T. G., Ed.; Wiley: New York, 1985; pp 221–284. Tuczek, F. Nitrogen Fixation in Nitrogenase and Related Small-Molecule Models: Results of DFT Calculations. In RSC Metallobiology Series 7; Hille, R., Schulzke, C., Kirk, M. L., Eds.; Royal Society of Chemistry: Cambridge, 2017; pp 223–274. Dance, I. Theoretical Investigations of the Mechanism of Biological Nitrogen Fixation at the FeMo Cluster Site. J. Biol. Inorg. Chem. 1996, , 581–586. Stavrev, K. K.; Zerner, M. C. Studies on the Hydrogenation Steps of the Nitrogen Molecule at the Azotobacter Vinelandii Nitrogenase Site. Int. J. Quantum Chem. 1998, 70, 1159–1168. Siegbahn, P. E. M.; Westerberg, J.; Svensson, M.; Crabtree, R. H. Nitrogen Fixation by Nitrogenases: A Quantum Chemical Study. J. Phys. Chem. B 1998, 102, 1615–1623. Lovell, T.; Li, J.; Liu, T.; Case, D. A.; Noodleman, L. FeMo Cofactor of Nitrogenase: A Density Functional Study of States MN, Mox, MR, and MI. J. Am. Chem. Soc. 2001, 123, 12392–12410. Xie, H.; Wu, R.; Zhou, Z.; Cao, Z. Exploring the Interstitial Atom in the FeMo Cofactor of Nitrogenase: Insights from QM and QM/MM Calculations. J. Phys. Chem. B 2008, 112, 11435–11439. Kästner, J.; Blöchl, P. E. Ammonia Production at the FeMo Cofactor of Nitrogenase: Results from Density Functional Theory. J. Am. Chem. Soc. 2007, 129, 2998–3006. Hallmen, P. P.; Kästner, J. N2 Binding to the FeMo-Cofactor of Nitrogenase. Zeitschrift fur Anorg. und Allg. Chemie 2015, 641, 118–122. Dance, I. Activation of N2, the Enzymatic Way. Zeitschrift für Anorg. und Allg. Chemie 2015, 641, 91–99. Varley, J. B.; Wang, Y.; Chan, K.; Studt, F.; Nørskov, J. K. Mechanistic Insights into Nitrogen Fixation by Nitrogenase Enzymes. Phys. Chem. Chem. Phys. 2015, 17, 29541–29547. Siegbahn, P. E. M. Model Calculations Suggest That the Central Carbon in the FeMoCofactor of Nitrogenase Becomes Protonated in the Process of Nitrogen Fixation. J. Am. Chem. Soc. 2016, 138, 10485–10495. McKee, M. L. A New Nitrogenase Mechanism Using a CFe 8 S 9 Model: Does H 2 Elimination Activate the Complex to N2 Addition to the Central Carbon Atom? J. Phys. Chem. A 2016, 120, 754–764. Rao, L.; Xu, X.; Adamo, C. Theoretical Investigation on the Role of the Central Carbon Atom and Close Protein Environment on the Nitrogen Reduction in Mo Nitrogenase. ACS Catal. 2016, 6, 1567–1577. Dance, I. Activation of N2, the Enzymatic Way. Zeitschrift für Anorg. und Allg. Chemie 2015, 641, 91–99. McKee, M. L. A New Nitrogenase Mechanism Using a CFe 8 S 9 Model: Does H 2 Elimination Activate the Complex to N 2 Addition to the Central Carbon Atom? J. Phys. Chem. A 2016, 120, 754–764. Weast, R. C. CRC Handbook of Chemistry and Physics, 54th ed.; Weast, R. C., Ed.; CRC Press: Cleavland, Ohio, 1974. Chang, C. M.; Wang, M. K. Linear Relationship for Acidity and Stability in Hexaaqua Metal Ions — Density Functional Studies. Chem. Phys. Lett. 1998, 286, 46–50. Femo-cofactor, N.; Harris, T. V; Szilagyi, R. K. Comparative Assessment of the Composition and Charge State of. Inorg. Chem 2011, 4811–4824. Hu, L.; Söderhjelm, P.; Ryde, U. Accurate Reaction Energies in Proteins Obtained by Combining QM/MM and Large QM Calculations. J. Chem. Theory Comput. 2013, 9, 24 ACS Paragon Plus Environment

Page 24 of 46

Page 25 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(34) (35) (36)

(37) (38) (39)

(40) (41) (42) (43) (44)

(45) (46)

(47)

(48) (49) (50)

(51) (52)

(53)

640–649. Rod, T. H.; Ryde, U. Quantum Mechanical Free Energy Barrier for an Enzymatic Reaction. Phys. Rev. Lett. 2005, 94, 1–4. Ryde, U.; Olsen, L.; Nilsson, K. Quantum Chemical Geometry Optimizations in Proteins Using Crystallographic Raw Data. J. Comput. Chem. 2002, 23, 1058–1070. Olsson, M. H. M.; Søndergaard, C. R.; Rostkowski, M.; Jensen, J. H. PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa Predictions. J. Chem. Theory Comput. 2011, 7 (2), 525–537. Schrödinger, L. Schrödinger Release 2015-4: Maestro. Schrödinger, LLC: New York, NY 2016. Furche, F.; Ahlrichs, R.; Hättig, C.; Klopper, W.; Sierka, M.; Weigend, F. Turbomole. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4, 91–100. Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Non-Empirical Meta-Generalized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. Lett. 2003, 91, 146401. Becke, A. D. Density-Functional Exchange-Energy Approximation With Correct Asymptotic-Behavior. Phys. Rev. A 1988, 38, 3098–3100. Lee, C.; Yang, W.; Parr, R. G. Into a Functional of the Electron Density F F. Phys. Rev. B 1988, 37, 785–789. Becke, A. D. A New Mixing of Hartree–Fock and Local Density-Functional Theories. J. Chem. Phys. 1993, 98, 1372. Schäfer, A.; Horn, H.; Ahlrichs, R. Fully Optimized Contracted Gaussian Basis Sets for Atoms Li to Kr. J. Chem. Phys. 1992, 97, 2571–2577. Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305. Eichkorn, K.; Treutler, O.; Ohm, H.; Haser, M.; Ahlrichs, R. Auxiliary Basis-Sets to Approximate Coulomb Potentials . Chem. Phys. Lett. 1995, 240, 283–289. Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R. Auxiliary Basis Sets for Main Row Atoms and Transition Metals and Their Use to Approximate Coulomb Potentials. Theor. Chem. Acc. 1997, 97, 119–124. Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104 (19 pages). Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456–1465. Szilagyi, R. K.; Winslow, M. A. On the Accuracy of Density Functional Theory for Iron—sulfur Clusters. J. Comput. Chem. 2006, 27, 1385–1397. Greco, C.; Fantucci, P.; Ryde, U.; Gioia, L. de. Fast Generation of Broken-Symmetry States in a Large System Including Multiple Iron–sulfur Assemblies: Investigation of QM/MM Energies, Clusters Charges, and Spin Populations. Int. J. Quantum Chem. 2011, 111, 3949–3960. Cao, L.; Ryde, U. Influence of the Protein and DFT Methods on the Broken-Symmetry and Spin States in Nitrogenase. Chem. - A Eur. J. 2017, submitted. Lukoyanov, D.; Pelmenschikov, V.; Maeser, N.; Laryukhin, M.; Yang, T. C.; Noodleman, L.; Dean, D. R.; Case, D. A.; Seefeldt, L. C.; Hoffman, B. M. Testing If the Interstitial Atom, X, of the Nitrogenase Molybdenum−Iron Cofactor Is N or C: ENDOR, ESEEM, and DFT Studies of the S ) 3/2 Resting State in Multiple Environments. Inorg. Chem. 2007, 46, 11437–11449. Klamt, A.; Schuurmann, G. Cosmo - a New Approach To Dielectric Screening in 25 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(54)

(55) (56) (57)

(58) (59)

(60)

(61)

(62)

(63) (64) (65) (66) (67)

(68)

(69) (70)

(71) (72)

Solvents With Explicit Expressions for the Screening Energy and Its Gradient. J. Chem. Soc. Perkin Trans. 2 1993, 799–805. Schäfer, A.; Klamt, A.; Sattel, D.; Lohrenz, J. C. W.; Eckert, F. COSMO Implementation in TURBOMOLE: Extension of an Efficient Quantum Chemical Code towards Liquid Systems. Phys. Chem. Chem. Phys. 2000, 2, 2187–2193. Klamt, A.; Jonas, V.; Bürger, T.; Lohrenz, J. C. W. Refinement and Parametrization of COSMO-RS. J. Phys. Chem. A 1998, 102, 5074–5085. Sigfridsson, E.; Ryde, U. Comparison of Methods for Deriving Atomic Charges from the Electrostatic Potential and Moments. J. Comput. Chem. 1998, 1, 377–395. Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Aqueous Solvation Free Energies of Ions and Ion-Water Clusters Based on an Accurate Value for the Absolute Aqueous Solvation Free Energy of the Proton. J. Phys. Chem. B 2006, 110, 16066–16081. Jensen, F. Introduction to Computational Chemistry, 3rd ed.; John Wiley & Sons, Ltd: Chichester, 2017. Case, D. A.; Berryman, J. T.; Betz, R. M.; Cerutti, D. S.; Cheatham, T. E.; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.; et al. AMBER 14. University of California: San Francisco 2014. Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696–3713. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926–935. Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. A Well-Behaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97, 10269–10280. Besler, B. H.; Merz, K. M.; Kollman, P. A. Atomic Charges Derived from Semiempirical Methods. J. Comput. Chem. 1990, 11, 431–439. Bartolotti, L. J.; Pedersen, L. G.; Charifson, P. S. Long Range Nonbonded Attractive Constants for Some Charged Atoms. J. Comput. Chem. 1991, 12, 1125–1128. Hu, L.; Ryde, U. Comparison of Methods to Obtain Force-Field Parameters for Metal Sites. J. Chem. Theory Comput. 2011, 7, 2452–2463. Seminario, J. M. Calculation of Intramolecular Force Fields from Second-Derivative Tensors. Int. J. Quantum Chem. 1996, 30, 1271–1277. Nilsson, K.; Lecerof, D.; Sigfridsson, E.; Ryde, U. An Automatic Method to Generate Force-Field Parameters for Hetero-Compounds. Acta Crystallogr. - Sect. D Biol. Crystallogr. 2003, 59, 274–289. Barney, B. M.; Mcclead, J.; Lukoyanov, D.; Laryukhin, M.; Yang, T.; Dean, D. R.; Hoffman, B. M.; Seefeldt, L. C. Diazene (HN=NH) is a Substrate for Nitrogenase : Insights into the Pathway of. Biochemistry 2007, 46, 6784–6794. Wu, X.; Brooks, B. R. Self-Guided Langevin Dynamics Simulation Method. Chem. Phys. Lett. 2003, 381, 512–518. Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684–3690. Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N -log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089–10092. Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327–341. 26 ACS Paragon Plus Environment

Page 26 of 46

Page 27 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(73)

Ryde, U. The Coordination of the Catalytic Zinc in Alcohol Dehydrogenase Studied by Combined Quantum-Chemical and Molecular Mechanics Calculations. J. Comput. Aided. Mol. Des. 1996, 10, 153–164. (74) Ryde, U.; Olsson, M. H. M. Structure, Strain, and Reorganization Energy of Blue Copper Models in the Protein. Internat. J. Quant. Chem. 2001, 81, 335–347. (75) Reuter, N.; Dejaegere, A.; Maigret, B.; Karplus, M. Frontier Bonds in QM/MM Methods: A Comparison of Different Approaches. J. Phys. Chem. A 2000, 104, 1720– 1735. (76) Hu, L.; Söderhjelm, P.; Ryde, U. On the Convergence of QM/MM Energies. J. Chem. Theory Comput. 2011, 7, 761–777. (77) Hu, L.; Farrokhnia, M.; Heimdal, J.; Shleev, S.; Rulíšek, L.; Ryde, U. Reorganization Energy for Internal Electron Transfer in Multicopper Oxidases. J. Phys. Chem. B 2011, 115, 13111–13126. (78) Kaukonen, M.; Söderhjelm, P.; Heimdal, J.; Ryde, U. QM/MM-PBSA Method to Estimate Free Energies for Reactions in Proteins. J. Phys. Chem. B 2008, 112, 12537– 12548. (79) Srinivasan, J.; Cheatham, T. E.; Cieplak, P.; Kollman, P. A.; Case, D. A. Continuum Solvent Studies of the Stability of DNA, RNA, and Phosphoramidate−DNA Helices. J. Am. Chem. Soc. 1998, 120, 9401–9409. (80) Genheden, S.; Ryde, U. The MM/PBSA and MM/GBSA Methods to Estimate LigandBinding Affinities. Expert Opin. Drug Discov. 2015, 441, 1–13. (81) Kuhn, B.; Kollman, P. A. Binding of a Diverse Set of Ligands to Avidin and Streptavidin: An Accurate Quantitative Prediction of Their Relative Affinities by a Combination of Molecular Mechanics and Continuum Solvent Models. J. Med. Chem. 2000, 43, 3786–3791. (82) Heimdal, J.; Kaukonen, M.; Srnec, M.; Rulíšek, L.; Ryde, U. Reduction Potentials and Acidity Constants of Mn Superoxide Dismutase Calculated by QM/MM Free-Energy Methods. ChemPhysChem 2011, 12, 3337–3347. (83) Kaukonen, M.; Söderhjelm, P. Proton Transfer at Metal Sites in Proteins Studied by Quantum Mechanical Free-Energy Perturbations. J. Chem. Theory Comput. 2008, 4, 985–1001. (84) Rod, T. H.; Ryde, U. Accurate QM/MM Free Energy Calculations of Enzyme Reactions: Methylation by Catechol O-Methyltransferase. J. Chem. Theory Comput. 2005, 1, 1240–1251. (85) Luzhkov, V. B.; Warshel, A. Microscopic Models for Quantum Mechanical Calculations of Chemical Processes in Solutions: LD/AMPAC and SCAAS/AMPAC Calculations of Solvation Energies. J. Comput. Chem. 1992, 13, 199–213. (86) Li, J.; Farrokhnia, M.; Rulíšek, L.; Ryde, U. Catalytic Cycle of Multicopper Oxidases Studied by Combined Quantum- and Molecular-Mechanical Free-Energy Perturbation Methods. J. Phys. Chem. B 2015, 119, 8268–8284. (87) Ullmann, G. M.; Knapp, E. W. Electrostatic Models for Computing Protonation and Redox Equilibria in Proteins. Eur. Biophys. J. 1999, 28, 533–551. (88) Åqvist, J.; Warshel, A. Computer Simulation of the Initial Proton Transfer Step in Human Carbonic Anhydrase I. J. Mol. Biol. 1992, 224, 7–14. (89) Olsson, M. H. M.; Hong, G.; Warshel, A. Frozen Density Functional Free Energy Simulations of Redox Proteins: Computational Studies of the Reduction Potential of Plastocyanin and Rusticyanin. J. Am. Chem. Soc. 2003, 125, 5025–5039. (90) Schutz, C. N.; Warshel, A. What Are the Dielectric “constants” of Proteins and How to Validate Electrostatic Models? Proteins Struct. Funct. Genet. 2001, 44, 400–417. (91) Boresch, S. The Role of Bonded Energy Terms in Free Energy Simulations - Insights 27 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(92)

(93) (94) (95) (96) (97)

(98)

(99) (100) (101) (102)

(103) (104)

(105) (106)

(107)

(108)

(109)

from Analytical Results. Mol. Simul. 2002, 28, 13–37. Riccardi, D.; Schaefer, P.; Cui, Q. pKa Calculations in Solution and Proteins with QM/MM Free Energy Perturbation Simulations : A Quantitative Test of QM / MM Protocols. J. Phys. Chem. B 2005, 109, 17715–17733. Böttcher, C. J. F. Theory of Electric Polarization; Elsevier: Amsterdam, 1973. Ryde, U.; Nilsson, K. Quantum Chemistry Can Locally Improve Protein Crystal Structures. J. Am. Chem. Soc. 2003, 125, 14232–14233. Kleywegt, G. J.; Jones, T. A. Databases in Protein Crystallography. Acta Crystallogr. Sect. D Biol. Crystallogr. 1998, 54, 1119–1131. Pannu, N. S.; Read, R. J. Improved Structure Refinement Through Maximum Likelihood. Acta Crystallogr. Sect. A Found. Crystallogr. 1996, A52, 659–668. Adams, P. D.; Pannu, N. S.; Read, R. J.; Brünger, A. T. Cross-Validated Maximum Likelihood Enhances Crystallographic Simulated Annealing Refinement. Proc. Natl. Acad. Sci. U. S. A. 1997, 94, 5018–5023. Brünger, A. T.; Adams, P. D.; Clore, G. M.; Delano, W. L.; Gros, P.; GrosseKunstleve, R. W.; Jiang, J.; Kuszewski, J.; Nilges, M.; Pannu, N. S.; et al. Crystallography & NMR System : A New Software Suite for Macromolecular Structure Determination. Acta Crystallogr. Sect. D Biol. Crystallogr. 1998, 54, 905– 921. Brunger, A. T. Version 1.2 of the Crystallography and NMR System. Nat. Protoc. 2007, 2 (11), 2728–2733. Engh, R. A.; Huber, R. Accurate Bond and Angle Parameters for X-Ray Protein Structure Refinement. Acta Crystallogr. Sect. A 1991, 47, 392–400. Kleywegt, G. J. Crystallographic Refinement of Ligand Complexes. Acta Crystallogr. Sect. D Biol. Crystallogr. 2006, 63, 94–100. Afonine, P. V; Ralf, W.; Headd, J. J.; Thomas, C. Towards Automated Crystallographic Structure Refinement with Phenix.refine. Acta Crystallogr. Sect. D Biol. Crystallogr. 2012, 68, 352–367. Tickle, I. J. Statistical Quality Indicators for Electron-Density Maps. Acta Crystallogr. Sect. D Biol. Crystallogr. 2012, 68, 454–467. Uranga, J.; Mikulskis, P.; Genheden, S.; Ryde, U. Can the Protonation State of Histidine Residues Be Determined from Molecular Dynamics Simulations? Comput. Theor. Chem. 2012, 1000, 75–84. Nilsson, K.; Ryde, U. Protonation Status of Metal-Bound Ligands Can Be Determined by Quantum Refinement. J. Inorg. Biochem. 2004, 98, 1539–1546. Nilsson, K.; Hersleth, H.-P.; Rod, T. H.; Andersson, K. K.; Ryde, U. The Protonation Status of Compound II in Myoglobin, Studied by a Combination of Experimental Data and Quantum Chemical Calculations: Quantum Refinement. Biophys. J. 2004, 87, 3437–3447. Warshel, A.; Dryga, A. Simulating Electrostatic Energies in Proteins: Perspectives and Some Recent Studies of pK As, Redox, and Other Crucial Functional Properties. Proteins Struct. Funct. Bioinforma. 2011, 79, 3469–3484. Alexov, E.; Mehler, E. L.; Baker, N.; M. Baptista, A.; Huang, Y.; Milletti, F.; Erik Nielsen, J.; Farrell, D.; Carstensen, T.; Olsson, M. H. M.; et al. Progress in the Prediction of pK a Values in Proteins. Proteins Struct. Funct. Bioinforma. 2011, 79, 3260–3275. Cheng, J.; Liu, X.; VandeVondele, J.; Sulpizi, M.; Sprik, M. Redox Potentials and Acidity Constants from Density Functional Theory Based Molecular Dynamics. Acc. Chem. Res. 2014, 47, 3522–3529.


Page 28 of 46

Page 29 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 1. Hydrogen bonds around the studied residues in the crystal structure and in the various MD simulations (the reference and the alternative protonation state for each of the investigated residues, A and C subunits; H-bonds with occupancies lower than 10% in all MD simulations were excluded; the full data is shown in Table S10). % is the percentage of snapshots in which the H-bond is found. A H-bond was considered to have a H–O/N distance of less than 2.5 Å and a H–O/N–X angle larger than 90º. A percentage larger 100% means that there are more than one such interaction (HOH) or that there is more than one hydrogen atom involved (e.g. Arg and Lys).

Atom1 His-195 HD1 HD1 ND1 HE2 His-274 ND1 HD1 NE2 HE2 His-362 ND1 HD1 HE2 HE2 HE2 HE2 His-451 HD1 HD1 HD1 NE2 HE2 HE2 Glu-153 OE1 OE1 OE2 HE2 OE2 OE1 Glu-380 OE1 OE2 OE1 OE2 OE2

Atom2 O Arg-277 O Wat H1 Wat S2B Wat Wat H Phe-299 Wat

Crystal structure A C

2.84 2.83 3.22 3.20

2.94 2.92 3.04 3.04

% Reference A C HIE

94 82 HID 3 96 83

84 78 3 98 69

Wat Wat NE2 His-274 O Tyr-297 NE2 His-451 Wat

2.83 2.83

HIE 75 74

2.88 2.88

4 68

4 78 11

OD2 Asp-228 O Asp-234 OG Ser-238 Wat NE2 His-274 Wat

2.85 2.88

8 HID 82 22 0 95

2.81 2.81

74 7 40 97

GLH HH Tyr-64 HG Ser-82 H Gly-84 O Pro-85 H Glu-153 Wat HE22 Gln-191 HE22 Gln-191 H Gly-356 H Gly-356 HE2 His-442

2

2.71 2.69 2.70 2.71

3.04 3.04

95

Alternative A C HIP 35 47 94 92 94 97 HIE 92 94

97 96 HIP 110 105 20 34 5 0 76 66 2 13 HIP 24 78 45 64 89 21 3 9 48 90 40 GLU 91 60 5 100 97

100

87 99 GLU 8 63 20 58 18

95 96 96 96 GLH 88 18


Page 30 of 46

Page 31 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


OE1

Wat

OE2 Wat HE2 Wat Glu-440 OE1 HZ Lys-428 OE2 HZ Lys-428 OE1 H His-442 OE2 H His-442 OE1 H Ser-443 OE2 H Ser-443 OE1 HG Ser-443 OE2 HG Ser-443 OE1 H Asp-445 OE2 HD2 Asp-445 HE2 OD1 Asp-445 OE1 HH Tyr-446 OE2 HH Tyr-446 HE2 OH Tyr-446 OE1 Wat OE2 Wat HE2 Wat Asp-445 OD1 HZ Lys-428 OD2 HZ Lys-428 OD2 HE/HH Arg-439 OD1 H Glu-440 OD2 H Glu-440 HD2 O Asp-445 OD1 Wat OD2 Wat HD2 Wat

2.58 2.58 3.02 3.01

226

227

10

185

203

6 4 98 99 GLU 1 13 3 27 54 44 0 30 96 50 1 86 84 44 85 1 24 73

GLH 17 0 20 12

2.82 2.81

0 9

40 14

2.87 2.87

26 21

99 4 5

2.98 3.00 2.76 2.77

2.83 2.79

2.83 2.83 2.77 2.78 2.88 2.87

21 71 137 29 25 50 2 ASP 1 67 1 24 26 19 61 24 99 99 212 111

90 106

71

10 11

3 1

149 116

59 111

ASH 6 99 50 0 61 17 102


58 66 20 22 98 18 10


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 2. Hydrogen bonds around the homocitrate ligand in the crystal structure and in the four MD simulations with varying protonation state of homocitrate (1Ha, 2H, 1Hc, and 0H). The entities in the table are the same as in Table 1 and full data is shown in Table S11. Atom 1 is always in homocitrate. Atoms

Crystal % structure 1Ha 2H 1Hc 0H 1 2 A C A C A C A C A C O1 HE21 Gln-191 2.89 2.86 83 O1 HH Arg-105B 18 1 1 26 37 9 72 O1 Wat 226 263 102 7 184 184 296 196 O2 HH Arg-105B 27 2 2 1 6 89 O2 Wat 2.74 2.75 243 285 16 92 73 101 317 229 2.83 2.84 HO2 Wat 91 1 59 72 O3 H Ile-425 2.83 2.84 38 14 10 30 34 2 61 O3 H Lys-426 2 28 O3 HZ Lys-426 13 79 44 1 10 31 12 90 O3 HG Ser-101B 2 17 29 1 24 60 0 O3 HE/H Arg-105B 44 11 25 146 48 10 53 0 O3 Wat 2.87 2.86 172 127 148 187 184 117 138 80 2.88 2.86 O4 H Ile-425 26 63 16 20 53 1 94 O4 HZ Lys-426 59 18 34 27 5 16 O4 HE2 Glu-440 37 O4 HG Ser-580 5 19 38 0 9 26 49 O4 HH Arg-105B 24 18 31 112 66 41 50 6 O4 Wat 2.84 2.85 135 163 165 210 173 175 155 123 2.82 2.81 2.72 2.72 O5 Wat 2.81 2.80 98 96 104 119 98 90 90 90 O6 Wat 2.83 2.84 177 211 191 186 194 213 203 209 2.73 2.73 O7 Wat 2 0 6 1 135 134 121 118 HO7 O1 homocitrate 15 59 HO7 Wat 72 97 94 86


Page 32 of 46

Page 33 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3. Calculated pKa values (in pKa units) and free-energy differences (∆G in kJ/mol) between the various protonation states of homocitrate from QM-cluster calculations of the isolated MoFe cluster. All geometries were optimised (opt) with the TPSS-D3/def2-SV(P) method, either in vacuum or in a COSMO continuum solvent with a dielectric constant (ε) of 4 or 80. Then, energies were calculated by single-point energy calculations with two DFT methods (TPSS-D3 or B3LYP-D3) and two basis set (SVP = def2-SV(P) or TZPD = def2TZVPD), in the COSMO solvent with ε = 4 or 80. ∆∆Gtd in the last column is the thermal corrections to the pKa values or free energies, obtained from frequency calculations in vacuum. It is added to all entries in the table. A pKa value larger than 7 indicates that the protonated state is predicted to be most stable at pH 7. The last row shows the most stable protonation state. Solvent opt Vacuum ε = 80 ε=4 energy ε = 80 ε = 80 ε=4 DFT, energy TPSS TPSS TPSS B3LYP B3LYP Basis set, energy SVP SVP TZPD TZPD TZPD pKa(2H/1Ha) 6.8 6.2 6.4 6.1 53.1 pKa(2H/1Hc) 10.6 11.6 11.6 11.7 58.8 pKa(1Ha/0H) 31.3 30.2 24.9 22.8 88.0 pKa(1Hc/0H) 27.5 24.8 19.7 17.1 82.3 ∆G(1Ha/1Hc) 21.9 30.9 29.9 32.5 32.5 Most stable 1Ha 1Ha 1Ha 1Ha 2H


Vacuum ∆∆Gtd TPSS SVP -5.3 -5.1 -7.0 -7.2 1.3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 4. Calculated pKa values (in pKa units) and free-energy differences (∆G in kJ/mol) between the various protonation states of homocitrate from QM/MM-2QM calculations of nitrogenase. The calculations were performed with two DFT methods and two basis set (SVP = def2-SV(P) or TZPD = def2-TZVPD). In the last four columns, solvation effects were treated with the QM/MM-PBSA approach with charges from either a vacuum (vac) or a polarised wavefunction (ptch). All values were corrected with the thermal corrections ∆∆Gtd from Table 3, with the term from the corresponding calculation for acetic acid subtracted. All pKa values were shifted with respect to the pKa value of the Glu-212A group (4.0). Therefore, a pKa value larger than 7 indicates that the protonated state is predicted to be most stable at pH 7. The last row shows the most stable protonation state. Method QM/MM QM/MM-PBSA-vac QM/MM-PBSA-ptch DFT TPSS TPSS B3LYP TPSS B3LYP TPSS B3LYP Basis SVP TZPD TZPD SVP TZPD SVP TZPD pKa(2H/1Ha) -35.7 -28.9 -30.2 5.1 6.7 10.5 12.1 pKa(2H/1Hc) -26.9 -26.2 -27.8 23.0 22.1 25.7 24.8 pKa(1Ha/0H) 32.1 28.5 27.0 61.6 56.5 65.1 60.0 pKa(1Hc/0H) 27.4 26.1 25.5 43.6 41.7 49.7 47.7 ∆G(1Hc/1Ha) 27.6 15.8 13.7 102.4 88.5 87.0 73.1 26.7 13.8 8.5 103.3 85.2 88.5 70.3 Most stable 1Ha 1Ha 1Ha 1Ha 1Ha 2H 2H


Page 34 of 46

Page 35 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 5. pKa values (in pKa units) and free-energy differences (∆G in kJ/mol) between the various protonation states of homocitrate in nitrogenase from the QTCP calculations. The raw results were obtained with the TPSS-D3/def2-SV(P) approach. They were then corrected with the thermal corrections ∆∆Gtd from Table 3 and acetic acid, as well with a correction to B3LYP-D3/def2-TZVPD level from Table 4. SE is the standard error obtained by error propagation from the standard deviations over the 200 snapshots in the individual FEP calculations. The pKa values were shifted with respect to the pKa value of Glu-212A (4.0). Therefore, a pKa value larger than 7 indicates that the protonated state is most stable at pH 7. The last row shows the most stable protonation state. Raw Corrected SE pKa(2H/1Ha) 0.1 2.5 1.3 pKa(2H/1Hc) 7.0 7.2 1.3 pKa(1Ha/0H) 24.9 19.1 1.7 pKa(1Hc/0H) 20.2 17.3 1.6 ∆G(1Hc/1Ha) 39.5 26.9 10.5 27.4 10.5 12.2 Most stable 1Ha 1Ha



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 6. Results of the quantum-refinement calculations for the four protonation states of homocitrate in the nitrogenase 3U7Q structure at the TPSS/def2-SV(P) level of theory. Three quality criteria are shown, Rfree, R, and RSZD. The corresponding results at the B3LYP/def2SV(P) level are given in Table S12.

2H 1Ha 1Hc 0H

Rfree 0.164 0.169 0.169 0.164

R RSZD 0.148 3.2 0.153 3.0 0.152 8.0 0.148 10.0


Page 36 of 46

Page 37 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. The (a) MoFe cluster, (b) P cluster, and (c) Ca site in nitrogenase, showing also the atom names8 and the models used to calculate charges for the MM force field.

a

b

c 37 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2. The four considered protonation states of homocitrate, 2H, 1Ha, 1Hc and 0H. Atom numbers are also shown. Non-polar H atoms are omitted


Page 38 of 46

Page 39 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. (a) The thermodynamic cycle employed in the QTCP calculations. (b) A schematic picture of the various perturbations performed for the proton-transfer reactions.86 The four letters represent the proton on the homocitrate ligand, charges and coordinates of both QM systems, and the proton on Glu-212A. The proton can be either modelled as a full proton, with zeroed dummy parameters (D), or be deleted (∅). The charges and coordinates can be either those of the starting state (A, in which the homocitrate is protonated and Glu-212A is deprotonated) or the final state (B, in which the homocitrate is deprotonated and Glu-212A is protonated). The number of perturbation steps are also indicated for each perturbation.

a

b



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4. The local surroundings of the considered residues in the crystal structure of nitrogenase: (a) His-195, (b) His-274, (c) His-362, (d) His-451, (e) Glu-153, (f) Glu-380, (g) Glu-440 and Asp-445, as well as (h) homocitrate. The central residue is shown in thick lines, residues mentioned in the text and in the tables in intermediate lines, other residues within 6 Å in thin lines (not labelled) and water molecules within 6 Å as red balls.

a

b 40 ACS Paragon Plus Environment

Page 40 of 46

Page 41 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


c

d



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

e

f 42 ACS Paragon Plus Environment

Page 42 of 46

Page 43 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


g

h 43 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5. RMSD for the various simulations: (a) involving His-192 (HIE state = reference or HIP = alternative state), 271 (HID or HIE), 359 (HIE or HIP) and 448 (HID or HIP), Glu-150 (reference protonated), 377 (reference deprotonated), 437 (reference protonated and Asp-445 deprotonated), or (b) homocitrate in the four protonation states 1Ha, 2H, 1Hc and 0H. Two RMSD values are given for each simulation, one for the studied residue alone and one for all residues with at least one atom within 3.5 Å of the studied residue. The RMSD values are the average over the two halves of the protein. The standard error of the RMSD values are 0.003– 0.015 Å, typically larger for the isolated residue and larger when the 3.5-Å surrounding is included (largest for subunit C of the simulation with the alternative protonation of His-362).

a

b


Page 44 of 46

Page 45 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. Electron-density maps of the (a) 1Ha and (b) 2H protonation states of the homocitrate ligand in the quantum refinement of the nitrogenase 3U7Q structure performed at the TPSS/def2-SV(P) level of theory. The 2mFo –DFc maps are contoured at 1.0 σ and the mFo – DFc maps are contoured at +3.0 σ (green) and –3.0 σ (red). The corresponding results at the B3LYP/def2-SV(P) level are given in Figure S1.

a

b 45 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC Graphic


Page 46 of 46

Protonation States of Homocitrate and Nearby Residues in

Recommend Documents