Understanding the Electronic Structure Basis for N2 Binding to FeMoco: A Systematic Quantum Mechanics/Molecular Mechanics Investigation

The FeMo cofactor (FeMoco) of Mo nitrogenase is responsible for reducing dinitrogen to ammonia, but it requires the addition of 3–4 e–/H+ pairs before N2 even binds. A binding site at the Fe2/Fe3/Fe6/Fe7 face of the cofactor has long been suggested based on mutation studies, with Fe2 or Fe6 nowadays being primarily discussed as possibilities. However, the nature of N2 binding to the cofactor is enigmatic as the metal ions are coordinatively saturated in the resting state with no obvious binding site. Furthermore, the cofactor consists of high-spin Fe(II)/Fe(III) ions (antiferromagnetically coupled but also mixed-valence delocalized), which are not known to bind N2. This suggests that an Fe binding site with a different molecular and electronic structure than the resting state must be responsible for the experimentally known N2 binding in the E4 state of FeMoco. We have systematically studied N2 binding to Fe2 and Fe6 sites of FeMoco at the broken-symmetry QM/MM level as a function of the redox-, spin-, and protonation state of the cofactor. The local and global electronic structure changes to the cofactor taking place during redox events, protonation, Fe–S cleavage, hydride formation, and N2 coordination are systematically analyzed. Localized orbital and quasi-restricted orbital analysis via diamagnetic substitution is used to get insights into the local single Fe ion electronic structure in various states of FeMoco. A few factors emerge as essential to N2 binding in the calculations: spin-pairing of dxz and dyz orbitals of the N2-binding Fe ion, a coordination change at the N2-binding Fe ion aided by a hemilabile protonated sulfur, and finally hydride ligation. The results show that N2 binding to E0, E1, and E2 models is generally unfavorable, likely due to the high-energy cost of stabilizing the necessary spin-paired electronic structure of the N2-binding Fe ion in a ligand environment that clearly favors high-spin states. The results for models of E4, however, suggest a feasible model for why N2 binding occurs experimentally in the E4 state. E4 models with two bridging hydrides between Fe2 and Fe6 show much more favorable N2 binding than other models. When two hydrides coordinate to the same Fe ion, an increased ligand-field splitting due to octahedral coordination at either Fe2 or Fe6 is found. This altered ligand field makes it easier for the Fe ion to acquire a spin-paired electronic structure with doubly occupied dxz and dyz orbitals that backbond to a terminal N2 ligand. Within this model for N2 binding, both Fe2 and Fe6 emerge as possible binding site scenarios. Due to steric effects involving the His195 residue, affecting both the N2 ligand and the terminal SH– group, distinguishing between Fe2 and Fe6 is difficult; furthermore, the binding depends on the protonation state of His195.


INTRODUCTION
Nitrogenase is the only type of enzyme capable of reducing dinitrogen gas to ammonia at ambient temperature and pressure and has been studied for decades. 1−5 The reaction has been proposed to follow the following stoichiometry: where in addition to the 2 NH 3 molecules, a molecule of H 2 is formed. This obligatory H 2 formation is now confidently known to be related to a critical H 2 elimination step involving hydrides during the reaction. Extensive kinetic studies 6−8 by Lowe and Thorneley (LT) led to a pioneering kinetic scheme that is still used to discuss the mechanism of the reaction. A simplified LT cycle is presented in Figure 1, where reduced intermediates are named E n , where the subscript n refers to the additional electrons/protons relative to the resting state (E 0 ). A recent study 9 discusses a steady-state kinetic model that establishes even more clearly the rate constants of several critical steps: the rate of H 2 revolution from the E 2 state was, e.g., found to be lower than the more reduced E 4 state, while H 2 formation from E 4 and N 2 binding to the E 4 state have equal rate constants. A previously unidentified slow process was also found, believed to be related to an N 2 activation process.
The FeMo cofactor (FeMoco) of the more active Modependent nitrogenase is generally accepted as the site where N 2 conversion to NH 3 takes place. 7,16,17 FeMoco contains eight open-shell metal ions (one Mo and seven Fe ions), with a complex electronic structure featuring antiferromagnetic coupling of high-spin Fe 2+/3+ ions, mixed-valence delocalization, and a strange non-Hund electron configuration at the Mo 3+ site. 18 The resting E 0 state has a spin of S = 3/2, and a charge of −1 (of the [MoFe 7 S 9 C] part) has been proposed as the most likely redox state based on calculations 10, 19 and experiments. 20 Together with the previous Mo 3+ assignment, 18 this suggests a combined Mo 3+ 3Fe 2+ 4Fe 3+ oxidation state, although calculations indicate that the electronic structure can be rather delocalized. 10, 19 Broken-symmetry density functional theory remains the primary computational methodology utilized to study the electronic structure and reactivity of the system. Ten brokensymmetry (BS) solutions of FeMoco were originally suggested based on approximate C 3 symmetry of the cofactor and where only Fe ions were considered as part of the magnetic coupling. 21−23 These solutions become 35 if symmetry is ignored. The most favorable BS7 class of solutions consist of the BS7−235, BS7−247, and BS7−346 solutions (abbreviated as BS235, BS247, and BS346 from now on) where the three numbers indicate which Fe ions are spin-down (X-ray structure numbering). These three solutions cannot be distinguished energetically (being within ∼1−2 kcal/mol of each other) but they do lead to distinct geometric Fe−Fe/Mo−Fe distances, with the BS235 solution giving the best agreement with the high-resolution X-ray structure. 10 Localized orbital analysis 10 of the cofactor indicates that these three BS7 solutions differ in their locations of mixed-valence pairs and localized Fe 2+ or Fe 3+ ions, which explains the specific geometric effects behind these almost isoenergetic solutions.
The integer-spin E 1 state of FeMoco has an additional e − / H + pair and is likely S = 1 or 2. 24 A recent Mo,Fe XAS and Mossbauer study 25 suggests that Fe-based reduction has occurred, while a joint Mo,Fe-EXAFS-QM/MM study 11 further suggests that the added proton is likely on either the bridging belt sulfide S2B or S5A. The calculations showed a preference for a M S = 2 spin state and the Fe reduction was found to be localized in the MoFe 3 subcubane. The S2B sulfide, bridging Fe2 and Fe6 (shown in Figure 1), appears most likely as the protonation site in E 1 , as X-ray crystallography 26−28 has revealed various structures indicating the lability of this sulfide position and computational studies have also suggested protonation of S2B in the E 1 state. 29,30 The E 2 state is known to possess an S = 3/2 spin state with two conformers present with distinct g-tensors (g = [4.21, 3.76, ∼1.97] and [4.69, ∼3.20, ∼2]) and a hydride likely present in at least one of them. 31−34 In a recent study 12 from our group, we proposed two models for the E 2 state: E 2 -hyd and E 2nonhyd. One contains a bridging hydride between Fe2 and Fe6 as well as a terminal sulfhydryl group (labeled E 2 -hyd here). The E 2 -hyd with a terminal SH − on Fe6, named as E 2hyd-SH − @Fe6, is more stable than the one with a SH − on Fe2, as also seen in a recent QM/MM study. 35 The E 2 -nonhyd model contains a doubly reduced Fe environment and two protons located on S2B and S5A. Both models are shown in Figure 1. As discussed both in ref 12 and more recently in ref 35, E 2 isomers are, however, quite sensitive to the DFT method used with TPSS, e.g., giving opposite trends to r 2 SCAN and TPSSh. A saddlepoint leading to H 2 formation from E 2 -hyd → E 0 was also reported in ref 12, with a fairly high activation barrier of ∼20 kcal/mol, which is consistent with a recent kinetics study, 9 indicating the E 2 → E 0 step to be slower than other H 2 formation pathways.
Little is unfortunately known experimentally about the E 3 state (with even few computational studies 29,30 dedicated to it) due to its EPR-silent nature and it will not be discussed further here. In contrast to the EPR-silent E 3 , the EPR-active E 4 state is much better characterized. Extensive EPR studies have revealed a spin-state of S = 1/2 in both the α-70 Val → lle mutant and native enzyme. 36,37 ENDOR studies have revealed that the state contains two bridging hydrides, bridging between the Fe ions (Fe-H-Fe) with two other protons bound to sulfides. 36,38 A complete structural model for E 4 remains controversial, however. The first proposed models based on the ENDOR data suggested hydrides bridging Fe ions (and protonated belt sulfides) without any associated structural changes suggested. 39 Computational studies have suggested models with either bridging or terminal hydrides as energetically favorable states. 13,15,30,38 Some of these models feature hydrides bridging different pairs (DP) of Fe ions (Fe2 and Fe6, Fe3 and Fe7) with closed belt sulfide bridges and were shown to be in good agreement with 1 H hyperfine tensor orientations, 15,38 here labeled E 4 -DP-Fe2/6(5) and E 4 -DP-Fe2/6(3), where Fe2/ 6(5) and Fe2/6(3) labels indicate that the hydrides bridging Fe2 and Fe6 points to either S5A and S3A, respectively. E 4 -DP-Fe2/6(3) was calculated to be lower in energy than E 4 -DP-Fe2/6(5), but only the latter agreed with ENDOR results. 15 A QM/MM study (using the TPSSh functional) from our research group in contrast suggested E 4 models with two bridging hydrides located on the same pair (SP) of Fe ions (Fe2 and Fe6) and with open belt sulfur bridges (terminal sulfhydryl groups) that were found to be energetically more favorable than models with closed belt sulfur bridges (including E 4 -DP-Fe2/6(5)). Such opening of protonated belt sulfides has also been discussed in computational work by Dance, 40 Raugei,13 and recently by Ryde and co-workers. 35 Kinetic and spectroscopic studies 7,37,41−44 indicate that dinitrogen will bind to the E 4 state, while no binding has been directly observed in earlier E n states (with the role of E 3 unclear). Alternative substrates and inhibitors like hydrazine, diazene, acetylene, and CO likely bind on the other hand to the E 1 or E 2 states. 4,5 The Fe2-Fe3-Fe6-Fe7 face has been proposed as a likely site for N 2 binding based on mutation studies 45,46 of α-70 Val where the N 2 reduction is dramatically decreased, with the Fe2 and Fe6 sites being the primary possibilities. Unfortunately, experimental characterization of dinitrogen-bound intermediates is scarce and the nature of the molecular and electronic structure of the cofactor prior to and after dinitrogen binding is unclear. Additionally, crystallographic studies 26−28 have shown that the bridging sulfide (S2B) between Fe2 and Fe6 can be displaced under certain conditions, implicating Fe2 and Fe6 as a site of ligand binding, but which could also be interpreted as sulfide lability being mechanistically relevant. 47,48 A recent study presenting an Xray crystal structure 28 apparently showing a missing sulfide and bound N 2 is controversial. 49−51 Some recent computational studies have investigated the protonation and dissociation of S2B and the mechanism of N 2 reduction after S2B dissociation, but it remains unclear whether S2B loss is part of the catalytic mechanism or not. 52−56 Other computational studies 13,57−60 have investigated N 2 binding to the whole cofactor; however, N 2 binding is rarely found to bind favorably to FeMoco. In a previous QM/MM study 14 from our group, however, we found that E 4 -SP type cofactor models will favorably bind N 2 .
H 2 evolution has for a long time been known to be a compulsory mechanistic step when reducing dinitrogen. 7 Unlike the unproductive side-reactions E 2 → E 0 and E 4 → E 2 (that proceed via a hydride-protonation, hp, mechanism 37,38 ) the obligatory H 2 formation (the H 2 formation part of eq 1) is known to occur via another mechanism, typically described as reductive elimination, re (though this nomenclature has also been challenged 61 ). This elimination involves a direct reaction of the two bridging hydrides as N 2 binds. 62 N 2 thus appears to act as a catalyst for the specific H 2 formation via the re mechanism (as established via isotope labeling 37,38 ). Hoffman and co-workers have classified the mechanistic possibilities for this N 2 /H 2 reaction as "dissociative", "concerted", "associative", and "intermediate". 63 A "dissociative" mechanism (i.e., H 2 leaving first) can be ruled out based on the isotope labeling experiments, while "concerted" (simultaneous H 2 dissociation and N 2 binding), "associative" (N 2 binding followed by H 2 dissociation), and "intermediate" (N 2 binds to a higher energy intermediate E 4 followed by H 2 leaving) remain possible scenarios. In the "concerted" scenario, H 2 evolution simply pays the energetic cost associated with N 2 binding. However, it seems hard to imagine such a selective mechanistic step taking place without N 2 binding occurring first. In the "associative" scenario, N 2 binding causes the elimination of H 2 , presumably leading to the formation of a more stable N 2 -bound state. This hypothesis, however, implies that N 2 binding to the cofactor has somehow become more favorable than before. Finally, in the "intermediate" scenario, H 2 formation occurs first, but with H 2 remaining bound in a higher energy intermediate, and next gets displaced by N 2 . A QM-cluster study 13 and a combination 64 of hybrid QM/MM metadynamics simulations and ENDOR measurements have suggested such a mechanism.
It is this initial binding of N 2 (primarily within the "associative" scenario defined above) that is the focus of our study. We will focus on the uncovering of both the electronic and molecular structure aspects of FeMoco that enable N 2 binding to occur. In our goal to understand why the E 4 state is capable of binding dinitrogen it is equally useful to understand why the earlier E 0 , E 1 , E 2 states are not capable of binding dinitrogen and to understand how both the cofactor redox state as well as single Fe oxidation state and local spin states correlate with binding energies. Diamagnetic substitution of the cofactor with a quasi-restricted orbital (QRO) transformation allows us to gain such insights, and we utilize this technique to shed light on the local ligand field of the individual Fe ions involved in N 2 binding in each of the redox, spin, and BS states calculated in this work. Finally, we discuss how the protonation state of His195 affects binding energies.

