Putative reaction mechanism of nitrogenase after dissociation of a sulfide ligand

We have investigated the implications of the recent crystallographic ﬁndings that the m 2 -bridging S2B sulﬁde ligand may reversibly dissociate from the active-site FeMo cluster of nitrogenase. We show with combined quantum mechanical and molecular mechanical (QM/MM) calculations that once S2B has dissociated, N 2 may bind in that position and can be protonated to two NH 3 groups by thermodynamically favourable steps. The substrate forms hydrogen bonds with two protein ligands, Gln-191 and His-195. For all steps, we have studied three possible protonation states of His-195 (protonated on either ND1, NE2 or both). We ﬁnd that the thermodynamically favoured path involves an end-on NNH 2 structure, a mixed side-on/end-on H 2 NNH structure, a side-on H 2 NNH 2 structure, a bridging NH 2 structure and a bridging NH 3 structure. In all cases, His-195 seems to be protonated on the NE2 atom. Dissociation of the NH 3 product is often unfavourable and requires either further reduction or protonation of the cluster or rebinding of S2B. In conclusion, our calculations show that dissociation of S2B gives rise to a natural binding and reaction site for nitrogenase, between the Fe2 and Fe6 atoms, which can support an alternating reaction mechanism with favourable energetics. (cid:1) 2020 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Nitrogenase (EC 1.18/19.6.1) is the only enzyme that can cleave the strong triple bond in N 2 [1][2][3].This is one of the most important reaction in nature: Even if the atmosphere of earth contains 78% N 2 , nitrogen is typically a limiting element for plant growth and a major component of artificial fertilizers.The reason for this is that the triple bond in N 2 is very strong, making the molecule inert [3,4].
The nitrogenase reaction is demanding, requiring 16 molecules of ATP, eight electrons and eight protons to cleave one NN triple bond: [ Nitrogenase is a large and complicated system [1][2][3][5][6][7][8][9].It consists of two components.The Fe protein contains a Fe 4 S 4 cluster and binds two molecules of ATP.In the ATP-loaded state, it docks to the MoFe protein and delivers electrons to the active-site FeMo cluster, via a Fe 8 S 7 Cys 6 site, called the P-cluster.After hydrolysis of the ATP molecules, the Fe protein is released and a new reduced and ATP-loaded Fe-protein may bind.The active site consists of a MoFe 7 S 9 C(homocitrate)CysHis cluster, [5][6][7][8][9] although alternative nitrogenases exist, in which the Mo ion is replaced by V or Fe [10].
Nitrogenase has been thoroughly studied with biochemical, kinetic, crystallographic and spectroscopic methods [1][2][3][11][12][13][14][15][16].The reaction is often described by nine intermediates, E 0 -E 8 , differing in the number of added electrons and protons, the Lowe-Thorneley cycle [17].It is currently believed that four electrons have to be added before N 2 can bind and that H 2 formation through reductive elimination is a prerequisite for the binding [1].It is also believed that Fe2 and Fe6 are involved in the binding of N 2 (the numbering of the Fe ions is shown in Fig. 1a) [3].
Many density-functional theory (DFT) studies have been performed to study the reaction mechanism of nitrogenase, but no consensus has been reached so far [3,12,[18][19][20][21][22][23][24][25][26][27][28][29][30][31][32].For example, it has been much disputed whether N 2 is protonated first by three protons on one of the N atoms, which leaves as NH 3 at the E 5 stage, before the second N atom is protonated, or if the protons are added alternatively to both N atoms, so that HNNH and H 2 NNH 2 are intermediates and the first NH 3 molecule does not dissociate until the E 7 stage [3,33].These two mechanisms are called the distal and alternating pathways, respectively, and they are illustrated in Fig. 2. The distal pathway was originally suggested by Chatt and has gained support from inorganic model complexes [34][35][36][37][38].
However, it has been suggested to apply also for nitrogenase by several authors [20,39].The alternating pathway is supported by the fact that nitrogenase can use hydrazine as a substrate and that hydrazine is released upon acid or base hydrolysis of the enzyme during turnover [1,3,40,41].Moreover, it has been shown that N 2 , N 2 H 2 , CH 3 N 2 H and N 2 H 4 all react via common intermediate, which strongly supports an alternative mechanism [3,33].
Most studies have assumed that the FeMo cluster remains intact throughout the reaction.However, in 2014, Rees et al. published a crystal structure of CO-inhibited Mo-nitrogenase, in which CO replaces one of the sulfide ions (S2B, bridging the Fe2 and Fe6 ions; cf.Fig. 1a) of the FeMo cluster and binds to the Fe2 and Fe6 atoms [8].They suggested that something similar might happen during the normal reaction mechanism of the enzyme, i.e. that the S2B ligand is labile and may dissociate reversibly during the reaction, forming a binding site for the N 2 substrate.This was recently supported by another crystal structure of V-nitrogenase, in which S2B had also dissociated from the cluster [39].In that structure, the nearby Gln-176 had changed its position, forming a pocket close to the active site, in which they found a prominent density, assigned to a HS -ion.This indicates that S2B may be protonated and displaced from the active site, but kept stored nearby inside the protein so that it can easily bind back again.Moreover, they found a light atom bound to the FeV cluster at the former position of S2B.They interpreted it as a nitrogen atom, which could represent the E 6 or E 7 intermediate (formally NH 2-or NH 2 -) of a dis-tal reaction mechanism.Yet, Bjornsson and coworkers showed by a combined quantum mechanics and molecular mechanics (QM/ MM) study that the crystal structure most likely contains a OH - ligand, rather than a N 2 -derived group, [42] although this study has not convinced the crystallographers [39].However, even if the crystal structure involves OH -, it is still possible that S2B reversibly dissociates from the cluster during the reaction mechanism and N 2 binds to that site.In fact, Nørskov et al. suggested such a mechanism for Mo-nitrogenase in 2015, based on QM-cluster calculations [20].In their mechanism, N 2 binds end-on to both Fe2 and Fe6 and is then sequentially protonated along a distal pathway.However, these calculations involved a oversimplified model of the cluster (with histidine modelled by NH 3 and homocitrate by two OH -ions), no parts of the surrounding enzyme were included and it was assumed that the E 0 state is doubly protonated, which recent quantum-refinement calculations have shown not to be the case [43].
Therefore, we here provide a more detailed study of the reaction mechanism of Mo-nitrogenase, assuming that the S2B group dissociates from the FeMo cluster and N 2 binds in its position.Our calculations employ QM/MM methods, including the entire solvated heterotetrameric structure of nitrogenase in the calculations.We thoroughly investigate the energetics and electronic structure of possible intermediates of both the alternating and distal mechanism and also different protonation states of His-195, which can form hydrogen bonds to the N 2 -derived ligand.Even if there is no conclusive evidence that dissociation of S2B is involved in the Fig. 1. a) The Mo-nitrogenase FeMo cluster with S2B dissociated (it was bridging Fe2 and Fe6 before dissociation) and atom names from the 3U7Q crystal structure [7].b) The QM system used in the QM/MM calculations (with H 2 NNH bound).Fig. 2. The alternating and distal pathways suggested for the second part (E 4 -E 8 ) of the nitrogenase mechanism, [3,33] showing the N 2 -derived intermediates.catalytic mechanism of nitrogenase, this study provides structures and energies for such a mechanism that can be compared to experiments and to similar studies of the reaction without dissociation of S2B.

