{'fraction_non_alphanumeric': 0.04373522458628842, 'fraction_numerical': 0.04265007944812619, 'mean_word_length': 4.138092393468738, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Thermodynamically Favourable States in the Reaction of Nitrogenase without Dissociation of any Sulfide Ligand

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


Introduction
Nitrogen is an essential element of all lifeforms, being a component of all amino acids and nucleic acids. Although the atmosphere of Earth contains 78 % of N 2 , nitrogen is still a limiting element for plant growth and a prominent component of fertilizers. The reason for this is the strong triple bond in N 2 , which makes it chemically inert. [1,2] Industrially, N 2 is converted to ammonia through the Haber-Bosch process, which requires high temperature and pressure. [2] Only a single group of enzymes can cleave the NÀ N bond in N 2 , the nitrogenases (EC 1.18/19.6.1), which work at ambient temperature and pressure. [1,3,4] Crystallographic studies have shown that the most active type of nitrogenase contains a MoFe 7 S 9 C(homocitrate) cluster (the FeMo cluster) in the active site, connected to the protein by a histidine and a cysteine residue at the opposite ends of the cluster (Figure 1). [5][6][7][8][9] There also exist alternative nitrogenases with the Mo ion replaced with either vanadium or iron, which have lower activities towards N 2 . [10] The nitrogenases catalyse the reaction The mechanism is normally discussed in terms of nine intermediates E 0 -E 8 , differing in the number of added electrons and protons, according to the Lowe-Thorneley scheme. [11] Thorough biochemical, kinetic and spectroscopic studies have indicated that the resting E 0 state needs to be reduced to the E 4 state before N 2 may bind. [1,3,4,[12][13][14][15][16][17] It has also been suggested that H 2 formation through reductive elimination is a prerequisite for the binding of N 2 , explaining why H 2 is a compulsory byproduct in the reaction. It is normally assumed that N 2 is directly reduced and protonated to N 2 H 2 upon binding to the enzyme. [1,18] It has long been debated whether the nitrogenases follow a sequential or alternating reaction mechanism. In the sequential mechanism, the first three protons bind to the same N atom of N 2 , which then dissociates as NH 3 from the E 5 intermediate, before the second N atom starts to be protonated. This mechanism was originally suggested by Chatt and has gained support from inorganic model complexes. [19][20][21][22][23] In the alternating mechanism, the protons are instead added alternatively to the two N atoms, so that HNNH and H 2 NNH 2 (hydrazine) are intermediates and the first NH 3 product does not dissociate until the E 7 state. It 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,24,25] Moreover, it has been shown that N 2 , N 2 H 2 , CH 3 NH 2 and N 2 H 4 all react via a common intermediate. [1,26] The nitrogenases have been thoroughly studied also by computational methods. [1,13,35-42,27À 34] Unfortunately, these studies have given very diverging and disparate suggestions. In fact, there is not even any consensus about the structure of the key E 4 intermediate. Important reasons for this are that different density-functional theory (DFT) methods give very different predictions of the relative stability of various intermediates, with differences of 600 kJ/mol [43] and that there are very many possibilities for the structures and electronic states of the intermediates. [44,45] Hoffman and coworkers have suggested a structure of the E 4 intermediate with two hydride ions bridging the Fe2 and Fe6 ions, as well as the Fe3 and Fe7 ions, and with two protons on the S2B and the S5A sulfides, all positioned on the same face of the FeMo cluster [15,46] (the name of the various Fe and sulfide ions are shown in Figure 1b). They have shown that this structure is lower in energy than a few other structures and that it may bind N 2 after reductive elimination of the two hydride ions, leaving the cluster in a doubly reduced state. [15,41,47,48] On the other hand, Siegbahn has suggested that the FeMo cluster needs to be reduced by four electrons from the resting state before the true E 0 state is reached, which involves a triply protonated central carbide and a strongly distorted cluster. [49] Then, this state is reduced by another four electrons to reach the E 4 state, from which H 2 dissociates and N 2 binds, bridging two Fe ions. It is successively protonated in a manner that is a mixture of the alternating and sequential mechanism, involving NNH 2 , HNNH 2 , H 2 NNH 2 , but also HNNH 3 . The first NH 3 dissociates at the E 7 level.
Dance has presented a mechanism in which E 4 contains two terminal hydride ions on Fe2 and Fe6, and two protons on S2B and S3B. N 2 then binds side-on to Fe6, without any dissociation of H 2 and is alternatively protonated to H 2 NNH 2 , at which level the NÀ N bond is cleaved, forming two NH 2 fragments on Fe2 and Fe6. [27,50,51] On the other hand, Nørskov and coworkers have suggested a mechanism in which the E 0 state is doubly protonated and a sulfide ligand dissociates from the cluster during the reaction mechanism. [37] This forms a binding site, where N 2 binds in an end-on fashion, bridging two Fe ions and it is then sequentially protonated on the outer N atom. The dissociation of the sulfide ion was inspired by several crystallographic studies of both Mo and V nitrogenase, showing that the S2B group can be replaced by several other ligands, for example CO, OH À and Se. [8,[52][53][54][55][56] Recently, we studied a similar mechanism, involving dissociation of S2B. [57] We used a larger and more realistic model system, which was studied with the combined quantum mechanical and molecular mechanics (QM/MM) approach. Our study indicated that the conversion of N 2 H 2 to two NH 3 molecules is thermodynamically favourable, but it follows a mainly alternating pathway (although the first intermediate involved a bridging NNH 2 group, which is normally connected to a sequential mechanism). Naturally, such studies do not prove that the nitrogenase mechanism actually involves a dissociated S2B group. To that end, it must be shown that the replacement of S2B by N 2 is energetically favourable, which has been questioned by Dance. [58] Moreover, it should be shown that a reaction without replacement of S2B is not possible or at least is less favourable. Here, we make an investigation of the latter reactions, that is, the formation of ammonia from bound N 2 for Mo nitrogenase without dissociation of S2B. We show that also such a reaction is possible and thermodynamic favourable, following an alternating mechanism.

The protein
The 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. [43,[59][60][61] The entire heterotetramer was considered in the calculations, because the various subunits are entangled without any natural way to separate them. The quantum mechanical (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, using a QM charge model. [59] The protonation states of all residues were the same as before: [59] 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 hydroxy 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 quantum-refinement study [59] and this protonation state is also supported by another QM/MM study. [62] 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 [59] ) to neutralise the protein and give an ionic strength of 0.2 M. [63] 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. [59] All MM calculations were performed with the Amber software. [64] For the protein, we used the Amber ff14SB force field [65] and water molecules were described by the TIP3P model. [66] For the metal sites, the MM parameters were the same as in our previous investigation. [59] The metal sites [45,59] were treated by a non-bonded model [67] and charges were obtained with the restrained electrostatic potential method, obtained at the TPSS/def2-SV(P) level of theory [68,69] and sampled with the Merz-Kollman scheme. [70] The FeMo cluster was modelled by MoFe 7 S 9 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 model, viz. Arg-96, Gln-191 and His-195 (sidechains), Ser-278 and Arg-359 (both backbone and sidechain, including the  Cα and C and O atoms from Arg-277), Gly-356, Gly-357 and  Leu-358 (backbone, including the Cα and C and O atoms from  Ile-355), as well as two water molecules. Finally, Phe-381 and Val-70 were also included because they are close to the putative N 2 binding site and therefore may affect the binding of the substrate. The QM system involved 183-190 atoms in total (depending on the number of added protons and N atoms) and is shown in Figure 1a. The net charge of QM region was À 3.

QM calculations
All QM calculations were performed with the Turbomole software (version 7.5). [71] All structures were studied with both the TPSS [68] and B3LYP [72][73][74] functionals with def2-SV(P) basis set. [69] The most stable states were examined also with the larger def2-TZVPD basis set. The calculations were sped up by expanding the Coulomb interactions in an auxiliary basis set, the resolution-of-identity (RI) approximation. [75,76] Empirical dispersion corrections were included with the DFT-D4 approach, [77] as implemented in Turbomole. All minima were fully optimized without any restraints. Transition states (for NÀ N cleavage and NH 3 dissociation) were determined as the highest point on the potential energy surface along the reaction coordinates, which were scanned with a step of 0.1 Å near the transition states.
Experiments have shown that the ground spin state of E 4 is a doublet, [1,78] and we used this state for E 4 models. For the other E n states, we used mainly the doublet or triplet states, but for the most interesting structures, 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: [79] Each of the seven Fe ions were modelled in the high-spin state, with either a surplus of α (four Fe ions) or β (three Fe ions) spin. Such a state can be selected in 35 different ways . [60] The various BS states were obtained either by swapping the coordinates of the Fe ions [80] or with the fragment approach by Szilagyi and Winslow. [81] The various BS states are named by listing the number in the Noodleman nomenclature (BS1-10), [79] followed by the numbers of the three Fe ions with minority spin (however, in the tables, only the latter three numbers are given). Most structures were studied in the BS10-147 state, i. e. with β spin on Fe1, Fe4 and Fe7, because it was found to be lowest in both this and in our previous study. [57] However, sometimes the calculations converged to other states (especially BS7-235). For twelve of the most stable structures, the relative stabilities of all 35 states were examined (with structures fully optimised for each BS state). Moreover, for all structures within 20 kJ/mol of the most stable structure at each E n level, the BS7-235 state was also studied.
As have been discussed before, [43,60] TPSS/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 distances, and 0.02 and 0.06 Å for metal-ligand distances, and a root-mean-squared-deviation (RMSD) of 0.06 Å for the metals and the first-sphere ligands. This is similar to the results obtained with the TPSSh [62] approach and appreciably better than with the B3LYP/def2-SV(P) method, which gives average and maximum deviations of 0.08 and 0.12 Å for the metalmetal and 0.04 and 0.11 Å metal-ligand distances, respectively and a RMSD of 0.08 Å. Therefore, we discuss primarily the TPSS/ def2-SV(P) results.

QM/MM calculations
QM/MM calculations were performed with the ComQum software. [82,83] In this approach, the protein and solvent are split into three subsystems: System 1 (the QM region) was relaxed by QM methods. System 2 contained all residues and water molecules with at least one atom within 6 Å of any residues and 39 water molecules). Finally, system 3 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 nonperiodic with 133 915 atoms. Most calculations were performed without relaxing system 2, but for the most interesting structures, calculations with relaxed surroundings were also performed. The effect of the relaxed surroundings are described in the Supporting Information.
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. [82,84] All atoms were included in the point-charge model, except the CL atoms. [85] The total QM/MM energy in ComQum is calculated as [82,83] in which E HL QM1þptch23 is the QM energy of the QM system truncated by HL atoms and embedded in the set of point charges modelling systems 2 and 3 (but excluding the selfenergy 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 MM123;q 1 ¼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 van der Waals link-atom corrections. [86] No cut-off is used for any of the interactions in the three energy terms in Equation (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.
QM/MM calculations give comparable energies only if they contain exactly the same number of electrons and atoms of each element in both the QM and MM systems. Therefore, we compare only structures within the same E n level. On the other hand, it means that we can study proton transfers within the QM system, for example from the homocitrate ligand or from His-195 to the substrate. For each transition from E n to E n + 1 , an electron and a proton is added to the QM system, and we compare the energies of structures with this proton in different positions.

Result and Discussion
In this investigation, we study the later part of the reaction mechanism of nitrogenase, assuming that the S2B ligand does not dissociate. We describe in separate sections states at different oxidation levels, from E 4 to E 8 .

N 2 -bound E 4 structures
We start with the N 2 -bound E 4 state. As in our previous study of the reaction mechanism with a dissociated S2B ligand, [57] we avoid the problem there is no consensus regarding the protonation of the E 4 state [27,38,43,59,61,[87][88][89] by starting from a state where N 2 has already bound to the cluster and is protonated to N 2 H 2 . The immediate protonation of the substrate upon binding is normally assumed, [1,18] although it has not been experimentally observed. Mutations and other studies have shown that the substrate most likely bind either the Fe2 or Fe6 ions of the FeMo cluster [1,90,91] and this is also supported by a systematic scan of all possible N 2 H 2 binding positions. [45] Therefore, we look for the best structure with N 2 H 2 bound to either Fe2 or Fe6, or to both. In the latter case, the substrate may bind on two different sides of the bridging S2B ligand and we call these structures Fe2/6(3) and Fe2/6(5), depending on whether it is on the same side as S3A or S5A. Likewise, when N 2 H 2 binds only to Fe2, the non-bonding N atom point either to the S3A or S5A sides of the cluster, which will be called Fe2(3) or Fe2(5) (and similar for binding to Fe6). We have considered three isomers of N 2 H 2 , viz. NNH 2 , cis-HNNH or trans-HNNH (the latter two are abbreviated cHNNH or tHNNH in the following). Moreover, we have tested three protonation states of His-195: with protons on either ND1 (HID) or NE2 (HIE) or on both (HIP; but this adds an extra proton to the system and was therefore studied for the next higher E n level). The results are collected in Table 1 and are shown Figure S1.
It can be seen that in nearly all structures, the HIE protonation state was more favourable than HID protonation by 21-132 kJ/mol. The only exception was the Fe2/6(3)-cHNNH state with B3LYP. All structures were studied in the BS10-147 state, but sometimes it shifted to the BS7-235 state (the spin on Fe6 is often small and may change sign; the latter state was also studied for all low-energy structures). However, BS10-147 was always 9-37 kJ/mol lower in energy than BS7-235 when both states were found, except for the Fe6-cHNNH(3) state, for which BS7-235 was 1 kJ/mol more stable at both the TPSS and B3LYP levels).
The most favourable structure has trans-HNNH end-on bound to Fe6 (Figure 2a). The Fe6-N distance is 1.91 Å and the NÀ N bond length is 1.26 Å, which is slightly longer than in isolated trans-HNNH, optimised with the same level of theory, 1.25 Å. The Mulliken spin populations are (in absolute terms) 3.2-2.7 e on the seven Fe ions, except Fe6, which has only 1.6 e (cf. Table S1 in the Supporting Information). Mo has a population of À 0.3 e. This structure is stabilised by a hydrogen bond from the HN group bound to Fe6 to the alcohol O7 atom of homocitrate (which coordinates to Mo), with a H···O distance of 1.95 Å (cf. Figure 2a). The other N atom of the substrate receives a hydrogen bond from the HE1 atom of Gln-191 (the HE1···N distance is 2.17 Å) and the other H atom of the substrate is directed towards S2B, with a H···S distance of 2.29 Å, but the N-H···S angle is only 124°. S2B receives another hydrogen bond from the HE2 atom of His-195 (2.13 Å, with a more ideal geometry). The corresponding structure in the quartet state is 30 kJ/mol less stable with the TPSS functional, but 7 kJ/mol

Chemistry-A European Journal
Research Article doi.org/10.1002/chem.202103933 more stable with the B3LYP. A full investigation of all 35 BS states, collected in Table S2, shows that the BS10-147 state is indeed the most stable BS state, 10-67 kJ/mol more stable than the other BS states (BS7-235 second lowest).
A structure with NNH 2 bound end-on to Fe6 is only 4 kJ/mol less stable with TPSS (also with the larger def2-TZVPD basis set and 3 kJ/mol with relaxed surroundings) and it is actually 6 kJ/ mol more stable with B3LYP. In this structure, the substrate has abstracted the alcohol proton from homocitrate, giving HNNH 2 (Figure 2b). It has a Fe6-N distance of 1.85 Å and a NÀ N bond length of 1.29 Å (1.21 Å in neutral NNH 2 and 1.23 Å in HNNH 2 + optimised with the same method). It is stabilised by hydrogen bonds from the Fe6-bound NH group to the alcohol O7 atom of homocitrate (H···O7 = 1.72 Å), from the NH 2 to one of the carboxylic groups of homocitrate (H···O1 = 1.67 Å) and from the other atom of NH 2 group to S2B (H···S2B = 2.31 Å; Figure 2b). The S2B atom also receives a hydrogen bond from His-195 (HE2···S2B = 2.12 Å) and the O1 atom of homocitrate receives another hydrogen bond from Gln-191 (HE1···O1 = 1.93 Å). This structure has also a low spin population on Fe6 (1.8 e). Again, a full investigation of all 35 BS states showed that BS10-147 is the most favourable BS (Table S2), although it is only 4 kJ/mol more stable than BS10-135 (14 kJ/mol more stable than BS7-235).
A structure with trans-HNNH binding to Fe2, directed to the S3A side (Figure 2c), is 20 kJ/mol less stable than Fe6-tHNNH structure (22 kJ/mol with the larger basis set, 34 kJ/mol with relaxed surroundings and 5 kJ/mol with B3LYP). The Fe2-N distance is 1.95 Å and the NÀ N bond length is 1.25 Å. The two H atoms of the substrate point towards S2B (H···S2B = 2.72 Å) and S1A (H···S1A = 2.19 Å), but the NÀ H···S angles are far from straight (83°and 131°, respectively). The Fe spin population on Fe2 (2.1 e) is only slightly less than on the other Fe ions (3.2-2.4 e in absolute terms). A full investigation of all 35 BS (Table S2) shows that BS10-147 indeed is the most stable state. The structure with the substrate directed to the S5A side is 41 kJ/ mol less stable.
Other structures studied are appreciably less stable. Structures with cis-HNNH bound side-on to Fe2 and Fe6 are 161- Table 1. Energies (in kJ/mol), NÀ N and FeÀ N distances (in Å) of the various structures of the E 4 states. All states were studied in the S = 1 = 2 state unless otherwise stated. His is the protonation state of His-195. The BS state is identified by the Fe ions with minority spin (e. g. 147 = BS10-147). Up to four different energies are listed: TPSS-D4/def2-SV(P) (TP), B3LYP-D4/def2-SV(P) (B3), TPSS-D4/def2-TZVPD (TZ), all based on TPSS-D4/def2-SV(P) structures obtained with fixed surroundings, as well as TPSS-D4/def2-SV(P) with relaxed surroundings (Rlx). When multiple FeÀ N distances are given, Fe2 comes before Fe6. We also considered structures with S2B dissociated from either Fe2 or Fe6 (but still bound to the other ion), because such structures have been suggested to be competitive by other authors. [35,42,48,58] However, with our methods, such structures were always high in energy, by 108-255 kJ/mol for structures with NNH 2 bridging Fe2 and Fe6 with the terminal N atom and by 152-204 kJ/mol for structures with cis-HNNH bridging Fe2 and Fe6.
The present results are somewhat different from those obtained in our previous study of the binding of N 2 H 2 to the FeMo cluster. [45] In that study, we found that the structure with trans-HNNH bound to Fe2 was 10 kJ/mol more stable than the Fe6-tHNNH state and 9 kJ/mol more stable than the Fe6-HNNH 2 state (19 and 3 kJ/mol with the larger def2-TZVPD basis set). The difference is most likely connected to the larger QM system used in the present study (models of Val-70, Gln-191 and Phe-381, all situated around the binding site, were not included in the previous study). The present results should be more reliable.
The results also differ from those obtained with a dissociated S2B ligand (and a rotated conformation of Gln-191), [57] for which a structure with NNH 2 bridging Fe2 and Fe6 was found to be most favourable. Such structures were at least 108 kJ/mol less stable than Fe6-tHNNH in this study and led to half-dissociation of S2B. Clearly, the active site with a bridging S2B group is so crowded that it disfavours structures with N 2 H 2 simultaneously bridging Fe2 and Fe6.
Bjornsson and coworkers studied structures with trans-HNNH bound to the FeMo cluster in the E 4 state. [42] They also found that binding to Fe6 was more favourable than to Fe2, in agreement with our results, but the difference was larger, 69 kJ/ mol. In their models, S5A was protonated and one of the protons on HNNH was abstracted from homocitrate.

E 5 structures
Next, we added an electron and a proton to the FeMo cluster (i. e. to the QM system) to obtain structures at the E 5 level. They were studied in the triplet state with BS10-147. The structures are described in Table 2 and the best are shown in Figures 3  and S2. As for the E 4 structures, HIE structures were always more stable than the corresponding HID structures, by 16-173 kJ/mol.
The most favourable state has H 2 NNH 2 (hydrazine) bound to Fe6, where the extra proton is abstracted from the hydroxy group of homocitrate (Figure 3a). The Fe6-N distance is 2.09 Å and the NÀ N bond length is 1.43 Å, which is the same as for isolated hydrazine, optimised at the same level of theory. The non-coordinating NH 2 group is directed to the S3A side of the cluster. The two H atoms of the Fe-bound NH 2 group forms hydrogen bonds to O7 of homocitrate (H···O7 = 1.94 Å) and S3B (H···S3B = 2.73 Å). The other two H atoms of hydrazine point towards O1 of homocitrate (H···O1 = 2.66 Å) and S2B (H···S2B = 2.34 Å). The spin density on Fe6 is 2.1 e (Table S3). The singlet state was 16 kJ/mol more stable than the triplet state with TPSS, but 38 kJ/mol less stable with B3LYP. The quintet was 4-33 kJ/mol less stable than the triplet. A full investigation of all BS states (Table S2) showed that BS7-235 is actually 9 kJ/mol lower in energy than BS10-147. In fact, nine different BS states were found within 11 kJ/mol of the lowest state. There are several structures with similar energies with slight variations in the hydrogen-bond lengths and the relative conformations of the two NH 2 groups. For example, a structure with the nonbonded NH 2 group directed towards S5A is only 1 kJ/mol less stable (6 kJ/mol more stable by B3LYP, but 38 kJ/mol less stable with relaxed surroundings).
Other structures are appreciably less stable. The secondbest structure had HNNH 3 bound end-on to Fe6 with the NH group and with the NH 3 group directed toward S5A (Figure 3b; again with a proton abstracted from homocitrate). It is 54-68 kJ/mol less stable than the Fe6-H 2 NNH 2 structure at the

Chemistry-A European Journal
Research Article doi.org/10.1002/chem.202103933 various levels of theory. A similar structure with the NH 3 group directed towards S3A is 10-18 kJ/mol less stable. A structure with HNNH 2 bound end-on to Fe6 with the NH group (i. e. without the proton transfer from homocitrate) is 66-84 kJ/mol less stable than the Fe6-H 2 NNH 2 structure. A structure with HNNH 2 bound to Fe2 with the NH 2 group directed towards S3A is 54-94 kJ/mol less stable than the best structure. The structure with the NH 2 group pointing in the opposite direction is 50 kJ/ mol less stable. Structures with NNH 3 bound to Fe2 are 184-196 kJ/mol less stable than the best Fe6-HNNH 3 structure.
We also studied a number of structures with HNNH 2 or NNH 3 bridging Fe2 and Fe6 on either side of S2B, but all of them were high in energy (182-471 kJ/mol less favourable than Fe6-HNNH 3 ). Moreover, we studied structures with a proton on His-195, instead of on the substrate, giving the HIP state. These structures were at least 233-430 kJ/mol higher than Fe6-H 2 NNH 2 -HIE, with Fe6-HNNH 2 lowest, 6 kJ/mol lower than Fe6-tHNNH. This may illustrate a possible path for the transfer of protons to the substrate and it is apparently strongly downhill.
For the best Fe6-HNNH 3 structure, we tested to cleave the NÀ N bond. However, this reaction turned out to prohibitive with an activation barrier of 119 kJ/mol and a reaction energy of 78 kJ/mol to a product with NH bound to Fe6 and NH 3 dissociated from the cluster, but forming hydrogen bonds to O1 of homocitrate, S2B and S3B. This indicates that nitrogenase does not follow a sequential reaction mechanism. On the other hand, the NÀ N bond can be cleaved in Fe2-NNH 3 , with a barrier of 49 kJ/mol, but the reactant and the Fe2-N product are 196 and 149 kJ/mol less stable than the Fe6-H 2 NNH 2 structure, showing that they are not expected to form during the reaction mechanism.

Chemistry-A European Journal
Research Article doi.org/10.1002/chem.202103933 Adding a proton and an electron to the previous structures gives intermediates at the E 6 level. These were studied primarily in the doublet BS10-147 state. The results are collected in

Chemistry-A European Journal
Research Article doi.org/10.1002/chem.202103933 The best structure at the TPSS level has H 2 NNH 2 bound endon to Fe6, with the non-coordinating NH 2 group pointing towards S3A (thus, the new proton in the E 6 state is added to homocitrate and not the substrate; Figure 4a). It has a Fe6-N distance of 2.10 Å and a NÀ N bond length of 1.45 Å, which is slightly longer than for the E 5 structures and the same molecule optimised in vacuum, 1.43 Å. The Fe ions have spin populations of 2.2-3.2 e, but 2.0 e on Fe2 and 1.4 e on Fe6 (Table S4). The spin on Mo is minor and slightly positive, 0.1 e. It is stabilised by a hydrogen bond from the Fe6-bound NH 2 group to the hydroxy O7 atom of homocitrate (H···O7 = 2.25 Å) and by a hydrogen bond from the other NH 2 group to S2B (H···S2B = 2.31 Å). The other H atom of the latter NH 2 group points in the direction of S1B, but the distance is long, 3.02 Å. It is also 2.94 Å from O1 of homocitrate. The fourth H atom of the substrate is 2.73 Å from S3B, but the geometry is far from ideal. It is also close to a methyl group of Val-70 (1.83 Å HÀ H distance). The quartet state was 16 kJ/mol less stable at the TPSS level, but 53 kJ/mol more stable with B3LYP. An investigation of all BS states (Table S2) showed that BS10-147 is most stable, 6 kJ/mol better than BS10-135. With B3LYP, a structure in which the two NH 2 groups are twisted with respect to each other (HÀ NÀ NÀ H torsion of 99°), as in the structure for the free hydrazine, is 38 kJ/mol more stable than the structure in Figure 4a, but at the TPSS level, the other structure is 43 kJ/mol more stable.
The corresponding structure with the non-coordinating NH 2 group pointing towards S5A (Figure 4b) is only 2-19 kJ/mol less stable with TPSS but 33 kJ/mol more stable B3LYP. It is stabilised by a hydrogen bond to the O1 atom of homocitrate (H···O1 = 2.04 Å), whereas the other H atoms interact with S2B, S3B and S1B as for the other conformation (H···S distances of 2.35, 3.13 and 2.79 Å). With B3LYP, this structure is further stabilised by 62 kJ/mol when studied in the BS7-235 state, whereas with TPSS the BS10-147 state is 48 kJ/mol more stable.
The corresponding structure with hydrazine bound end-on to Fe2 (Figure 4c) is 17-33 kJ/mol less stable. The two H atoms of the NH 2 group bound to Fe2 form hydrogen bonds to SG of Cys275 (H···SG = 2.43 Å) and S2B (H···S2B = 2.84 Å). The other two H atoms interact with ND1 of His-195 and S1A (H···ND1 = 2.53 and H···S1A = 3.15 Å), but with poor geometries. A full investigation of all BS states (Table S2) showed that this structure is most stable in the BS6-157 state, which is 16 kJ/mol more stable than the BS10-147 state. End-on bound HNNH 3 structures are 105-282 kJ/mol less stable, whereas side-on (Fe2/ 6) structures are 117-309 kJ/mol less stable and those with HNNH 3 dissociate to NH and NH 3 .
We have also studied the same structures as for E 5 , but with HIP and an extra electron. They were all high in energy, 156-547 kJ/mol less stable than the best Fe6-H 2 NNH 2 structure, showing that proton transfer from His-195 is strongly favourable. In several of this type of structures, the substrate automatically abstracts a proton from His-195, forming HID states instead (which are less stable than the HIE structures, as has already been discussed).

E 7 structures
After adding yet another proton and electron, we reach the E 7level intermediates. They were studied in the triplet BS10-147 state. The structures are listed in Table 4 and the most stable structures are shown in Figures 5 and S4. As usual, the HID structures were less stable than the corresponding HIE structures by 3-176 kJ/mol.
Only a few structures were obtained with H 2 NNH 3 bound to the FeMo cluster. They had a NÀ N distance of 1.42-1.43 Å. The most stable one had H 2 NNH 3 bound end-on to Fe6, with the NH 3 directed towards the S5A side (Figure 5a), but binding to Fe2 was only 4 kJ/mol less favourable with TPSS (30 kJ/mol by B3LYP). An investigation of all BS states (Table S2) showed that this complex is most stable in the BS2-234 state, which is actually 36 kJ/mol more stable than the BS10-147 state. The NÀ N bond can readily be cleaved in this structure with an activation barrier of only 32 kJ/mol and an exothermic reaction energy of À 154 kJ/mol. The product has NH 2 bound end-on to Fe6 and NH 3 dissociated, but hydrogen bonded to the cluster (Figure 5b). The two H atoms of NH 2 form hydrogen bonds to O1 of homocitrate (H···O1 = 2.47 Å) and S2B (H···S2B = 2.53 Å). The dissociated NH 3 molecule forms hydrogen bonds to NH 2 (H···N = 2.06 Å) and O2 of homocitrate (H···O2 = 2.23 Å), whereas the third H atom does not form any favourable interaction, but instead is quite close to a methyl group of Val-70 (2.03 Å HÀ H distance). A structure with H 2 NNH 3 dissociated from the Fe ions, but still hydrogen-bonded to the cluster is 83 kJ/mol more stable than the bound Fe6-H 2 NNH 3 structure (Figure 5c). When it is dissociated, it is appreciably harder to cleave the NÀ N bond -the calculated barrier is 91 kJ/mol.
Structures with a cleaved NÀ N bond and both NH 2 and NH 3 coordinated to the cluster are up to 84 kJ/mol more stable than the bound Fe6-H 2 NNH 3 structure, but 21 kJ/mol less stable than the structure with NH 3 dissociated. The best one has NH 2 bridging Fe2 and Fe6 on the side facing S5A, whereas NH 3 binds to Fe6 (Figure 5d). The two Fe-NH 2 distances are 1.93 and 1.95 Å, whereas the Fe6-NH 3 distance is 2.27 Å. S2B has moved considerably, but it still binds to Fe2 and Fe6, and it receives a hydrogen bond from His-195 (HE2···S2B = 2.36 Å). NH 3 forms a hydrogen bond to O1 of homocitrate (H···O1 = 2.00 Å), whereas the second H atom points towards S2B (H···S2B = 2.65 Å). The third H atom does not form any favourable interactions. The two H atoms of NH 2 point towards S3A and S2B (H···S3A = 3.21 and H···S2B = 2.63 Å). The spin populations on Fe2 and Fe6 are relatively low, 2.1 and 2.2 e, respectively, but that on Fe7 is even lower, 1.5 e (Table S5). NH 3 may dissociate from this structure, but the activation barrier is rather high, 78 kJ/mol.
There are several other structures with comparable energies (cf. Table 4), for example with NH 2 bridging Fe2 and Fe6 on either side of S2B and with NH 3 either on Fe2 or Fe6. The relative energies sometimes differ rather much between TPSS and B3LYP.
As for the other E n states, we tested also structures with HIP protonated, but all these were at least 183 kJ/mol less stable.

E 5 -E 8 structures with only one N atom
We studied also structures with only a single N atom, i. e. after NÀ N bond cleavage and dissociation of a NH 3 product. These were studied at four levels of oxidation and protonation (E 5 -E 8 ), even if the results in the previous subsections indicate that only the E 7 and E 8 states are involved in the reaction mechanism. The results are collected in Table 5 and the best structures are shown in Figures 6, 7 and S5.
The best E 5 structure has N bound end-on to Fe6 with a FeÀ N distance of 1.60 Å (Figure 6a). The N atom receives a hydrogen bond from HE1 of Gln-191 (2.55 Å, but this hydrogen also forms a hydrogen bond to O1 of homocitrate with a H···O1 distance of 2.40 Å). We tested also the singlet and quartet states for this structure. The latter was 38 kJ/mol less stable at the TPSS level (33 kJ/mol with B3LYP). However, the singlet was 39 kJ/mol more stable with TPSS, but 11 kJ/mol less stable with B3LYP. A structure with NH bound to Fe6 (with the proton abstracted from homocitrate) is only 2 kJ/mol less stable (38 kJ/ mol by B3LYP), but it is 14 kJ/mol more stable with the larger basis set and 9 kJ/mol more stable if the surroundings are relaxed.
The corresponding structure with N bound end-on to Fe2 (Figure 6b) is also 2 kJ/mol less stable (35 kJ/mol with B3LYP), but 1 kJ/mol more stable with the larger basis set, 41 kJ/mol more stable if the surroundings are relaxed. It has an even shorter Fe2-N bond length of 1.54 Å. The N atom does not receive any polar hydrogen bond, but it is 2.00 Å from a HB atom of Ser-278. The corresponding structures with N bridging Fe2 and Fe6 are 30 kJ/mol (on the S3A side) and 57 kJ/mol (S5A side) less stable. In both cases, S2B moves to a position where it interacts with more Fe ions than Fe2 and Fe6. Moreover, N receives the hydrogen bond from His-195 (instead of S2B; 2.44 and 2.50 Å, respectively). The corresponding HID structures are 78-126 kJ/mol less stable. Interestingly, for the three most stable structures BS7-235 was found to be 5-37 kJ/mol more Table 4. Energies (in kJ/mol), NÀ N and FeÀ N distances (in Å) of the various structures of the E 7 states. All states were studied in the S = 1 state unless otherwise stated. The entries are the same as in Table 1 stable than the BS10-147 state. However, the Fe6-N structure was still the best structure. The most stable E 6 structure has NH 2 bound end-on to Fe6 with the extra proton on the substrate abstracted from homocitrate ( Figure 6c). It has a Fe6-N bond length of 1.85 Å and the H atoms of NH 2 point toward O1 and O7 of homocitrate (H···O1 = 2.01 Å, H···O7 = 2.17 Å) and S2B (H···S2B = 2.59 Å). The spin population on Fe6 (2.0 e) is slightly lower than that on Fe2 (2.1 e; Table S6). This structure is 126 kJ/mol more stable (63-90 kJ/mol with B3LYP, the larger basis set or relaxed surroundings) than a structure with NH binding end-on to Fe6 (Fe6-N bond length of 1.73 Å). The corresponding HID structures are 51-221 kJ/mol higher in energy. Structures with HIP are 305-373 kJ/mol higher in energy, showing that proton transfer from His-195 is favourable.
The most stable E 7 structure has NH 3 bound end-on to Fe6 (with a proton abstracted from homocitrate). The Fe6-N bond length is 2.06 Å (Figure 6d). One of the three H atoms forms hydrogen bonds to both O1 and O7 of homocitrate (H···O1 = 2.01 and H···O7 = 2.03 Å). The other two are rather close to S2B (H···S2B = 2.76 and 2.86 Å), but also to S1B and S3B (H···S1B = 2.96 Å and H···S3B = 3.09 Å). The spin population of Fe6 is 2.1 e (Table S6). NH 3 cannot dissociate from this structure (the energy keeps rising by more than 120 kJ/mol when the Fe6-N bond is elongated). A full investigation of all BS states (Table S2) showed that BS6-156 is actually lowest, but the BS10-147 state is only 8 kJ/mol higher in energy and there are five BS states within 10 kJ/mol of BS6-156. This Fe6-NH 3 structure is 99-113 kJ/mol more stable than a structure with NH 2 bound to Fe6 (and a proton on homocitrate). The structure with NH 2 bound to Fe2 is 95-119 kJ/mol less stable. Structures with NH 2 bridging Fe2 and Fe6 are 191-350 kJ/mol higher in energy than the best structure. Structures with HID are~90 kJ/mol less stable than the corresponding HIE structures. Structures with HIP are 281-418 kJ/mol less stable, again indicating that proton transfer from His-195 is strongly favourable.
Finally, we studied also E 8 states. Again, the most stable structure has NH 3 bound to Fe6 with a FeÀ N distance of 2.13 Å (Figure 7a). NH 3 forms a hydrogen bond to O1 of homocitrate with a HÀ O distance of 1.93 Å. One of the other H atoms is 2.68 Å from S2B, but the last atom does not form any  (Table S6). The structure was studied in the doublet state and the corresponding quartet state is 23 kJ/mol less stable with TPSS, but 35 kJ/mol more stable with B3LYP. A Table 5. Energies (in kJ/mol), NÀ N and FeÀ N distances (in Å) of the various structures of the E 5 -E 8 states with a single N atom. All states were studied in the S = 1 (E 5 or E 7 ) or S = 1 = 2 (E 6 or E 8 ) state unless otherwise stated. The entries are the same as in Table 1 [92,93] A structure, in which NH 3 has abstracted the proton from O7 of homocitrate, forming NH 4 + , is actually 2 kJ/mol more stable than the Fe6-NH 3 structure. NH 4 + has dissociated from Fe6 and the four H atoms form hydrogen bonds to O1 of
A structure with NH 3 bound to Fe2 is 8-11 kJ/mol less stable than the Fe6-NH 3 structure (but 28 kJ/mol more stable with relaxed surroundings). It has a slightly shorter FeÀ N bond (2.06 Å; cf. Figure 7c). Two of the H atoms of NH 3 points approximately towards SG of Cys-275 and S2B (H···SG = 2.75 Å and H···S2B = 2.89 Å). A full investigation of all BS states showed that BS10-147 indeed is the most stable state, but only 2 kJ/mol more stable than the BS6-157 state. We also tried to find structures with NH 3 bridging Fe2 and Fe6, but they were at least 326 kJ/mol less stable. Structures with HID were 81-120 kJ/mol less stable than the corresponding HIE structures. Structures with HIP were 152-292 kJ/mol less stable, showing that proton transfer from His-195 is strongly favourable.

Conclusions
In this investigation, we have studied possible reaction intermediates of nitrogenase, assuming that the S2B remains bound to the FeMo cluster. To avoid the problem that the structure of the E 4 intermediate is not known and that different DFT functionals give very different relative stabilities of various protonation states, [43,44] we started our study after H 2 has dissociated and N 2 has bound to the cluster and has become doubly protonated to N 2 H 2 , [1,18] so that no protons remain bound to the cluster. Based on the accumulated experimental evidence, [1,90,91] as well a systematic study of the binding of N 2 H 2 to the FeMo cluster, [45] and in agreement with most previous computational studies, [27] we have assumed that N 2 H 2 binds either to Fe2 or Fe6.
Our study has led to the following conclusions: * For the E 4 state, Fe6-tHNNH, Fe2-tHNNH and Fe6-HNNH 2 structures are all competitive (within 5 kJ/mol with at least one of the four levels of theory included in Table 1). Based on these results, we suggest the reaction mechanism in Figure 8. In this mechanism, the substrate binds to Fe6. In the E 5 state, it is protonated to H 2 NNH 2 , whereas in E 6 , the proton is added to homocitrate, so that the ligand remains H 2 NNH 2 . In the E 7 state, the substrate is protonated to H 2 NNH 3 , in which the NÀ N bond is readily cleaved and NH 3 automatically dissociates. The resulting NH 2 group remains bound to Fe6 and is protonated to NH 3 . In the E 8 state, NH 3 dissociates and the resting E 0 state is formed.
It has been much discussed whether the reaction mechanism of nitrogenase follows a sequential or alternating reaction mechanism. [1] Our results suggest that the enzyme follows an alternating mechanism, with HNNH and H 2 NNH 2 as intermediates (although the former may be protonated by homocitrate to HNNH 2 ). Moreover, the NÀ N bond is cleaved in the E 7 state and that the NH 3 products dissociate at the E 7 and E 8 levels, in accordance with an alternating mechanism.
Between the various E n levels in Figure 8, we have assumed that an electron and a proton are delivered to the FeMo cluster. The electrons are provided by the Fe protein, via the Pcluster. [1,4] The protons ultimately come from the solvent. Dance Figure 8. Suggested reaction mechanism for nitrogenase, assuming that S2B remains bound to the cluster.

Chemistry-A European Journal
Research Article doi.org/10.1002/chem.202103933 has investigated different proton paths from the solvent to the FeMo cluster. He suggested that protons are transferred from the surface via a chain of water molecules to a water molecule close to homocitrate and the S3B atom. [94][95][96][97] He has also studied how the proton may be transferred within the FeMo cluster, starting on the S3B atom and ultimately ending up close to the substrate-binding sites. [96] His-195 may also provide protons to the substrate, as has often been assumed. [38,49] In fact, our calculations show that such a transfer (via the HIP state of His-195) is always downhill. However, this gives the HID state of His-195, which is nearly always unfavourable. Moreover, Dance has shown that rotation of the sidechain of His-195 is unlikely in the protein. [94,98] Therefore, His-195 can probably not provide more than one proton to the substrate (at the most).
However, it should be remembered that QM studies of nitrogenase are extremely complicated. Our survey of the most stable structures at the various E 4 -E 8 levels have indicated that several structures often have quite similar energies and that here are many conformations of the bound intermediates, depending on the direction of the non-coordinated N-group, the HÀ NÀ NÀ H dihedral angle and the formation of hydrogen bonds and other polar interactions. Moreover, sometimes other BS states than BS10-147 are most stable and a full BS investigation of every structure is currently too demanding. This makes it harder to settle the most stable structures and it cannot be excluded that we might have overlooked some lowenergy structures.
For all states, we have optimised structures with both TPSS and B3LYP. In general, the two methods give relative energies that agree reasonably. However, in some cases, differences > 50 kJ/mol are observed, without any significant differences in the geometry. This is often observed for the various spins states (B3LYP typically prefers higher spin states than TPSS, as expected). However occasionally, such differences also occur for structures with differences only in the hydrogen-bond pattern or coordination mode, making interpretation of the results harder.
For the most interesting states, we have also recalculated energies with the larger def2-TZVPD basis set. As in our previous studies, [43,[59][60][61] this has typically only a minor effect on the energies, < 20 kJ/mol. However, calculations with relaxed surroundings occasionally have larger effects, up to 26 kJ/mol, which may indicate that the structure need to relax more than is allowed in the rather large QM system or problems with local minima in the surroundings. However, a more detailed study of the relaxation of the surroundings (in the Supporting Information) indicates that the large energy differences are connected to major movements of water molecules and other groups far from the substrate-binding site, suggesting that it reflects more occurrences of multiple local minima, rather than important relaxation of the surroundings.
For many of the E n states, we have observed that a transfer of the hydroxy proton of homocitrate to the N 2 -intermediate is favourable. This indicates homocitrate may constitute a proton buffer, which can be used to stabilise certain intermediates of the reaction, especially H 2 NNH 2 and NH 3 . This provides an attractive explanation why homocitrate is mandatory for the function of nitrogenase. [88,99] Bjornsson and coworkers have also suggested that the hydroxy proton of homocitrate is employed to form trans-HNNH from N 2 in the E 4 intermediate, although the reaction was uphill. [42] For all E n levels, binding of the intermediates to Fe6 seems to be preferred and therefore we have suggested such binding in the mechanism in Figure 8. This site has more hydrogenbond possibilities (besides sulfide ions), involving His-195, Gln-191 and the homocitrate ligand. It can also employ the suggested proton buffer of homocitrate and it is also closer to the end of the suggested proton-transfer path, involving a chain of water molecules, ending close to homocitrate and S3B. [94][95][96][97] However, for the E 4 and E 8 states, binding of the N 2intermediates to Fe2 are competitive. Moreover, it has been suggested that a likely N 2 -binding channel ends at Fe2. [27] We have previously studied the reaction mechanism of nitrogenase, assuming that S2B dissociates from the FeMo cluster, opening up for an obvious binding site of the substrate. [57] This gave rise to a mainly alternating reaction mechanism, in which the substrate and the intermediates bound in a bridging mode (with one or both N atoms) between Fe2 and Fe6. In the present study, bridging intermediates were always found to be much less stable than end-on intermediates, except for E 7 structures with both NH 2 and NH 3 . The reason for this is that the FeMo cluster is too crowded if both the substrate and S2B bridges Fe2 and Fe6. In fact, it is often observed that S2B moves to other positions or reacts with the substrate or other sulfide ions in such high-energy structures with a bridging substrate. Thus, it seems clear that bridging substrate structures are unlikely when S2B remains bound.
In conclusion, this study shows that the second part of the nitrogenase reaction (after binding of N 2 ) is possible also if the S2B ligand has not dissociated. Such a reaction follows an alternating mechanism with the substrate and intermediates binding to Fe6. It has also pointed out an important role for the homocitrate ligand as a proton buffer. In future studies, we study the binding of N 2 and dissociation of H 2 to the E 4 state of the cluster.