COMPUTATIONAL DETAILS
The QM/MM model in this work builds on a previous model setup from our group 10 that has successfully been used for several E n states. 10−12,14 However, our previous studies employed a truncated spherical droplet model for the actual QM/MM modeling. While the spherical droplet QM/MM model is known to be a highly economical yet accurate strategy for QM/MM reaction profile studies on proteins, 65,66 a recent study from our group on VFe protein of V nitrogenase 67 demonstrated some sensitivity to the size of the model, suggesting that bulk electrostatics have an effect on the redox properties of the cofactor. In this work, we utilize the full box from the original classical MM setup of MoFe protein 10 without truncating the system.
A new QM/MM code, ASH, 68 developed in our group, interfaced to the OpenMM library 69 and the ORCA quantum chemistry code 70 was used in this study. ASH has been designed with the aim of having full flexibility to perform QM/ MM calculations on large metalloproteins, with an interface to the OpenMM molecular mechanics library and a flexible interface to the ORCA quantum chemistry code with specific functionality to conveniently control BS-state convergence. The CHARMM36 force field 71 was used to describe the MM region as in previous papers in our group. 10−12,14,67,72−74 Geometry optimizations were performed using an interface to the geomeTRIC optimization library 75 using HDLC internal coordinates. 76 Standard electrostatic embedding QM/MM with link atoms and a charge-shifting scheme was used. 77 A comparison between ASH and Chemshell 78 results was performed on the spherical droplet model to demonstrate that previous QM/MM results were reproducible (see Section 17 in the SI). The full solvated model of MoFe protein has 320,814 atoms in the E 0 state (see Figure S1a in the SI). An active region of 990 atoms was defined, centered on FeMoco (see Figure S1b in the SI).
In this work, geometry optimizations were first obtained for all tested spin states and BS solutions using the smaller QM-I region and the larger QM-V region was used to perform geometry optimizations on the lowest-energy isomers (according to the QM-I data). Finally, single-point calculations using the large QM-VI region were performed on the QM-V geometries. Section 2 in the SI compares energy differences obtained with different QM regions.
2.1. Quantum Chemical Details. The r 2 SCAN 81 density functional was used in QM and QM/MM calculations using the implementation in the LibXC library. 82 This functional emerged as the most accurate in a recent benchmarking study from our group on the structural properties of spin-coupled iron−sulfur dimers as well as FeMoco (compared to highresolution X-ray structures), 74 giving on average slightly better molecular structures than the TPSSh functional that we have previously favored in our FeMoco research. In particular, the Fe−Fe and Mo−Fe distances of the spin-coupled metal dimers were found to be highly sensitive to the treatment of the covalency of the Fe−S bond, suggesting the r 2 SCAN functional to better treat the electronic structure of complex Fe−S clusters. The r 2 SCAN functional has furthermore been shown to be a highly reliable density functional for both maingroup as well as closed-shell and open-shell transition metal reaction energies (primarily organometallic systems tested) and transition metal spin-crossover energies according to benchmark studies. 83−85 The r 2 SCAN functional has additionally the practical advantage of not including HF exchange, unlike the TPSSh functional. This results in less costly QM steps, which allowed us to explore more possibilities in this work as well as using large QM-regions while making no sacrifice with respect to overall electronic structure treatment of FeMoco. SCAN, 86 the predecessor to r 2 SCAN, was known to suffer from numerical problems, which led to the development of r 2 SCAN, which suffers less from numerical problems. Nonetheless, we investigated the numerical grid dependence of the r 2 SCAN exchange−correlation integrals, as shown in the Section 20 of the SI. The "defgrid2" integration grid in ORCA, used in most calculations in this work, was found to give satisfactory numerical precision; we estimate numerical grid uncertainty in reaction energies to be on the order of ∼0.2 kcal/mol.
In this work, we continue to utilize an all-electron scalar relativistic treatment that we have shown to be close to the allelectron relativistic basis set limit of iron−sulfur compounds. 74 The D4 dispersion correction 84,87 was used in all calculations. The Split-RI-J approximation 88 as implemented in ORCA was used to speed up the Coulomb integral evaluations utilizing a decontracted general Coulomb fitting auxiliary basis set 89,90 (named SARC/J in ORCA). The ZORA scalar relativistic Hamiltonian 91,92 was used with the relativistically recontracted ZORA-def2-TZVP 89,93 basis set on Fe, S, carbide, hydrides, and S-bound protons of FeMoco as well as the dinitrogen ligand in N 2 -bound states. The all-electron SARC-ZORA-TZVP 94 basis set was used for Mo. The smaller ZORA-def2-SVP basis set was used for other atoms (including protein residues). Test calculations indicated that differences in using ZORA-def2-SVP or ZORA-def2-TZVP on the homocitrate were tiny (see Section 2 in the SI). The smaller basis set (ZORA-def2-SVP) on homocitrate was used for optimizations with the QM-I region but was increased to ZORA-def2-TZVP when doing single-point calculations using the QM-VI region. Localized orbital calculations used the Pipek−Mezey localization method 95 and was used to derive electron configuration diagrams of the most relevant states of FeMoco (shown in the SI). The CPCM solvation model 96 97 No additional non-electrostatic term was included.
Hirshfeld population analysis 98 was used to calculate atom charges and spin populations. Vibrational frequencies were calculated numerically using a partial Hessian approach 99 (N 2 and the binding Fe defining the partial Hessian). The calculated vibrational frequencies of N 2 were not scaled. The stretching vibration of free N 2 at the r 2 SCAN level is 2432 cm −1 , which is an overestimate compared to the experimental gas value (2330 cm −1 ). In this study, we focus on the shift (Δ) in N 2 frequency relative to free N 2 at the same level. The unbound N 2 was calculated in vacuum for all QM/MM binding energy calculations, while in CPCM-cluster binding energy calculations, CPCM was included (ε = 4).
Binding energy in this work will be used to refer to the energy of binding N 2 to the most favorable isomer of each redox state while we will occasionally use the term "single-step N 2 binding energy" when referring to N 2 binding to a specific isomer. Finally, we note that reported N 2 binding energies in this work are based on electronic energies at 0 K: ΔE 0K (without zero-point vibrational energy, ZPVE). The focus of this study is not to give realistic estimates of free energies of binding of N 2 but rather to understand the trends in the electronic part of the binding energy. While we expect ΔE 0K to be a decent approximation to ΔH 298K (ZPVE being the largest missing contribution), we do not include entropic contributions, which will have considerable contributions to the free energies of binding. Unfortunately, translational and rotational entropic contributions from gas phase statistical mechanics (10−15 kcal/mol in favor of the reactant in an A + B → C reaction) may also overestimate entropy effects; in this work, we will report ΔE 0K binding energies as this is sufficient to discuss the relevant trends, but it should be kept in mind that free energies of binding will differ.