The protein
The QM/MM calculations were based on the 1.0-Å crystal structure of Mo-nitrogenase from Azotobacter vinelandii (PDB code 3U7Q) [7].The setup of the protein is identical to that of our previous studies of the protein [44][45][46][47][48].The entire heterotetramer was included in the calculations, because the various subunits are entangled without any natural way to separate them.The QM calculations were concentrated on the FeMo clusters in the C subunit because there is a buried imidazole molecule from the solvent rather close to the active site (~11 Å) in the A subunit.The two P clusters and the FeMo cluster in subunit A were modelled by MM in the fully reduced and resting states, respectively [44].
The protonation states of all residues were the same as before [44]: All Arg, Lys, Asp and Glu residues were assumed to be charged, except Glu-153, 440 and 231D (a letter ''D" after the residue number indicates that it belongs to that subunit; if no letter is given, it belongs to subunit C; subunits A and B are identical to the C and D residues).Cys residues coordinating to Fe ions were assumed to be deprotonated.His-274, 451, 297D, 359D and 519D were assumed to be protonated on the ND1 atom, 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.The homocitrate ligand was modelled in the singly protonated state with a proton shared between the hydroxyl group (which coordinates to Mo) and the O1 carboxylate atom.This protonation state was found to be the most stable one in an extensive QM/MM, molecular dynamics and quantumrefinement study [44] and this protonation state is also supported by another recent study [49].
The protein was solvated in a sphere with a radius of 65 Å around the geometrical centre of the protein.160 Cl -and 182 Na + ions were added at random positions (but not inside the protein [44]) to neutralise the protein and give an ionic strength of 0.2 M [50].The final system contained 133 915 atoms.The added protons, counter ions and water molecules were optimised by a simulated annealing calculation (up to 370 K), followed by a minimisation, keeping the other atoms fixed at the crystal-structure positions [44].
All MM calculations were performed with the Amber software [51].For the protein, we used the Amber ff14SB force field [52] and water molecules were described by the TIP3P model [53].For the metal sites, the MM parameters were the same as in our previous investigation [44].The metal sites [44,48] were treated by a non-bonded model [54] and charges were obtained with the restrained electrostatic potential method, obtained at the TPSS/ def2-SV(P) level of theory [55,56] and sampled with the Merz-Kollman scheme [57].
The FeMo cluster was modelled by MoFe 7 S 8 C(homocitrate) (CH 3 S)(imidazole), where the two last groups are models of Cys-275 and His-442.In addition, all groups that form hydrogen bonds to the FeMo cluster were also included in the QM system, viz.Arg-96 and His-195 (sidechains), Ser-278 and Arg-359 (both backbone and sidechain, including the Ca and C and O atoms from Arg-277), Gly-356, Gly-357 and Leu-358 (backbone, including the Ca and C and O atoms from Ile-355), as well as two water molecules.Following the crystal structure of V-nitrogenase with the putative N-intermediate, [39] S2B was removed from the FeMo cluster.Moreover, the sidechain of Gln-191 was included in the QM system and this group was flipped so that OE1 points towards the former position of S2B and His-195, whereas the NE2 atom of Gln-191 still forms a hydrogen bond to one of the carboxylate groups of the homocitrate ligand (1.8-1.9Å H . . .O distance).Finally, Phe-381 was also included, because it is close to the putative N 2 binding site.The QM system involved 170-178 atoms in total (depending on the number of added protons and N atoms) and is shown in Fig. 1b.The charge was -1, except in a few cases when electrons and protons were not added together.
Experiments have shown that the ground spin state of E 0 and E 2 are quartets with a surplus of three a electrons, whereas E 4 is a doublet [3,12].Consequently, we used these spin states for these three states.For the other oxidation states, we checked which of the two or three lowest spin states has the most favourable energy at the TPSS and B3LYP/def2-SV(P) levels of theory.
The electronic structure of all QM calculations was obtained with the broken-symmetry (BS) approach [30]: Each of the seven Fe ions were modelled in the high-spin state, with either a surplus of a (four Fe ions) or b (three Fe ions) spin.Such a state can be selected in 35 different ways ( 7! 3!4! ) [45].A starting wavefunction was obtained by first optimising the all-high-spin state with 35 unpaired electrons and then changing the total a and b occupation numbers to the desired net spin.This gave one of the BS states.The other BS states were obtained by simply swapping the coordinates of the Fe ions [66].In some cases, we instead used the fragment approach by Szilagyi and Winslow to obtain a proper BS state [67].The various BS states are named by listing the number in the Noodleman nomenclature (BS1-10) [30], followed by the numbers of the three Fe ions with minority spin, e.g.BS7-235, indicating that Fe2, Fe3 and Fe5 have b spin.For ligand-bound structures, previous studies have shown that the spin on some Fe ions is often quite low and sometimes states with only two Fe b spin might become competitive [46].Such states were also considered and they are named by just giving the two Fe ions with negative spin, e.g.BS-14, indicating that Fe1 and Fe4 have b spin.We have thoroughly studied the 35 BS states for the resting state, the protonated resting state and the reduced state, and how their energies vary with the QM method, the size of the basis set, the geometry and the influence of the surroundings [45].For each E n state, we studied the relative energy of all BS states with both TPSS and B3LYP/ def2-SV(P).
As have been discussed before [45,47], TPSS-D3/def2-SV(P) calculations give geometries that reproduce the crystal structure of the resting state of nitrogenase excellently with average and maximum deviations of 0.05 and 0.09 Å for the metal-metal and 0.02 and 0.06 Å metal-ligand distances, respectively and a root-meansquared-deviation (RMSD) of 0.06 Å for the metals and the firstsphere ligands.This is similar to the results obtained with the TPSSh approach [49] and appreciably better than with the B3LYP-D3/def2-SV(P) method, which gives average and maximum deviations of 0.08 and 0.12 Å for the metal-metal and 0.04 and 0.11 Å metal-ligand distances, respectively and a RMSD of 0.08 Å.Therefore, we discuss primarily the TPSS-D3/def2-SV(P) results.

QM/MM calculations
The QM/MM calculations were performed with the COMQUM software [68,69].In this approach, the protein and solvent are split into two subsystems: System 1 (the QM region) was relaxed by QM methods.System 2 contained the remaining part of the protein and the solvent and it was kept fixed at the original coordinates (equilibrated crystal structure to avoid the risk that different calculations end up in different local minima).The total system was spherical and non-periodic with 133 915 atoms.
In the QM calculations, system 1 was represented by a wavefunction, whereas all the other atoms were represented by an array of partial point charges, one for each atom, taken from the MM setup.Thereby, the polarisation of the QM system by the surroundings is included in a self-consistent manner (electrostatic embedding).When there is a bond between systems 1 and 2 (a junction), the hydrogen link-atom 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 [68,70].All atoms were included in the point-charge model, except the CL atoms [71].
The total QM/MM energy in ComQum was calculated as [68,69] where E HL QM1þptch2 is the QM energy of the QM system truncated by HL atoms and embedded in the set of point charges modelling system 2 (but excluding the self-energy of the point charges).E HL MM1;q 1 ¼0 is the MM energy of the QM system, still truncated by HL atoms, but without any electrostatic interactions.Finally, E CL MM12;q1¼0 is the classical energy of all atoms in the system with CL atoms and with the charges of the QM region set to zero (to avoid double-counting of the electrostatic interactions).Thus, ComQum employs a subtractive scheme with electrostatic embedding and van der Waals linkatom corrections [72].No cutoff is used for any of the interactions in the three energy terms in Eqn. 3.
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.

Result and discussion
Inspired by the recent crystal structures of Mo and Vnitrogenase, showing that the S2B sulfide group may dissociate reversibly from the active site [8,39], we investigate whether N 2 may bind to the empty position between Fe2 and Fe6 ions, previously occupied by S2B, and react to form two molecules of NH 3 .Since the great majority of previous QM studies of nitrogenase have been performed on the Mo enzyme [3,12,[18][19][20][21][22][23][24][26][27][28][29][30][31][32], we decided to study that enzyme.Therefore, we removed the S2B ligand from the FeMo cluster and replaced it with N 2 in the bridging position between Fe2 and Fe6.Moreover, we flipped the Gln-191 residue in the same way as observed for Gln-176 in the crystal structure of V-nitrogenase, i.e. so that the OE1 atom points towards N 2 and His-195 (corresponding to His-180 in V-nitrogenase), whereas the NE2 atom still forms a hydrogen bond to the acetate group of the homocitrate ligand (cf.Fig. 1b).
Dissociation of S2B opens a natural binding site for N 2 , which strongly facilitates the computational modelling of the reaction mechanism (in a previous study, we investigated ~60 different binding modes, assuming that S2B remains bound [48]).Moreover, it avoids the problem of defining the E 4 state before H 2 dissociation and N 2 binding, for which there are even more possibilities and for which there has been strong disagreement between various computational studies [14,21,43,47,[73][74][75][76].On the other hand, the protonation of His-195, which is within hydrogen-bonding distance to the substrate, can be expected to strongly affect the results.Therefore, we investigated the three possible protonation states of this residue, viz.singly protonated on either the ND1 or NE2 atoms, or protonated on both atoms.These states are called HID, HIE and HIP in the following.
As mentioned in the introduction, the nitrogenase mechanism is normally described by the nine E 0 -E 8 states in the Lowe-Thorneley cycle [17].E 0 is the resting state, which is in the Mo III  We then add successively four more electrons and protons, and investigate where these preferably go, when the N-N bond is cleaved and when NH 3 dissociates (taking three electrons, three protons and one nitrogen with it).
We try to determine the most stable conformation of the E 4 -E 8 states (i.e.protonation and binding mode of the N 2 -derived intermediate).Thereby, we get information whether a reaction mechanism of Mo-nitrogenase with a dissociated S2B group is thermodynamically feasible and whether it follows a distal or alternating mechanism (cf.Fig. 2).

N 2 -bound E 4 structures, e 2 h 2 n 2
We start with the N 2 -bond E 4 state of the protein.To avoid speculating about the structure of the E 4 or E 2 states of the FeMo cluster, for which there are many possible protonation and electronic states, and for which there is still no consensus about the structure, [3,13,14,21,32,43,46,[74][75][76]78,79] because different DFT functionals give widely different results [47], we first assume that the two added protons have bound to N 2 , forming either HNNH or NNH 2 [48].
Thus, we started to optimise with QM/MM the structure of cis-HNNH bound side-on to Fe2 and Fe6 (with His-195 in the HIE state).The structure is shown in Fig. 3 (under HNNH).The two Fe-N distances are quite similar, 1.80 and 1.88 Å, as can be seen in Table 1.The N-N bond length is 1.29 Å, which is slightly longer than in free cis-HNNH (1.24 Å), optimised at the same QM level of theory.The structure is stabilised by a short hydrogen bond from HNNH to the carbonyl oxygen of Gln-191 (the H . . .O distance is 1.49 Å).The other proton of HNNH points against the NE2 atom of His-195 (2.18 Å), but this atom is already protonated so the interaction is probably not very favourable.The BS investigation of this structure, presented in Table S2, showed that it is most stable in the BS8-345 state.It has a low spin population on Fe6 (0.5), whereas that on Fe2 (2.7) is similar to that on the other Fe ions (2.3-2.7, in absolute terms), except that on Fe1 (3.3).The Mo ion has a low spin population (-0.1).
The corresponding HID structure has similar N-N and Fe-N bond lengths (1.28, 1.90 and 1.83 Å) and a similar strong hydrogen bond to Gln-191 (1.46 Å), but the hydrogen bond to His-195 is appreciably stronger (1.73 Å) because the NE2 atom is now not protonated.However, this HID structure is 35 kJ/mol less stable than the corresponding HIE structure (but only 3 kJ/mol with B3LYP), reflecting that the HIE structure is intrinsically more stable (e.g. by 41 kJ/mol for the E 0 state; 39 kJ/mol with B3LYP).
Fig. 3 also shows the structure of NNH 2 bound end-on to nitrogenase with His-195 in the HIE state.The unprotonated N atom interacts symmetrically with both Fe2 and Fe6 (the Fe-N bonds are 1.79 and 1.80 Å).The N-N bond is 1.27 Å, i.e. slightly longer than for free NNH 2 (1.21 Å).This structure is also stabilised by a strong hydrogen bond to the carbonyl oxygen of Gln-191 (1.48 Å), whereas the other proton points towards His-195 (2.19 Å), without making any favourable interaction.It turned out to be most stable in the BS10-147 state (Table S2; BS10-135 with B3LYP).It has low and negative spin populations on both Fe2 and Fe6 (-1.5 and -1.6), whereas the other Fe ions have larger spin populations (3.1 on Fe1 and 2.6-2.9 on the other ions).The spin on Mo is small but positive (0.2). Interestingly, this structure turned out to be 35 kJ/mol more stable than the side-on HNNH structure (31 kJ/mol with the best B3LYP states), although isolated NNH 2 is 64 kJ/mol less stable than cis-HNNH at the same level of theory (61-73 kJ/mol with different DFT functionals, basis sets and with or without a water-like continuum solvent).This agrees with the structure suggested by Nørskov and coworkers [20].The corresponding HID structure is similar (Table 1), but has a hydrogen bond between NE2 of His-195 and the second hydrogen on NNH 2 (1.64 Å).It is 37 kJ/mol less stable than the HIE structure.
To put these structures in a perspective, we also performed some investigations of the corresponding N 2 -bound structures (although far from exhaustive regarding possible positions of the two protons).An example with the two protons on S3B and His-195 (i.e. a HIP structure) is shown in Fig. 3, showing that N 2 is asymmetrically bound to Fe2 and Fe6.One N atom is 1.80 Å from Fe6 and 2.02 Å from Fe2, whereas the other is 2.10 Å from Fe2.The N-N distance is 1.17 Å, i.e. 0.06 Å longer than for free N 2 optimised at the same level of theory.N 2 is not stabilised by any hydrogen bonds.Instead, the doubly protonated His-195 forms a hydrogen bond with Gln-191 (2.01 Å).In this structure, Fe6 has a low spin population (2.0), whereas that on Fe2 (2.7) is similar to that of the other Fe ions (2.7-3.0, but 3.2 on Fe1).This structure is 212 kJ/mol less stable than the NNH 2 structure (117 kJ/mol with B3LYP), showing that protonation of the substrate is strongly favourable.
A structure with His-195 in the HIE state and S2A protonated instead is slightly less stable (13 kJ/mol), whereas the corresponding HID structure is strongly unfavourable (by 116 kJ/mol).However, a more stable N 2 structure was unexpectedly found if the OE1 atom of Gln-191 was protonated.As can be seen in Fig. 3, it is stabilised by hydrogen bonds from both Gln-191 and His-195 to the distal atom of N 2 .It has Fe-N bonds of 1.78 (Fe6) and 2.00 (Fe2) Å and a slightly longer N-N bond, 1.18 Å.It is actually 81 kJ/mol more stable than the structure protonated on His-195, but still 131 kJ/mol less stable than the NNH 2 structure (49 kJ/mol with B3LYP).
We also studied some structures with NNH bound to the cluster.An example with S3B protonated and His-195 in the HIE state is shown in Fig. 3. NNH binds side-on to the two Fe ions with Fe-N bonds lengths of 1.84 and 1.88 Å and a N-N bond of 1.24 Å.The protonated N atom (coordinated to Fe6) forms a strong hydrogen bond to OE1 of Gln-191 (1.30 Å), whereas the other N atom receives a hydrogen bond from His-195 (1.93 Å).This structure is 17 kJ/mol less stable than the best N 2 structure (with Gln-191 protonated; 36 kJ/mol with B3LYP) and 148 kJ/mol less stable than the best N 2 H 2 structure (85 kJ/mol with B3LYP).In the corresponding HID structure, the proton on NNH is spontaneously transferred to Gln-191, giving rise to the N 2 -S3B-Gln structure.In conclusion, our best N 2 -bound E 4 state involves end-on NNH 2 with HIE (cf.Fig. 3) and strongly favourable protonation of N 2 .

E 5 structures with N 2 H 3
Next, we added an electron and a proton and studied the various e 3 h 3 n 2 structures.Adding a proton to the side-on HNNH structure, gives the side-on HNNH 2 structure, shown in Fig. 4 (bottom left; HIE state).It has a somewhat longer Fe2-N bond for the doubly protonated N atom than the HNNH structure, 1.98 Å, whereas the Fe6-N bond is slightly shorter, 1.78 Å.The N-N bond length is 1.41 Å.There is still a hydrogen bond between the ligand and the carbonyl of Gln-191 (1.58 Å).We studied this complex in both the triplet and the quintet states and the two states were essentially degenerate (within 2 kJ/mol with both basis sets).For the triplet, BS7-247 was lowest in energy, whereas for the quintet, BS-17 Energies (TPSS-D3/def2-SV(P) and B3LYP-D3/def2-SV(P) energies in kJ/mol) and geometries (Å) of the various structures.The N-N and N-Fe bond lengths are given, together with possible hydrogen bond lengths to OE1 of Gln-191 (NH . . .OE1) and to NE2 or HE2 of His-195 (NH . . .NE2 or N . . .HE2).The last column (Q-H) gives the hydrogen-bond distance between OE1 of   (only two Fe ions with negative spin) was found to be lowest in energy (but BS7-247 was only 4 kJ/mol less stable).The spin densities are similar to those of the HNNH complex and there is no spin density on the HNNH 2 ligand.A slightly (6 kJ/mol) more stable structure can be obtained with the opposite protonation of His-195 (HID structure, also shown in Fig. 4).It involves a mixture of end-on and side-on binding with the singly protonated N atom bridging Fe2 and Fe6 with Fe-N distances of 1.89 and 1.88 Å, whereas the doubly protonated N atom is 2.00 Å from Fe2.The two protons of the NH 2 group form hydrogen bonds to Gln-191 and His-195, 1.70 and 1.68 Å, respectively.
However, the most stable structure is obtained by protonating the other N atom of HNNH (the one closes to .This gives the structure in the centre of Fig. 4 4).However, in contrast to the situation for the E 4 intermediates, the end-on NNH 3 structure turned out to be less stable than the best H 2 NNH structure by 51-83 kJ/mol (34-78 kJ/mol with B3LYP).This is in contrast to the study by Nørskov et al. [20], which suggested a NNH 3 intermediate.
Interestingly, the dissociation of NH 3 from the NNH 3 intermediate is not automatic.The energy for the dissociation reaction (assuming that NH 3 ends up in the water solvent and that the N atom remains bound to the cluster in the e 0 n 1 state) is highly unfavourable in the HID state (170 kJ/mol).In the HIE state, it is still unfavourable by 47 kJ/mol, but this is close to the expected gain in translational and rotational entropy from the released NH 3 .However, if the cluster is protonated (on His-195, without any addition of an electron, i.e. an e 3 h 4 n 2 state), the dissociation becomes favourable by 48 kJ/mol and if an electron is also added (e 4 h 4 n 2 , i.e. a E 6 state), the reaction is even more favourable, 78 kJ/mol.
The N 2 H 3 structures are appreciably more stable than the corresponding N 2 H 2 structures with a protonated His-195: The HIP-HNNH structure is 137 kJ/mol less stable than the HIE-H 2 NNH structure and the HIP-NNH 2 structure automatically reorganises to the HID-NNH 3 structure.This shows that protonation of the N 2 H 2 ligand is strongly favourable.
In conclusion, our calculations indicate that the most stable E 5 state is HIE-H 2 NNH, i.e. indicating an alternating reaction mechanism, in contrast to the suggestions by Nørskov and coworkers [20].Moreover, we show that the preferred binding mode is a mixture between end-on and side-on binding (cf.Fig. 4).

E 6 (e 4 h 4 n 2 ) structures with N 2 H 4
As already mentioned, we can obtain the HIP-NNH 3 structure, from which NH 3 can dissociate with a reaction energy of -78 kJ/mol.In fact, we can also find an end-on HNNH 3 complex with an intact N-N bond of 1.45 Å and two Fe-N bonds of 1.91 Å (HIE state).It has a very short hydrogen bond to Gln-191 (1.40 Å).It is most stable in the singlet BS10-147 state.The corresponding HID structure is similar (Fig. 5, right side; Fe-N = 1.94 and 1.92 Å; H . . .OE1 = 1.45 Å).In fact, it is 33 kJ/mol more stable than the HIE structure (36 kJ/mol with B3LYP), because it is stabilised by the hydrogen bond to His-195 (1.64 Å).From the HIE-HNNH 3 structure, the reaction energy for the dissociation of NH 3 is almost thermoneutral (+1 kJ/mol), so that the reaction is expected to be favourable if the reaction entropy is included.However, from the HID complex, the dissociation is strongly unfavourable (by 93 kJ/mol).
Moreover, we can also obtain a structure with a cleaved N-N bond, one N bound to the cluster and NH 4 + hydrogen bonded to it (1.54Å).The latter group forms also hydrogen bonds to Gln-191 and His-195 (both 1.62 Å; HID structure).It was 64 kJ/mol less stable than the HNNH 3 structure.However, all these structures are less stable than structures in which H 2 NNH 2 (i.e.hydrazine) binds side-on to the cluster.With HIE, it has a N-N bond length of 1.45 Å and the Fe-N bonds are 1.96 and 1.98 Å (Fig. 5, left side).It has a strong hydrogen bond between H 2 NNH 2 and Gln-191 (1.54 Å), but His-195 also forms a hydrogen bond to the OE1 atom of Gln-191 (2.24 Å).The structure is most stable in the BS10-135 state (but several other BS states are close in energy).It has low spin populations on both the Fe2 and Fe6 ions (1.7 and 1.4).This structure is 13-46 kJ/mol more stable than the HNNH 3 structures (78-114 kJ/mol with B3LYP).
The corresponding HID-H 2 NNH 2 structure has a similar geometry, with Fe-N bonds of 1.97 and 1.98 Å. H 2 NNH 2 forms hydrogen bonds to both Gln-191 (1.50 Å) and His-195 (1.72 Å).Still, it is 35 (83) kJ/mol less stable than the HIE structure, reflecting that intrinsically HIE is more stable than HID.The corresponding HIP-H 2 NNH and HIP-HNNH 2 structures are 139 and 180 kJ/mol less stable than the best HIE-H 2 NNH 2 structure (189 and 233 kJ/mol with B3LYP), again showing that protonation of the N 2 H 3 ligand is favourable.
In conclusion, our results indicate that the HIE-H 2 NNH 2 hydrazine complex is the most stable E 6 state, still supporting an alternating mechanism.

E 7 structures with N 2 H 5
Next, we study the E 7 states.When adding an extra electron and proton to the H 2 NNH 2 structure, a H 2 NNH 3 structure can be formed.It has an intact N-N bond of 1.46 Å, but the triply protonated N atom has dissociated from Fe2, whereas the other N atom still binds to Fe2 (Fe-N = 1.99 Å; left side of Fig. 6; HID state).It is stabilised by hydrogen bonds both to Gln-191 (1.63 Å) and His-195 (1.56 Å).
However, it is strongly favourable to cleave the N-N bond, forming a NH 2 + NH 3 complex (it can be done with an activation barrier of only 31 kJ/mol).In the HIE structure (Fig. 6, centre), the NH 2 group bridges Fe2 and Fe6 with Fe-N distances of 1.96 and 1.87 Å, whereas NH 3 binds terminally to Fe2 with a Fe-N distance of 2.05 Å. NH 2 forms a hydrogen bond to the OE1 atom of Gln-191 with a H . . .O distance of 1.71 Å.The HE2 atom of His-195 also forms a hydrogen bond to the OE1 atom with a distance of 2.15 Å.The singlet and triplet states were close in energy (4 kJ/mol).The former is most stable in the BS10-135 state, whereas the latter is more stable in the BS6-157 and BS7-235 states.
The corresponding HID complex is similar, with Fe-N bonds of 1.85, 1.99 and 2.07 Å.However, the hydrogen bond between NH 2 and the side-chain OE1 atom of Gln-191 is slightly shorter 1.66 Å.In addition, it has a hydrogen bond between NH 3 and the NE2 atom of His-195, but it is long, 2.09 Å, and with a suboptimal geometry.This structure is 101 kJ/mol less stable than the HIE structure with both TPSS and B3LYP.The corresponding HIP-H 2 NNH 2 structure is 196 kJ/mol less stable and therefore 297 kJ/mol less stable than the HIE-NH 2 + NH 3 structure, again illustrating the strongly favourable protonation of the ligand.We have also found a HIE structure with NH 3 on Fe6 and NH 2 bridging Fe2 and Fe6.It was 23 kJ/mol higher in energy than the best HIE structure.A similar structure in which NH 2 binds only to Fe2, but with a hydrogen bond from NH 3 to NH 2 (1.74 Å) is 100 kJ/mol less stable than the best HIE structure (both structures are shown in the right-hand side of Fig. 6. As for the NNH 3 structure, the reaction energy for the dissociation of NH 3 from the HIENH 2 + NH 3 complex is unfavourable by 105 kJ/mol.However, if another electron and proton are added, this energy goes down to 21 kJ/mol, which could be overcome by entropic effects from the released NH 3 molecule.In fact, NH 3 can be released from Fe6 in this complex with a small barrier (14 kJ/mol) and a favourable reaction energy (16 kJ/mol; giving NH 3 hydrogen bonding to Cys-275, S1A and S2A).

E 5 -E 8 structures with a single N atom
Finally, we studied structures with only one N atom bound.To make the study complete, we included also the E 5 -E 6 states, even if our results favour an alternating mechanism (which involves a single N atom only in the E 7 and E 8 states, cf.Fig. 2).The structure of the N 3-complex (E 5 ; e 0 n 1 ) is shown in Fig. 7  forms a hydrogen bond to .It was found to be most stable in the quartet state with two negative spins on Fe1 and Fe5 (BS-15).The spin population on Fe6 is still very low, 0.4.The corresponding HID structure has also Fe-N distances of 1.81 and 1.72 Å.The hydrogen-bond distance is 1.70 Å. NE2 of His-195 is 3.10 Å from the hydrogen atom of NH.It is 59 kJ/mol less stable than the HIE complex with TPSS, but actually 18 kJ/mol more stable with B3LYP.The corresponding HIP-N complex is 125 kJ/mol less stable than the HIE complex (160 kJ/mol with B3LYP), showing that protonation of N 3-is favourable.
Adding another electron and proton to the NH complex, gives the E 7 NH À 2 complex (e 2 h 2 n 1 ), which is also shown in Fig. 7 (with HIE).The ligand is still bound to both Fe ions with Fe-N bonds of 1.92 and 1.96 Å. NH 2 forms a hydrogen bond to Gln-191 (1.76 Å).The structure is most stable in the triplet state (BS7-235), but the quintet state is close in energy (4 kJ/mol less stable).Fe6 has a rather low spin population (1.6), whereas Fe2 has a normal spin population (-2.8) compared to the other Fe ions (2.3-3.2).The corresponding HID complex is similar, with Fe-N bonds of 1.91 and 1.95 Å and the same hydrogen-bond length.It is 55 kJ/mol less stable than the HIE structure, but it is 75 kJ/mol more favourable than the HIP-NH structure (54 and 118 kJ/mol with B3LYP, respectively).
Next, we added the eighth electron and proton to the NH À 2 complex, giving a NH 3 complex (e 3 h 3 n).It is shown in the right-hand side of Fig. 7 (with HIE).NH 3 is still symmetrically bound to Fe2 and Fe6 with Fe-N distances of 2.18 and 2.21 Å.The hydrogen bond to Gln-191 is 1.58 Å.It was found to be most stable in the doublet state with BS10-147.The spin population on Fe2 (2.3) is only slightly smaller than that on the other Fe ions (2.5-3.2).The corresponding HID complex is similar, with Fe-N bonds of 2.08 and 2.33 Å, and a hydrogen-bond length of 1.60 Å.However, the NH 3 group forms a hydrogen bond also to His-195, 1.83 Å.Therefore, it is only 5 kJ/mol less stable than the HIE complex.It is 45 kJ/mol more stable than the HIP-NH À 2 complex (19 and 64 kJ/mol with B3LYP), i.e. protonation of NH À 2 is still downhill.This is the only complex in this section that is involved in our preferred alternating reaction path.
Finally, we tried to dissociate NH 3 from this E 8 state, which was not facile.The reaction energy for the dissociation of NH 3 was strongly unfavourable, by 160 and 195 kJ/mol for the HIE and HID structures, respectively.The dissociation energy of NH 3 was only partly improved if the complex was protonated (i.e. to the e 3 h 4 n 1 state), 150 kJ/mol.If it was also reduced (e 4 h 4 n 1 ), the reaction energy decreased to 111 kJ/mol, which is still too large.The problem is most likely that that the dissociation of NH 3 leaves Fe2 and Fe6 coordinatively unsaturated.Therefore, we investigated whether the reaction energies were improved if SH -binds when NH 3 dissociates.We then assumed that SH -comes from the pocket formed by the flip of Gln-191 and that Gln-191 rotates back to its original position when S2B rebinds (it gains 67 kJ/mol by that back-rotation in our calculations).Quite satisfactory, this changed the reaction energy for dissociation of NH 3 , accompanied by the rebinding of SH -to the cluster (giving the protonated e 0 h 1 state) so that it became favourable by 153 and 177 kJ/mol for HIE and HID, respectively.

Conclusions
We have investigated putative reaction paths for N 2 replacing S2B in Mo-nitrogenase using QM/MM calculations, inspired by two recent crystal structures, showing that the S2B group may reversibly dissociate from the catalytic cluster [8,39].We have studied all possible protonation states of the N 2 -derived ligand and of His-19 5 for the E 4 -E 8 intermediates.Based on these calcula-tions, we suggest the reaction mechanism summarised in Fig. 8.The reaction starts with N 2 binding to the empty S2B site, between Fe2 and Fe6, in an irregular half-bridging mode.N 2 is then protonated by the two protons, already present in the cluster in the E 4 state.The first protonation is slightly upwards (17 kJ/mol), whereas the second protonation is strongly favourable by 148 kJ/mol.However, it should be noted that these energies are very uncertain, as we did not perform any full investigation of all (~2500) possible positions of the two protons in the FeMo cluster [43].On the other hand, our calculations quite conclusively suggest that the most stable N 2 H 2 structure has NNH 2 bound end-on with one N atom binding Fe2 and Fe6 (Fig. 3; rather than side-on binding HNNH).
For the E 5 state, many different structures are possible and we find that the H 2 NNH structure in Fig. 4 is most stable.It shows an irregular structure with the singly protonated N atom bridging Fe2 and Fe6, whereas the doubly protonated N binds only to Fe6.
For the E 6 state, a rather symmetric H 2 NNH 2 (hydrazine) structure is most stable (Fig. 5).When it is reduced and protonated to the E 7 state, the N-N bond is easily cleaved, giving NH À 2 bridging Fe2 and Fe6, whereas NH 3 is bound to only Fe2 (Fig. 6).Finally, after the final reduction, NH 3 dissociates and NH À 2 is protonated to NH 3 .The latter group cannot dissociate until S2B rebinds to the cluster.
In most of these structures, the protonated N 2 -ligand forms hydrogen bonds to the OE2 atom of the rotated Gln-191 residue.The sidechain of His-195 is also close to the ligand and we have tested three possible protonation states of it, HID, HIE and HIP.In all cases, the HIP structure is strongly unfavourable and may easily protonate the ligand, suggesting a possible protonation path of the ligand.However, as has been discussed by Dance,[80] it is unlikely that His-195 may provide more than one proton to the substrate.It would also give the HID state after the proton transfer, whereas we find that the HIE state is preferred for all structures in the suggested reaction path.
Thus, we suggest a mechanism that is mainly alternating: For E 5 -E 7 , the protons are added alternatingly to the two N atoms.Moreover, the calculations suggest that hydrazine (H 2 NNH 2 ) is a reaction intermediate and the first NH 3 product does not dissociate until the E 7 state.However, for the E 4 intermediate, we actually suggest that the NNH 2 structure, rather than HNNH, is most stable, although the former intermediate traditionally is connected to the distal reaction mechanism.
In the latter aspect, our calculations agree with those in the previous study by Nørskov and coworkers [20].However, they then suggested a distal mechanism, forming NNH 3 after the addition of the third proton.Our calculations agree that such a mechanism is also possible.However, for the key E 5 state, we find that the H 2 NNH intermediate is significantly more stable.It is likely that Nørskov did not consider all possible N 2 H 3 intermediates and their QM-cluster model was quite small, excluding the crucial Gln-191 and His-195 residues.Therefore, we strongly believe that the present QM/MM results with a large QM system and the explicit consideration of the surrounding protein and solvent are more reliable than Nørskov's QM-cluster calculations with a minimal cluster model of the active site.
Kästner and Blöchl have suggested a mechanism of nitrogenase in which S5A only dissociated from Fe7 but remained bound to Fe3 [32,81].The calculations were based on an incorrect model with a central nitride ion, without considering the formation of H 2 (N 2 binds to the E 2 state) and were performed on a minimal QMcluster without any account of the surrounding.They suggested that N 2 binds bridging Fe3 and Fe7 after dissociation of the protonated S5A from Fe7 (but not from Fe3), and follows a mainly alternating mechanism.However, they suggested that HNNH and HNNH 3 are involved in the mechanism, rather than NNH 2 and H 2 NNH 2 in our mechanism.Again, we believe that our calculations with a correct FeMo model and full account of the surroundings are more accurate.
Recently, we have considered possible structures for N 2 H 2 binding to the intact FeMo cluster (i.e. with S2B still bound and Gln-191 in the non-flipped orientation) [48].The calculations still indicated that Fe2 and Fe6 are the most probable binding sites for N 2 H 2 .However, the most stable structure had trans-HNNH bound terminally to Fe2.This is in contrast to the present result, where end-on bound NNH 2 is most stable.Structures involving NNH 2 were also tested for the intact cluster, but they were found to be at least 28 kJ/mol less stable than the best structure.However, a structure with HNNH 2 bound to Fe6, where the third proton was transferred automatically from the homocitrate ligand, was quite low in energy, 3-29 kJ/mol less stable than the best structure.
In conclusion, we have shown that Mo-nitrogenase may bind N 2 in the empty site formed by dissociation of the S2B ligand and this molecule can then be reduced and protonated to two molecules of ammonia along a favourable alternating mechanism, which is 51 kJ/mol more favourable than a distal mechanism at the E 5 intermediate.However, this does of course not prove that this is the actual mechanism followed by Mo-nitrogenase.To that end, similar studies are needed for the corresponding reaction without dissociation of S2B (or with S2B dissociated from only one of Fe2 or Fe6 as has recently been suggested both computationally and experimentally [79,82,83]), as well as detailed studies of the dissociation mechanism and the binding of N 2 .These will be investigated in future studies.
(HIE state) and it has a similar asymmetrical coordination with Fe2-N = 1.95 Å, Fe6-N = 2.00 Å for the singly protonated N atom and Fe6-N = 1.93 Å for the doubly protonated N atom.It has a short hydrogen bond between the NH 2 group and Gln-191 (1.54 Å).His-195 also points towards the NH group with a H . . .N distance of 2.62 Å.It is 53 kJ/mol more stable than the other HIE-HNNH 2 structure.The corresponding HID structure has Fe2-N = Fe6-N = 1.92 Å for the singly protonated N atom and Fe6-N = 1.96Å for the other N atom, H . . .O E1 = 1.53 Å and H . . .N E2 = 3.37 Å (opposite polarity), but it is 39 kJ/mol less stable than the HIE structure (55 kJ/mol with B3LYP) but 8 kJ/mol more stable than the HID-HNNH 2 structure.The end-on NNH 3 complex was most stable in the quintet state.It has still an intact N-N bond of 1.40 Å (HIE state).The two Fe-N bonds are almost unchanged compared to the NNH 2 structure (1.78 and 1.81 Å).It has a strong hydrogen bond to Gln-191 (1.38 Å).In fact, when His-195 is instead protonated on the ND1 atom, the NNH 3 complex is stabilised by 32 kJ/mol, owing to a favourable hydrogen bond between one of the protons on NNH 3 and His-195 (1.66 Å; shown in the right side of Fig.

Fig. 8 .
Fig. 8. Suggested reaction mechanism of Mo-nitrogenase with a dissociated S2B, based on the QM/MM calculations.
[12,21,49,77]  4oxidation state[12,21,49,77].Inthis study, we have considered the E 4 -E 8 states, after the dissociation of H 2 and binding of N 2 .To clarify exactly what state has been studied, we will use a more detailed nomenclature.We denote each state by e i h j n k , where i is the number of extra electrons, j is the number of extra protons and k is the number of extra nitrogen atoms relative to a E 0 resting state with S2B dissociated (i.e.MoFe 7 S 8 C in the Mo III Fe II 3 Fe III 4 oxidation state, without any bound hydrogen atoms).The starting, N 2bound E 4 state is then e 2 h 2 n 2 , because four electrons and protons have been added, but two of each have dissociated again with H 2 .

Table 1
and HE2 of