Diamagnetic Substitution and Ligand-Field Diagrams.
Due to the complex delocalized electronic structure of Fe and Mo in FeMoco, understanding the local electronic structure as N 2 binds to Fe is challenging. Localized orbital analysis has been used in previous QM/MM studies from our group 10,12,14,67,73 as it allows us to derive approximate electron configurations diagrams of all Fe and Mo ions in each FeMoco state, something not possible by decomposition of the spin density. Unfortunately, localized orbitals prevent us from discussing a ligand-field model for each local Fe ion as localized orbitals have no well-defined energies. Diamagnetic substitution (DS) has in previous studies proved useful to get insights into the local electronic structure of FeMoco. 18,73 Here, all open-shell metal ions, except a chosen one, are replaced by diamagnetic ions with a similar atom radius.
Diamagnetic substitution was used in this work in two ways: (i) to study N 2 binding in a diamagnetically substituted version of the resting state cofactor, varying the Fe oxidation state and spin state of a single Fe ion (Fe6). During optimization, only Fe6 and N 2 were allowed to move. (ii) to derive approximate ligand-field diagrams of specific Fe sites of FeMoco intermediates throughout the study by utilizing quasi-restricted orbitals of a diamagnetically substituted model. In each diamagnetically substituted calculation, Mo 3+ was substituted with In 3+ and Fe 3+/2+ centers were substituted with Ga 3+ . The Ga 3+ ion nuclear charges were modified to have a value of 30.5 rather than 31 to make the atomic charges of the sulfides closer to that of the full FeMoco in the E 0 state as previously used. 73 Because a diamagnetically substituted model contains only a single open-shell Fe ion, we can derive approximate ligand-field diagrams based on the quasi-restricted orbital transformation of the canonical orbitals. Ligand-field diagrams were obtained in a three-step procedure: (i) determine the local redox state and spin state of the Fe center of interest according to localized orbital analysis of the full FeMoco using the QM-VI region. (ii) diamagnetically substitute all other Fe ions except one ion (kept in the redox state and spin state suggested by the first step). (iii) perform a single-point calculation and derive quasirestricted orbitals, which were used to make a ligandfield diagram. A step-by-step example is shown in Section 15 of the SI. The ligand-field diagrams derived in this way are of course approximate due to both (i) the spin-pairing enforced by the quasi-restricted orbital transformation of an unrestricted Kohn−Sham determinant with different orbitals for different spins and (ii) the diamagnetic substitution changing the chemical system. For the case of some high-spin metal ion cases, non-Aufbau configurations are found, which is an artifact of the approach. Future work will consider an ab initio ligandfield approach 100 instead.

Broken-Symmetry
States. Ten classes of BS determinants were originally proposed to describe FeMoco (in its E 0 state) assuming C 3 symmetry. 21,22 If the approximate C 3 symmetry is ignored, however, this results in 35 ) BS solutions (flipping all 7 Fe sites). This is now known to be a simplification as the Mo ion is an open-shell ion as well. 18 However, as the Mo ion preferentially adopts an unusual spin-coupled non-Hund configuration, 18 it is not necessary to consider Mo as part of the spin-flipping problem. The BS determinants can be described in terms of which Fe ions have predominantly local spin-down vectors (according to spin population analysis of the SCF solution). The most favorable broken-symmetry solutions for the E 0 state are of the BS7 class and are here labeled as follows: BS235, BS247, and BS346. While energetically very similar (∼1 kcal/mol), these BS determinants differ electronically in terms of locations of mixed-valence delocalized Fe pairs and localized Fe 2+ or Fe 3+ ions. 10 The localized orbital analysis of these three BS solutions is shown in Figure S5 of the SI. The BS235 solution has Fe2, Fe3, and Fe5 spin-down, which results in mixedvalence delocalized pairs between Fe2 and Fe3 (being both spin-down), between Fe1 and Fe4 (being both spin-up), and between Fe6 and Fe7 (being both spin-up), as shown in Figures S5 and S6. For reduced states of FeMoco, especially when the hydride or N 2 is bound, the nature of the spin coupling problem changes, apparently due to the stability of intermediate or low-spin local configurations on H/N 2 -binding Fe ions. Thus, the BS147 (of BS10 class) was, e.g., found to be the most stable for E 4 -SP. 14 In this study, we have used four BS solutions (BS235, BS247, BS346, and BS147) in initial explorations of all models of the E n states. Additionally, we present additional calculations on more BS states for different E 4 isomers that contain a more complicated electronic structure (including Fe ions with low or intermediate spin states) than the other E n states in Section 11 of the SI.

RESULTS
N 2 binding was systematically explored to Fe2 or Fe6 sites in QM/MM calculations of the E 0 , E 1 , E 2 , and E 4 states and is discussed in Sections 3.2, 3.3, 3.4, and 3.5, respectively. Section 3.6 shows the results of N 2 binding to E 4 for different protonation states of α-His195 on N 2 binding to the E 4 -SP model. We start our discussion about N 2 binding to FeMoco, however, by showing the results for a diamagnetically substituted version of FeMoco with Fe6 in multiple oxidation and spin states. In our study, we primarily focused on end-on N 2 binding to Fe2/Fe6 sites trans to the Fe−C bond. As discussed in the SI, other N 2 binding modes are generally found to be unfavorable. Our study is also focused on N 2 binding energies (calculated as 0 K electronic energies) and are not intended to be realistic estimates of the free energies of binding at room temperature (entropic effects would shift the binding energies considerably).

N 2 Binding to Diamagnetically Substituted FeMoco: Fe(III) vs Fe(II) vs Fe(I). Substituting all Fe and
Mo ions for diamagnetic ions (Ga 3+ and In 3+ ) in FeMoco except Fe6 in a cluster model of FeMoco allows us to explore how N 2 binding to a local Fe ion in the FeMoco coordination environment is affected by local oxidation states and spin states without the complications arising from the complex spin coupling and mixed-valence delocalization present. The diamagnetically substituted QM-cluster model, shown in Figure S1, was calculated utilizing a CPCM continuum solvation model. The following eight possibilities for Fe6 were studied: Fe 3+ (S = 5/2, 3/2, 1/2), Fe 2+ (S = 2, 1, 0), and Fe 1+ (S = 3/2, 1/2). We note that the belt Fe ion positions are essentially equivalent in a diamagnetically substituted model of E 0 with no protein environment present.
Prior to N 2 binding, the lowest-energy isomer for each redox state has the Fe ion in a high-spin state. Dinitrogen was found to bind to five of the eight states calculated: Fe 3+ (S = 1/2), Fe 2+ (S = 1, 0), and Fe 1+ (S = 3/2, 1/2). N 2 otherwise spontaneously dissociated for the other cases. The low spin states found, imply that N 2 -bound states of FeMoco will favor Inorganic Chemistry pubs.acs.org/IC Article a change in spin state upon binding. The QRO ligand-field diagrams of these five states (see Figure 2c and Figure S3 in the SI) reveal the common theme of d xz and d yz orbitals being doubly occupied (where the z-axis is defined along the Fe−N 2 bond). The d xz and d yz orbitals have the correct symmetry to overlap with the π* orbitals of dinitrogen, and considerable N 2 π* orbital character can in fact be seen in the QRO orbitals (see Figure 2c for the Fe(II) S = 1 case). Further spin-pairing by double occupation of the d xy orbital as in a Fe 2+ S = 0 state offers no additional advantage, in fact the binding energy becomes less favorable (see Figure 2a). As N 2 binds to the cofactor, an approximate trigonal bipyramidal geometry might be expected; however, such a geometry is in fact not stable as the Fe−C distance is strongly elongated from 2.08 to 2.72 Å (see Figure 2), suggesting a completely broken Fe−C bond. This Fe−C "bond cleavage" occurred for all five N 2 -bound isomers (see Figure S3). Clearly, this suggests that the coordination change of the binding site Fe is likely to accompany N 2 binding and we note that a flexible Fe−C interaction has indeed been suggested as a role for the interstitial carbide, with model compounds from Peters being particularly relevant in this context. 101 However, recent 13 C ENDOR data of various states of FeMoco are more indicative of a stabilizing rather than flexible role for the interstitial carbide. 102 In the diamagnetically substituted model, only the Fe and N 2 atoms were allowed to move and with no protons present either, and elongation of the Fe−C bond constitutes the only structural rearrangement possible.
The N 2 binding energy becomes more favorable in going from the Fe 3+ state, via Fe 2+ , to the Fe 1+ state. While for the Fe 2+ state, the N 2 binding is still endothermic by +13.9 to +18.9 kcal/mol, for the Fe 1+ state, the N 2 binding has become convincingly exothermic by −15.1 to −30.8 kcal/mol. A more reduced Fe obviously makes it easier to generate double occupations of d xz and d yz orbitals (that overlap with N 2 π*), and a more reduced Fe can even more effectively push electron density into the N 2 π* orbitals as the Fe d-orbital energies should be shifted to higher energy (closer to N 2 π* levels). The vibrational frequency of the N−N bond (a well-known indicator of N 2 bond activation 103 ) relative to free N 2 , ν NN , correlates well with the trend in binding energies. ν NN decreases from −217 to −476 cm −1 as the Fe oxidation state goes from Fe 3+ to Fe 1+ , which indicates unsurprisingly that a more reduced Fe is positively correlated with the activation degree of N 2 .
Overall, the results using the diamagnetically substituted FeMoco model demonstrate that double occupations of d xz and d yz orbitals, a reduced Fe, and the coordination change of the Fe binding site are relevant to favorable terminal N 2 binding. Such local spin-state changes of the multi-Fe spincoupled FeMoco seem likely to occur as N 2 binds to the cofactor; however, structural changes (including those associated with protonation) are likely to play an equally important role, which are not realistically captured in this simple model.  states. Calculated ν NN frequencies are relative to free N 2 (2432 cm −1 at the r 2 SCAN level of theory). Each binding energy is relative to the lowestenergy isomer (always high spin) for each redox state without N 2 . Crosses in red indicate that N 2 does not bind to Fe of that state. (b) Lowestenergy structure of DS FeMoco with a high-spin Fe 2+ at Fe6 before N 2 binding. (c) Optimized structure of N 2 bound to DS FeMoco with a Fe 2+ (S = 1) and the associated ligand-field diagram of the Fe 2+ based on quasi-restricted orbitals. Also shown are d xz and d yz quasi-restricted orbitals. The full cluster model is shown in Figure S1 of the SI.  Figure S5 in the SI shows the electronic structure diagrams of all three BS solutions based on the localized orbital analysis. Sixteen N 2 -bound states were studied by binding N 2 to Fe2 or Fe6 in the E 0 state in two spin states (M S = 3/2, 1/2) and four BS solutions using the QM-I region, as shown in Figure 3, but only three of these states were found to give N 2 -bound minima: E 0 -N 2 @Fe6-BS147, E 0 -N 2 @Fe2-BS147, and E 0 -N 2 @ Fe2-BS346, all of which are in M S = 1/2 spin state. The N 2 binding energies to Fe6 and Fe2 are found to be 16.6 and 21.4 kcal/mol, respectively, using the larger QM-VI region and the BS147 solution. The difference in N 2 binding energy between Fe2 and Fe6 sites appears to be due to the steric clashing of the N 2 ligand with the His195 residue when bound to Fe2. This preference is reversed when a cofactor-only cluster-continuum model is used, as shown in Figure S7 in the SI.
Localized orbital analysis (see, e.g., Figure S29) and Hirshfeld spin populations are consistent with local spin state changes occurring for most of the N 2 -bound states. For the N 2 @Fe6 state highlighted in Figure 3c, the low Hirshfeld spin population of 1.53 on Fe6 is highly suggestive of a local spin state change and localized orbital analysis confirms spinpairing occurring for d xz and d yz orbitals, as seen in Figure S29 of the SI. The QRO-based ligand-field diagram of the Fe6 ion (using diamagnetic substitution and choosing Fe 2+ S = 1 based on localized orbital analysis) reveals a changing ligand field where the d z 2 orbital is destabilized due to the approximate trigonal bipyramidal geometry in the N 2 -bound state.
As shown in the FeMoco structures before and after N 2 binding to Fe6, the Fe6-C and Fe6-S2B distances become longer after N 2 binding (Figure 3), although the Fe6-C elongation is far from the complete cleavage that was seen for the diamagnetically substituted model in Section 3.1. Clearly, however, this implies that N 2 binding to Fe will lead to structural changes that are likely to depend on the protonation state of the cofactor.
As higher E n states of the cofactor are known to involve protonated belt sulfides (the S2B sulfide bridging Fe2 and Fe6 in particular), we also carried out E 0 calculations with an additional H + present on S2B, labeled E 0 H + here (but without an electron as in the E 1 state). We note in this context that Rees and co-workers have suggested based on X-ray structures (1.85 and 2.30 Å resolution) of MoFe proteins (both Azotobacter vinelandii and Clostridium pasteurianum) under pH = 4.5−6 conditions, that an E 0 H + state may rather be protonated at either S5A or S3A. 105 As our study is restricted   Figure 4 shows the binding of N 2 to the FeMoco E 0 H + model. The results reveal that the availability of an additional proton on S2B introduces additional flexibility for N 2 binding, the proton enabling spontaneous full belt sulfide-bridge opening to give a partially or even fully terminal sulfhydryl group, and additional favorable N 2 -bound states could be found for both N 2 @Fe2 and N 2 @Fe6 (compared to regular E 0 ). This hemilability of protonated belt sulfur bridges in reduced E n states has been found to be a common structural theme in previous reduced and ligand-bound studies by us 12,14,73 but had not previously been explored at the E 0 redox level. However, as for the E 0 state, the N 2 binding energy of N 2 remains endothermic and high (>16.8 kcal/mol), suggesting structural rearrangement alone via sulfur protonation not to be sufficient for favorable N 2 binding to FeMoco. This lack of N 2 binding to the resting state is of course fully consistent with experimental observations. A local spin-state change on Fe6   Figure S9 of the SI suggests that Fe6 is best described as a high-spin Fe 2+ ion. The results shown in (b) and (c) were calculated using the QM-VI region. The N 2 binding energies are relative to the E 1 -BS346 with M S = 2 (b) when using both QM-I and QM-VI regions.

Inorganic Chemistry
pubs.acs.org/IC Article was found for the E 0 H + -N 2 @Fe6 state but not for the E 0 H + -N 2 @Fe2 state (see Figure 4).

N 2 Binding to FeMoco: The E 1 Redox State.
While there is no indication of the E 1 state binding N 2 , this state has been implicated in binding alternative substrates and inhibitors such as CO and acetylene, 5,73 indicating that it possesses increased reactivity compared to the resting state. Previous computational studies 29,30 and a recent joint EXAFS-QM/MM study 11 suggest the state to be best described as an Fe-reduced cofactor with a protonated S2B belt sulfide bridge. The spin state is known to be integer spin, 24 and an M S = 2 brokensymmetry state was the most favorable according to the QM/ MM calculations. 11 The added electron was found to localize in the MoFe 3 subcubane, on either Fe5, Fe6, or Fe7, depending on whether the BS235, BS346, or BS247, respectively, were calculated, with the spin-down Fe in the MoFe 3 sub-cubane, thus determining the site of the redox event (Fe 3+ → Fe 2+ ).
Twenty-four isomers were systematically explored with N 2 bound to Fe2 or Fe6 in three spin states (M S = 0, 1, 2) and four BS solutions. As shown in Figure 5, considerably more states were found to give N 2 -bound local minima at the E 1 redox level compared to the E 0 redox level. As found for E 0 H + , all N 2 -bound states resulted in the protonated belt-sulfide (S2B) becoming terminal on the other Fe as N 2 bound to Fe2 or Fe6 . The hemilability of the belt-sulfide thus acts to preserve the (distorted) tetrahedral coordination of both Fe ions as N 2 binds.
Localized orbital analysis was performed for all N 2 -bound states, which showed that most states showed spin pairing of d xz and d yz orbitals at the N 2 -binding Fe site, with some showing spin-pairing of only one of these orbitals. For the most favorable N 2 -bound state, E 1 -N 2 @Fe6, shown in Figure 5c, however, only the d xz orbital on Fe6 showed spin-pairing, as did the state without N 2 according to similar Hirshfeld spin populations of Fe6 (see ligand-field diagrams in Figure 5c). The slightly more favorable N 2 binding energies in the E 1 state for E 1 -N 2 @Fe6 and E 1 -N 2 @Fe2 states of 9.9 and 11.2 kcal/ mol, respectively, compared to the E 0 /E 0 H + states are likely to be primarily due to the overall more reduced Fe environment of the cofactor. This also results in increased activation of the N 2 ligand with a ν NN = −287 cm −1 shift of the N−N stretching frequency compared to free N 2 . The N 2 binding energy is still found to be rather endothermic, which is fully consistent with the lack of experimental evidence for N 2 binding in the E 1 state. This result can be contrasted with the more favorable calculated binding energy for CO in this same redox state, which was found to be −8.3 kcal/mol (TPSSh-QM/MM level of theory) in our previous study. 73 CO binding to Fe6 of the E 1 state using the same QM level as N 2 binding to FeMoco Inorganic Chemistry pubs.acs.org/IC Article (r 2 SCAN and the QM-I region) was also studied, shown in Section 7 of the SI, indicating a similar binding energy (−9.6 kcal/mol) as our previous study. Thus, despite an almost identical coordination geometry of the E 1 -CO@Fe6 species and the E 1 -N 2 @Fe6 species, favorable binding of N 2 to FeMoco clearly requires some missing element.

N 2 Binding to FeMoco: The E 2 Redox State.
The E 2 redox state has been proposed to contain a hydride based on the state's ability to evolve H 2 (with slow kinetics as recently revealed 9 ) as well as its photochemical properties. 34 The E 2 state thus would appear to share a hydride in common with the N 2 -binding E 4 state; however, in contrast, the E 2 state has not been found to bind N 2 .
Unlike the E 1 state, where computational studies and most experimental studies are in good agreement on the most likely model, the E 2 state offers more structural possibilities. In a recent QM/MM study, 12 we performed a detailed exploration of the energy landscape of the E 2 redox state using a TPSSh-QM/MM protocol. Two isomers were found to be the most energetically favorable: (i) a structure with a hydride bridging Fe2 and Fe6 with a terminal sulfhydryl group on Fe6, E 2 -hyd-SH − @Fe6 (see Figure 6b) and (ii) a nonhydride, doubly beltsulfur protonated state, E 2 -nonhyd (see Figure 6c). The hemilability of the protonated belt sulfur bridge (S2B) allows it to become a terminal sulfhydryl group that stabilizes a bridging hydride between Fe2 and Fe6. Additionally, we included an alternative E 2 isomer, E 2 -SH − @Fe2 (see Figure 6b). This isomer is predicted to be relatively high in energy (11.7 kcal/ mol using the QM-VI region) compared to E 2 -SH − @Fe6, likely due to steric repulsion involving the terminal SH − group and His195. This makes it an unlikely candidate for the E 2 state (previous results 12,30,35,40 also found a similar trend); however, as the results for E 0 and E 1 have revealed, stronger binding is found to Fe6 than Fe2, making N 2 binding to such a geometry with an open coordination site to Fe6 of interest.
These three types of E 2 isomers were thus tested for their ability to bind dinitrogen. Figure 6 shows N 2 binding energies and relative N−N frequencies for E 2 -hyd and E 2 -nonhyd isomers calculated using the QM-I region. We note that binding energy refers here to the energy of binding N 2 to the lowest-energy isomer of either the E 2 -hyd state or the E 2 -nonhyd state. The results reveal a considerably changed picture with respect to N 2 binding compared to the E 0 /E 1 data. There is now a much stronger tendency for binding N 2 according to the E 2 -hyd isomer data. As shown more clearly in Figure 7, binding N 2 to E 2 -hyd-SH − @Fe6 (the most favorable E 2 isomer) at the Fe2 site is uphill by +5.4 kcal/mol (a considerable shift compared to E 0 /E 1 ). Even more favorable is to bind N 2 directly to the alternative isomer E 2 -hyd-SH − @Fe2, giving a single-step exothermic binding energy of −9.8 kcal/mol, implying a starkly different affinity for binding N 2 . Calculating the binding energy relative to the more favorable E 2 isomer (E 2 -hyd-SH − @Fe6) reduces this binding energy down to +1.9 kcal/ mol, which still implies that E 2 -hyd isomers behave quite differently to E 0 and E 1 with respect to N 2 binding.
In sharp contrast to E 2 -hyd models, the data for the E 2nonhyd models in Figure 6 might appear at first glance to be of little interest. The binding energies are more endothermic, closer to the E 1 data in magnitude: being +8.5 (E 2 -nonhyd-N 2 @Fe6) and +6.2 (E 2 -nonhyd-N 2 @Fe2) kcal/mol using the largest QM-VI region (see Figure S11). However, the endothermic binding energies of E 2 -nonhyd models are contrasted with the relative N−N frequencies that are now larger in magnitude than seen before, especially for the Fe6 site, and curiously do not correlate well with binding energies. These E 2 -nonhyd-N 2 @Fe6 states clearly activate N 2 more than other FeMoco states, yet the binding energies are unfavorable. A possible interpretation of this conundrum is that these states are favorable with respect to the Fe-N 2 interaction but N 2 binding leads to such an unfavorable Inorganic Chemistry pubs.acs.org/IC Article electronic structure for the rest of the cofactor that the binding energy remains unfavorable. We note that in our previous study, the total QM/MM energy difference between E 2 -hyd and E 2 -nonhyd was reported to be 0.1 kcal/mol. 12 The energy difference in our present study is predicted to be much larger: 12.1 kcal/mol (QM-VI region). As discussed in the SI, the difference arises primarily due to different DFT methods used (r 2 SCAN vs TPSSh) and also due to a slightly different QM/MM setup. A recent QM/MM study using different functionals (TPSS, r 2 SCAN, TPSSh, and B3LYP) but small basis set (def2-SV(P)) also reports a similar energy difference. 35 Such a large difference between DFT methods that actually predict very similar FeMoco geometries (and energies of E 4 isomers as shown in the SI) appears at first troubling. It should be pointed out, however, that it is not too surprising that the energy difference between two states with electrons either localized in Fe-hydrides vs in Fe d-orbitals for this complicated multimetal cofactor would be challenging and highly sensitive to the DFT method. At this point, the relevance of E 2 -nonhyd as a model for the E 2 state is thus unclear, being so sensitive to the functional. The possible relevance of a N 2 -bound E 2 -nonhyd state in the overall mechanism will be discussed later.
The more favorable N 2 binding energies (although still endothermic) for the E 2 -hyd FeMoco states would seem to be connected to the presence of a hydride bound to the N 2binding Fe ion (Fe2 or Fe6). The presence of a hydride ligand in the E 2 -hyd-N 2 @Fe2 or E 2 -hyd-N 2 @Fe6 leads to unique distorted trigonal bipyramidal geometries. Such a ligand-field is expected to lead to destabilization of the d z 2 orbital, and this can be seen in the QRO-based ligand-field diagrams for these two states in Figure 7b,c. The raised d z 2 orbital level would then no longer be populated, and this appears to lead to the stabilization of intermediate-spin Fe(II) for E 2 -hyd-N 2 @Fe2 and low-spin Fe(II) for E 2 -hyd-N 2 @Fe6 according to localized orbital analysis. The Hirshfeld spin populations of Fe2 and Fe6 when N 2 is bound are 2.10 and 0.08 in agreement with the localized orbital analysis. These two lower spin states arise due to spin-pairing or double occupation of the d xz and d yz orbitals that, as mentioned, are capable of π* backbonding to the N 2 ligand. The changed ligand-field in these states due to the presence of the hydride likely makes these electron configurations more favorable than the otherwise weak-field tetrahedral ligand-field that favors high-spin states.  38 have revealed that the state contains two Fe-hydrides and two S-based protons with the hyperfine tensors consistent with bridging rather than terminal hydrides. The precise nature of the E 4 geometry with respect to the location of the hydrides and protons is controversial. A model 13 with all hydrides and protons on the same Fe2-Fe3-Fe6-Fe7 face, suggested by Raugei and Hoffman, was put forward as being the most consistent with the hyperfine tensor orientations, while Cao and Ryde reported that such a model but with two protons pointing to opposite directions has slightly lower energy, 15 here labeled E 4 -DP-Fe2/6(5) (DP-Fe2/6(5): hydrides on different pairs of Fe ions and the hydride bridging Fe2 and Fe6 pointing to S5A). Cao and Ryde also found another energetically more favorable model 15 than E 4 -DP-Fe2/6(5), here labeled E 4 -DP-Fe2/6(3) (DP-Fe2/6(3): hydrides on different pairs of Fe ions and the hydride bridging Fe2 and Fe6 pointing to S3A), which appeared also to be consistent with the ENDOR data. In a study 14 from our group, we put forward different models, here labeled E 4 -SP, that consist of a terminal sulfhydryl group on either Fe2 or Fe6 and bridging hydrides between Fe2 and Fe6. These E 4 -SP models were energetically favored over other models tested, apparently due to the energetic stabilization associated with opening the protonated S2B belt sulfide and stabilizing bridging hydrides between Fe2 and Fe6. As discussed in the previous study, there are multiple E 4 isomers of this type, i.e., with an open belt sulfide bridge, but with hydrides either bridging, partially bridging, terminal, or even with an H 2 adduct. The most favorable E 4 model found in this work is E 4 -SP-SH − @Fe6 (M S = 3/2) shown in Figure 9c. While this model is predicted to have the wrong spin state (and features one bridging and one terminal hydride), it is only 1.7 kcal/mol lower in energy than the next-lowest energy model E 4 -SP-SH − @Fe6 (M S = 1/2) model, shown in Figure 9d, which is the same model discussed in our previous work 14 (labeled E 4 -o in that study). The E 4 -SP-SH − @Fe6 (M S = 1/2) and E 4 -SP-SH − @Fe2 (M S = 1/2) isomers are the ones of primary interest, as they feature bridging hydrides (and thus in qualitative agreement with the ENDOR data) but also because they feature a more accessible N 2 binding site on either Fe2 or Fe6, respectively.

N 2 Binding to
In this study, we do not seek to come to a clear conclusion regarding the nature of the E 4 state. In the SI, we present additional energetic data that compares the energies of the E 4 -DP-Fe2/6(3), E 4 -DP-Fe2/6(5), and the E 4 -SP models with four different density functionals (r 2 SCAN, TPSSh, 106,107 B97-D3, 92,108,109 and B3LYP* 110,111 ), as shown in Figure S17, which were the four best performing DFT methods in a recent structural benchmarking study 74 of iron−sulfur dimers and Figure 8. N 2 binding energies and relative N−N frequencies (relative to free N 2 , ν NN,free = 2432 cm −1 ) for three types of E 4 models calculated using the QM-I region. All binding energies are relative to the lowest-energy E 4 isomer E 4 -SP-SH − @Fe6 with M S = 3/2 (Figure 9c).

Inorganic Chemistry pubs.acs.org/IC Article
FeMoco E 0 from our group. The results reveal that the four functionals show an overall impressive consistency and they all predict the strong preference of E 4 -SP models over E 4 -DP-Fe2/6(3) and E 4 -DP-Fe2/6(5) models by 15−20 kcal/mol. However, there is clearly a disconnect between the interpretation of the ENDOR data and the calculated energetics that is not easily resolved and additional experiments are likely required for a firm conclusion on the nature of the E 4 state.
Here, on the other hand, we are interested in the binding of N 2 to models of the E 4 state, and E 4 -SP, E 4 -DP-Fe2/6(3), and E 4 -DP-Fe2/6(5) models are of interest. Figure 8 shows N 2 binding energies to each of the E 4 -DP-Fe2/6(3), E 4 -DP-Fe2/ 6(5), and the E 4 -SP models relative to the most stable E 4 model calculated at the r 2 SCAN-QM/MM level of theory. E 4 -SP models and especially models with N 2 bound to Fe6: E 4 -SP-N 2 @Fe6 show considerably different N 2 binding energies. As shown in Figure 9 (presenting the data with a large QM-VI region), a highly favorable single-step N 2 binding energy of −15.2 kcal/mol is found for the E 4 -SP-SH − @Fe2 → E 4 -SP-N 2 @Fe6 binding. However, calculating the energy with respect to the most stable E 4 isomer found, E 4 -SP-SH − @Fe6 (M S = 3/ 2, BS346) (shown in Figure 9c), the energy becomes −4.1 kcal/mol. A similar binding energy for E 4 -SP-SH − @Fe6 → E 4 -SP-N 2 @Fe2 is −4.0 kcal/mol, which shows that Fe2 and Fe6 binding sites cannot be so easily distinguished. In contrast, N 2 binds to Fe2 and Fe6 in E 4 -DP-Fe2/6(3) and E 4 -DP-Fe2/ 6(5) models with single-step N 2 binding energies close to zero kcal/mol (see Figure S18), which, however, is not an improvement over the E 2 models in Section 3.4, perhaps as the number of hydrides at the Fe2/Fe6 binding site remains the same.
Prior to N 2 binding, these models feature local spin states of the 5-coordinated Fe ion of either intermediate-spin or highspin states according to localized orbitals (see Figures S19 and S20). After N 2 binding, creating effectively a 6-coordinated Fe ion, the ligand-field diagrams reveal an approximate octahedral ligand-field with a sizeable gap between t 2g and e g orbitals. The d xz and d yz orbitals are again doubly occupied in the E 4 -SP-N 2 @Fe6 and E 4 -SP-N 2 @Fe2 states. The Hirshfeld spin populations are in good agreement with the orbital analysis, showing small spin populations on the N 2 -binding Fe ion, suggesting a local low-spin state, clearly due to the change in coordination enabled by the hydrides.

N 2 Binding to E 4 -SP in Two α-His195
Protonation States. The protonation state of α-His195 is plausibly assigned as being protonated on the N ε atom in the E 0 state according to the high-resolution X-ray crystal structure. However, it is unknown whether this protonation state persists for different E n states or whether the N δ position could be protonated instead. This histidine residue has in fact been suggested to have a proton transfer role during the reduction of N 2 according to mutation studies. 112 If a proton transfers from N ε to the cofactor (becoming either a hydride or S-bound proton), a simultaneous Grotthuss-type N δ protonation could in fact occur, possibly via a water molecule hydrogen-bonding to N δ in the X-ray structure. This would give an inverse protonation state of His195 (His195-N δ (H)) that could persist depending on its stability.
Our results show that the N ε protonation state is more energetically favorable than the N δ one by >15 kcal/mol for E 4 -SP-SH − @Fe2, E 4 -SP-SH − @Fe6, E 4 -SP-N 2 @Fe2, and E 4 -SP-N 2 @Fe6 cases, shown in Table S6, which is consistent with previous results 12 for E 0 , E 1, and E 2 states; however, the energy difference of these protonation isomers of His195 would be solely determined by the specific protein and cofactor environment, which is currently biased toward the E 0 state in our model (as our model is based on the resting state X-ray structure). It is thus not unreasonable to imagine that a slightly different environment present in a reduced cofactor state could change this protonation state preference. A future QM/MM dynamics study may be able to shed more light on the plausibility of this idea.
Here, we explore how assuming a His195-N δ (H) protonation state instead would affect the binding of N 2 to the E 4 -SP isomers. A doubly protonated His195 residue was also Inorganic Chemistry pubs.acs.org/IC Article briefly explored (see the SI); however, seeing as spontaneous proton transfer from His195 to the cofactor occurs in that case, such a protonation state may not be realistic. Interestingly, Dance has previously described reduced FeMoco isomers containing a doubly protonated His195; the use of a cluster model (instead of a QM/MM model) with a different functional (BLYP instead of r 2 SCAN) is likely the reason for the different prediction. 113 Figure 10 compares how the relative energies of E 4 -SP-SH − @Fe6 and E 4 -SP-SH − @Fe2 isomers as well as N 2 binding energies are affected by the His195 protonation state. First of all, the energy difference between the E 4 -SP-SH − @Fe6 and E 4 -SP-SH − @Fe2 isomers is reduced from 11.1 to 6.7 kcal/mol in going from His195-N ε (H) to His195-N δ (H). Second, the binding preference of N 2 to either Fe2 or Fe6 sites is quite strongly affected by the His195 protonation state, with N 2 @Fe2 and N 2 @Fe6 states being about equally stable for His195-N ε (H), while for His195-N δ (H), stronger N 2 binding is generally found as well as a clearer preference for binding to Fe6. The reasons for these trends are not entirely clear but appear to arise due to steric and H-bonding effects involving SH − and N 2 with the histidine residue, depending on whether the Fe2 site has a terminal N 2 or SH − ligand bound.
Finally, we note that mutation studies where α-Val70 is mutated into α-Ile70 have shown to prevent N 2 reduction from occurring. 45,46 The crystal structure reveals that the additional methyl group in isoleucine compared to valine is positioned directly above the Fe6 site, introducing steric hindering. While  Inorganic Chemistry pubs.acs.org/IC Article the lack of N 2 reduction in this variant suggests Fe6 as a likely potential binding site for N 2 , it is also possible that the bridging → terminal sulfhydryl group conversion (on either Fe2 or Fe6) seen in the E 4 isomers would be similarly sensitive to such steric hindering with the protonation state of His195 further complicating matters.

DISCUSSION
The results of this study highlight the fact that favorable N 2 binding is actually rather difficult to achieve in calculations of FeMoco. The diamagnetic substitution results (Section 3.1) show that a sufficiently reduced Fe ion in the Fe(I) redox state will favorably bind and activate N 2 but such a localized reduced Fe(I) state may not easily form in FeMoco with the strong tendency of the cluster to delocalize electrons over multiple Fe ions. Figure 11 shows the most favorable calculated N 2 -bound states for each redox state of FeMoco according to the results of the previous sections. FeMoco at the E 0 level will only barely give stable N 2 -bound local minima and only in a seemingly unfavorable distorted trigonal bipyramidal geometry of the N 2 -bound Fe ion (with elongated Fe−C and Fe−S bonds). The addition of a proton to the bridging S2B sulfide allows more flexibility in stabilizing N 2 -bound geometries that together with an additional electron in the E 1 state leads to a terminal N 2 -bound structure in a distorted tetrahedral/seesaw geometry. While the E 1 -N 2 structure shows moderate N 2 activation, the binding energy remains unfavorable. However, adding yet another electron and localizing the two electrons in a bridging hydride between Fe2 and Fe6 in the E 2 state changes the picture, with a considerably more favorable N 2bound state, albeit still slightly endothermic. The addition of a single hydride as a ligand appears to enable the formation of a low-spin electron configuration on the N 2 -binding Fe that is now more favorable than before, even though the E 2 -hyd-N 2 state does not activate N 2 more than the E 1 -N 2 state does. As the ligand-field diagrams in Figure 11 show, there is a common theme of spin-pairing in the d xz and d yz orbitals at the N 2bound Fe ion, even in the E 0 -N 2 state, that should enable backbonding to N 2 . However, such a spin-paired electron configuration (beneficial for N 2 back-bonding) may well come at the cost of destabilizing the spin-coupled electronic structure of FeMoco (an Fe ion in a non-high-spin state will have much less favorable antiferromagnetic coupling interactions). The addition of strong σ-donating hydrides to FeMoco may be what makes such low-spin electron configurations more favorable, despite the cost of losing antiferromagnetic interactions. In fact, adding another hydride to FeMoco, i.e., going to the E 4 state, but importantly adding the hydride to the same pair of Fe ions (rather than different pairs of Fe ions as in previously proposed E 4 models), enables the formation of a distorted octahedral N 2 -bound state in a low-spin electronic configuration (again with d xz and d yz orbitals doubly occupied). Interestingly, this N 2 -bound state features the N 2 -bound Fe ion in a S = 1/2 Fe 3+ oxidation state (i.e., more oxidized than the Fe 2+ ions in the other redox states), but unlike the other redox states, the N 2 binding has now become thermodynamically favorable (based on the electronic energy alone, not free energy). It seems that the addition of two hydrides in the E 4 state have now allowed a favorable low-spin electron configuration to form, with occupation of d xz and d yz orbitals, which is useful for N 2 binding, but this time presumably without the associated energy penalty of forming such a configuration like in the other redox states. In this context, it is noteworthy that Peters and co-workers have synthesized two low-spin Fe 2 (μ-H) 2 complexes (either 2Fe(II) or Fe(II)Fe(I)) capable of binding N 2 with a local octahedral geometry. 114 The 2Fe(II) complex with one bound N 2 gives a similar N 2 frequency shift of Δ = −264 cm −1 (relative to free N 2 ), as found in our E 4 -SP-N 2 @Fe6 model (Δ = −254 cm −1 ). There is nonetheless a difference in the coordination environment: sulfide vs phosphine ligation, carbon vs silicon axial ligand, but also the hydride geometries differ (more asymmetric Fe−H bridges in the cofactor model, more acute <H-Fe-H angles and a non-planar Fe 2 (μ-H) 2 unlike the complexes). It is also worth noting the rich literature describing synthetic Fe complexes with N 2 ligands with increasingly more biomimetic ligand environments. 114 A common theme in the results as presented in Figure 11 is the general preference of N 2 binding to Fe6 rather than Fe2. This effect is seemingly at least partly due to Fe6 being more sterically accessible; however, this is further complicated by the fact that both N 2 -bound states and hydride-containing states stabilize an open S2B belt sulfide bridge with a terminal sulfhydryl group on either Fe6 or Fe2, which causes a similar steric effect. The closeness of the His195 residue is the key here, a residue that has been implicated as both a protontransfer residue and being important to N 2 binding. With the results also rather dependent on the protonation state of His195, as shown in Section 3.6, it thus seems difficult to convincingly distinguish between Fe2 and Fe6 at this stage. Fe6 is arguably, however, a more interesting binding site, being close to the Mo ion as well as the homocitrate ligand with a protonated alcohol group aptly positioned for possible proton transfer to the N 2 ligand. Another QM/MM study furthermore has found protonated N 2 -substrates such as diazene to be thermodynamically more stable in the Fe6 position, 54 revising a previous study 124 that favored binding to Fe2. While the use of a QM/MM model in our calculations is important to describe a realistic electrostatic and steric environment around the cofactor, we note that much simpler cofactor-only (same size as QM-I region) calculations with a continuum model (using the same functional and basis set) were also performed for N 2 binding to E 0 , E 1 , E 2 -hyd, and E 4 -SP models in Figure  11. This data, shown in Figure S24 in the SI, reveal the same trend as the QM/MM data, giving further evidence for the importance of hydride ligation for favorable N 2 binding energies. The differences between cluster and QM/MM data Inorganic Chemistry pubs.acs.org/IC Article (primarily affecting Fe2 binding) can be attributed to steric effects. Experimentally, it is now well established that a molecule of H 2 is formed as N 2 binds to the E 4 state and isotope labeling has shown that H 2 is derived from the two hydrides, a kinetic step that only occurs when N 2 is present (without N 2 , E 4 will instead relax by H 2 evolution to the E 2 state via a regular hydride + proton mechanism). 37,38,62 Figure 12 shows the results of N 2 binding to our E 4 -SP models, in both Fe2 and Fe6 positions, followed by the energy of performing the reductive elimination step. Consistent with experiments, the H 2 evolution via reductive elimination is neither strongly endothermic nor exothermic but rather close to being thermoneutral instead (in both Fe2 and Fe6 positions). While the overall reaction energy for E 4 + N 2 → E 4 -N 2 ′ + H 2 is slightly more favorable at the Fe2 position, N 2 is in contrast more activated at the Fe6 position.
The nature of the calculated E 4 -N 2 ′ state after reductive elimination has occurred (whether at Fe2 or Fe6 position) is interesting. As 4 e − and 4 H + are present in FeMoco in the E 4 state of FeMoco (relative to E 0 ) but reductive elimination removes 2 e − and 2 H + via the two hydrides as H 2 , the redox state is now identical to a nonhydridic E 2 state with N 2 bound to Fe2 or Fe6 with an open belt sulfide bridge. Such a state was in fact previously calculated and discussed in Section 3.4 when we considered N 2 binding to E 2 -nonhyd models. The E 2nonhyd state may be related to the experimentally described E 4 (2H)*, a state formed after photoinduced formation of H 2 from two hydrides. 37 We note that E 2 -nonhyd models are in fact not interesting from the point of view of showing any favorable N 2 binding; they are in fact strongly endothermic for N 2 binding, but counterintuitively, the E 2 -nonhyd-N 2 models instead showed the most activated N 2 substrates of all the states calculated in this work. The most activated N 2 substrate is in fact found in the E 2 -nonhyd-N 2 @Fe6 state, which is shown in Figure 12, and data for several electronic states of this model is shown in Figure 6 and discussed in Section 3.4. The relevance of this state, E 2 -nonhyd-N 2 @Fe6, now seems clearer: it is not favorable to bind N 2 to the E 2 redox state directly (even though a nonhydridic E 2 state leads to the strongest N 2 activation); instead, we must first proceed to the E 4 redox state where the thermodynamically favored doublehydride geometries both enable favorable N 2 binding via stabilizing low-spin electronic configurations and also allow access to the N 2 -activating E 2 -nonhyd state via the reductive elimination step involving the hydrides.
It is unclear what happens next. The E 2 -nonhyd-N 2 state having such low affinity for N 2 according to our calculations should dispel the substrate (unless trapped behind a barrier). Perhaps, more likely the N 2 ligand is promptly converted to a more favorable protonated form (N 2 H or N 2 H 2 ), though this would seemingly conflict with the known reversible N 2 /H 2 exchange step in the E 4 state. 44 A recent kinetic study 9 suggests that a slow kinetic step after N 2 binding/H 2 evolution occurs, which might be assigned to a slow N 2 activation/protonation process. Previous calculations from our group suggested that diazene formation (via SH and OH) might already form, although it was calculated to be slightly uphill (+5.6 kcal/ mol). 14 Clearly, mysteries remain about the nature of this key mechanistic step.

CONCLUSIONS
We have systematically investigated N 2 binding to the FeMo cofactor of nitrogenase using a QM/MM model to unravel the nature of how and why N 2 binds and gets activated for protonation. Critical to this work is the use of detailed electronic structure analysis of each broken-symmetry DFT state calculated and the correlation with calculated N 2 binding energies to each redox state of the cofactor. Using a combination of localized orbital analysis, Hirshfeld spin populations, and ligand-field diagrams derived from quasirestricted orbitals (made possible via diamagnetic substitution), we could derive cofactor electron configurations but also obtain a ligand-field level insight into the local Fe spin state changes occurring as a function of hydride and N 2 coordination in the models calculated.
This study builds on models for the E 2 and E 4 states previously introduced by our research group 12,14 but are here analyzed in more detail and compared with the E 0 and E 1 states. The results strongly suggest that it is double-hydride coordination to Fe2 and Fe6 that might be responsible for an increased binding affinity of N 2 to the cofactor by stabilizing a local low spin state in the N 2 -bound state. Overall, this chemically plausible model offers an explanation for why hydride coordination in E 4 allows N 2 binding to occur and furthermore how reductive elimination could be triggered by N 2 binding. Importantly, we also show that the trend of more favorable N 2 binding energy with the E n redox state is reproduced not only with our QM/MM model but also with a QM-I cluster model of FeMoco (lacking the protein environment). This suggests that it is primarily the electronic structure of the cofactor that explains the most basic mechanism of how N 2 might bind to cofactor while the protein environment is likely to influence precisely to which Fe ion N 2 will bind.
There are important caveats to our calculations, however. The results reveal favorable N 2 binding based on electronic energies (approximate enthalpies). It is unclear whether the calculated binding is able to overcome the loss of translational entropy (may also be offset by H 2 evolution), which, based on gas-phase estimates, would be on the order of 10−15 kcal/mol. Recent work by Dance, however, has suggested a lower ∼4 kcal/mol entropic penalty based on N 2 binding from a Inorganic Chemistry pubs.acs.org/IC Article diffusible position within the protein. 59,125 Future free energy simulations will hopefully be able to clarify this. The problem of DFT-method dependency on the results of FeMoco calculations is still far from being fully understood. In this work, we used the r 2 SCAN functional that our previous benchmarking study 74 identified as the best performing density functional for the geometric features of iron−sulfur and iron− molybdenum−sulfur spin-coupled dimers as well as FeMoco itself. Future benchmarking of relevant reaction energies would be desired to estimate better the uncertainty associated with the energies. We also note that the results of this work have revealed local spin-state changes to be especially critical to the stability of N 2 -bound species. Encouragingly, r 2 SCAN is known to be an accurate DFT method according to benchmarking studies of Fe spin crossover complexes. 83 The calculated N 2 -bound species in this work show only weak-to-moderate N 2 activation rather than strong N 2 activation, as defined by Studt and Tuczek. 103 While strong N 2 activation (Δν NN = ∼1000 cm −1 ), typically requiring bridging M-N 2 -M interactions, might be expected to be necessary to activate N 2 for protonation, successful stoichiometric and catalytic N 2 fixation has been achieved for complexes with only moderately activated N 2 . 126 In a recent study 127 by Peters and co-workers on N 2 reduction catalysis with molecular Fe/Mo compounds, the use of a concerted proton-electron transfer mediator, which enables favorable Hatom transfer to the substrate, even avoids strongly reduced metal ion species (that are known to activate N 2 more strongly) and allows H-addition to metal-bound N 2 seemingly without ever activating N 2 strongly. The details of concerted proton-electron transfer to FeMoco are, however, not well understood. The level of N 2 activation in the most favorable N 2 -bound FeMoco species discussed in this work is (Δν NN = 302 cm −1 ) on par with the level of N 2 activation found in, e.g., the Mo-bound catalytic Schrock complex (Δν NN = 340 cm −1 ), 103,128 previously classified as an example of moderate N 2 activation. Much more work is required before the mechanism of nitrogenase can be claimed to be understood, but it is our hope that this study represents a step toward a more detailed understanding of the electronic structure basis of biological N 2 fixation. ■ ASSOCIATED CONTENT
Additional information on QM-regions, electronic structure analysis, orbitals, geometries, broken-symmetry states, and relative energies (PDF)