Introduction

The atmosphere of the earth contains 78% N2, but nitrogen is still a limiting element for plant growth. The reason for this is the strong triple bond in N2, which makes it highly inert. There is only one enzyme that can cleave this bond, nitrogenase (EC 1.18/19.6.1), which is present in a few bacteria and archaea [1,2,3]. The nitrogenases catalyse the chemical reaction.

$${\text{N}}_{2} + 8e^{-} + 8{\text{H}}^{ + } + 16{\text{ ATP}} \to 2{\text{NH}}_{3} + {\text{H}}_{2} + 16{\text{ADP}} + 16{\text{P}}_{i}$$
(1)

The electrons are provided by the Fe protein, which contains a [4Fe4S] cluster and also binds and hydrolyses ATP [1,2,3]. The Fe protein binds to the MoFe protein, which contains the P-cluster, used for electron transfer, and the catalytic FeMo cluster [4,5,6,7,8]. The latter is a complicated MoFe7S9C(homocitrate) cluster, although in some proteins, the Mo ion is replaced by V or Fe [9].

The nitrogenase reaction has been extensively studied with various spectroscopic and kinetic methods [1,2,3]. The reaction is normally described by nine states, E0–E8, differing in the number of delivered electrons (and probably also protons) [10]. It is currently believed that it is the E4 state that binds N2 and that one molecule of H2 dissociates in this binding process [1,2,3]. ENDOR experiments indicate that E4 contains two hydride ions that bridge between two pairs of Fe ions [3, 11, 12].

Many computational studies have also been published with the hope of giving a detailed atomistic and energetic picture of the nitrogenase reaction [3, 13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33]. Unfortunately, they have not reached any consensus so far. On the contrary, different studies have suggested either sequential (distal) or alternative protonation of N2 (i.e. that the first three protons are first added to one of the nitrogen atoms, which dissociates as NH3 before the second nitrogen atom is protonated, or the protons are added alternatively to the two nitrogen atoms so that HNNH and H2NNH2 become intermediates) [14, 15, 17, 18]. Likewise, it has been suggested that N2 binds either end-on or side-on and the binding can be to either one or several metals [14, 15, 20, 33]. Some groups suggest that N2 binds instead to the central carbide ion [17, 18]. There is not even any agreement regarding the structure of E4 or where protons are added beyond the E1 state [16, 20, 23, 25, 34].

We have recently shown that a prime problem of the computational studies is that different density functional theory (DFT) methods give widely different results [23, 25]. In particular, hybrid functionals prefer protonation of the central carbide ion, whereas pure functionals prefer the formation of metal-bound hydride ions. Therefore, the former methods suggest that E4 contains a triply protonated carbide ion [16, 35, 36], whereas the latter functionals suggest that it instead involves two or three hydride ions [20, 23, 25, 26]. The latter methods give results that agree better with the experimental data [3, 11, 12] and they also give geometries of the resting state closer to the crystal structure and a correct spin state on the Mo ion [25]. On the other hand, hybrid functionals seem to give better H2 binding energies [25].

In a series of publications, we have performed a systematic study of the reaction mechanism of nitrogenase [21,22,23, 25, 26]. We have used the combined quantum mechanics and molecular mechanics (QM/MM) approach [37, 38], in which the whole MoFe protein is included in the calculations. Moreover, we have tried to systematically address the problems involved in the computational study of nitrogenase, including the effect of the DFT functional, the basis set, the surroundings and the broken-symmetry state. Our aim is to systematically go through all possibilities or at least to employ some well-defined heuristic approach when the number of alternatives becomes too large. The working hypothesis is that the structures with the lowest QM/MM energy are those involved in the reaction mechanism, as it has also strongly been argued by Siegbahn [35, 39]. We have suggested structures of the E0–E4 states and shown that the previously suggested protonation states, are often quite high in energy [23, 26]. Here, we continue this investigation by studying the binding of N2 to the cluster.

As discussed before [23, 34], this is a formidable task: There are more than 50 possible positions where protons can bind to the FeMo cluster. This gives more than 6 million structures of the E4 state (504) and for each there are at least 35 possible broken-symmetry states. N2 can also bind to a large number of sites (we tested ~ 60 structures in this study). This would give a very large number of possible structures to test (on the order of 1010), which is out of the reach of today’s computational resources. Fortunately, the problem can be strongly simplified. First, it is believed that H2 dissociates when N2 binds [3]. This removes two protons and two electrons, i.e. bringing the FeMo cluster to the same redox and protonation level as in the E2 state, for which there are “only” around 2500 possible structures. Second, it is normally assumed that the two first protonations of the N2 substrate take place directly after the binding [3, 40]. Again, this consumes two protons and two electrons, bringing the FeMo cluster to the same redox level as in the resting E0 state, with no extra protons. The structure of the resting state is known from several crystal structures [4,5,6,7,8].

Therefore, we study here systematically the binding of either HNNH (diazene) or NNH2 to the FeMo cluster of nitrogenase with QM/MM methods. Once this structure is found, it may be possible to work backwards to find possible structures also of the N2-bound conformation and the E4 state.

Methods

The protein

All calculations were based on the 1.0-Å crystal structure of nitrogenase from Azotobacter vinelandii (PDB code 3U7Q) [6]. The setup of the protein is identical to that of our previous studies of the protein [21,22,23, 25, 26]. The entire heterotetramer was included in the calculations, because the four subunits are entangled without any natural way to separate them. The QM calculations were concentrated on the FeMo clusters in the C subunit because there is a buried imidazole molecule from the solvent rather close to the active site (~ 11 Å) in the A subunit. The P-clusters and the FeMo cluster in subunit A were modelled by MM in the fully reduced and resting states, respectively [21].

The protonation states of all residues were the same as before [21]: 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 was modelled in the singly protonated state with a proton shared between the hydroxyl group (which coordinates to Mo) and the O1 carboxylate atom. This protonation state was found to be the most stable one in an extensive QM/MM, molecular dynamics and quantum-refinement study [21] and this protonation state is also supported by another study [41].

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 (not inside the protein [21]) to neutralise the protein and give an ionic strength of 0.2 M [42]. 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 [21].

All MM calculations were performed with the Amber software [43]. For the protein, we used the Amber ff14SB force field [44] and water molecules were described by the TIP3P model [45]. For the metal sites, the MM parameters were the same as in our previous investigation [26]. The metal sites [21, 26] were treated by a non-bonded model [46] and charges were obtained with the restrained electrostatic potential method, obtained at the TPSS/def2-SV(P) level of theory [47, 48] and sampled with the Merz–Kollman scheme [49].

QM calculations

All QM calculations were performed with the Turbomole software (versions 7.1 and 7.2) [50]. We employed two DFT methods, TPSS [47] and B3LYP [51,52,53], and two different basis sets of increasing size, def2-SV(P) [48] and def2-TZVPD [54]. The calculations were sped up by expanding the Coulomb interactions in an auxiliary basis set, the resolution-of-identity (RI) approximation [55, 56]. Empirical dispersion corrections were included with the DFT-D3 approach [57] and Becke–Johnson damping [58], as implemented in Turbomole.

The FeMo cluster was modelled by MoFe7S9C(homocitrate)(CH3S)(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 in the crystal structure [6] were also included, viz. Arg-96 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 (backbones including the Cα and C and O atoms from Ile-355), as well as two water molecules. Moreover, all models included either HNNH or NNH2 binding to the cluster, giving a total of 151 atoms (shown in Fig. 1a). Following extensive Mössbauer, anomalous dispersion and QM investigations [16, 19, 41, 59], we used the oxidation-state assignment \({\text{Mo}}^{{{\text{III}}}} {\text{Fe}}_{3}^{{{\text{II}}}} {\text{Fe}}_{4}^{{{\text{III}}}}\) of the metal ions, as in the resting state, giving a net charge of – 3 for the QM system.

Fig. 1
figure 1

a The quantum system with second-sphere residues marked and b the FeMo cluster with atom names indicated. The view of the cluster is the same as in all other figures in this article

Full size image

Experiments have shown that the ground spin state of E0 is quartet with a surplus of three α electrons [3, 19]. Consequently, we used this spin state for this work (previous studies have shown that structures and energies obtained with different spin are similar and that DFT calculations are not accurate enough to decide the most stable spin state [23]). The electronic structure of all QM calculations was obtained from the fragment approach by Szilagyi and Winslow to obtain a starting BS state [60]. 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 to reach the desired spin state. Such a state can be selected in 35 different ways (\(\frac{7!}{{3!4!}}\)) [22]. The other BS states were obtained by simply swapping the coordinates of the Fe ions [61].

We have thoroughly studied the 35 BS states for several binding modes and how their energies vary with the QM method, the size of the basis set, the geometry and the influence of the surroundings [22, 23]. The conclusion was that the effects of the basis set and the surroundings were restricted (up to 7–11 kJ/mol), the effect of geometry intermediate (up to 37 kJ/mol, but the correlation, R2, was 0.92–0.98) and that the effect of the DFT functional (TPSS or B3LYP) was large (up to 58 kJ/mol). Therefore, we first studied all systems with the same BS state (although it changed during the geometry optimisation for some states). For the best N2H2 binding modes, we performed a systematic study of all 35 BS states (obtained by simply swapping the Fe ions [22, 61]) with both the TPSS-D3 or B3LYP-D3 methods with the def2-SV(P) basis set, using optimised structures. If the BS calculations did not lead to the expected state, we assumed that this BS state is high in energy and it was not further studied. The various BS states are named by giving their number in the Noodleman nomenclature (BS1–10) [30], followed by the numbers of the three Fe ions with minority spin [26, 62], e.g. BS7-346, indicating that Fe3, Fe4 and Fe6 have β spin (the latter three numbers unambiguously define the state; the numbering of the Fe ions is taken from the 3U7Q crystal structure [6] and is shown in Fig. 1b).

For the free N2H2 ligands, the QM system was immersed into a continuum solvent, employing the conductor-like screening model (COSMO) [63, 64], implemented in Turbomole. The default optimised COSMO radii were employed and a water solvent radius of 1.3 Å [65], whereas a radius of 2.0 Å was used for the metals [66]. The dielectric constant was 80 (water).

QM/MM calculations

The QM/MM calculations were performed with the ComQum software [67, 68]. 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 atom in system 1 and it was optionally relaxed by MM. Thus, it included all atoms in residues 59, 61, 62, 65–74, 92, 95–98, 191–199, 226–231, 234, 235, 253–255, 273–282, 300, 353–355, 358–364, 377–383, 385, 386, 401 422–427, 438, 440–444, 450 and 450 from subunit C and residues 97, 98, 101 and 105 from subunit D, in total 87 residues and 35 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). The total system was spherical and non-periodic with 133 919 atoms.

In the QM calculations, system 1 was represented by a wavefunction, whereas all the other atoms were represented by an array of partial point charges, one for each atom, taken from the MM setup. Thereby, the polarisation of the QM system by the surroundings is included in a self-consistent manner (electrostatic embedding). When there is a bond between systems 1 and 2 (a junction), the hydrogen link-atom approach was employed: the QM system was capped with hydrogen atoms (hydrogen link atoms, HL), the positions of which are linearly related to the corresponding carbon atoms (carbon link atoms, CL) in the full system [67, 69]. All atoms were included in the point-charge model, except the CL atoms [70].

The total QM/MM energy in ComQum was calculated as [67, 68]

$$E_{{\text{QM/MM}}} = E_{{\text{QM1+ptch23}}}^{{{\text{HL}}}} + E_{{{\text{MM123,q}}_{{1}} { = 0}}}^{{{\text{CL}}}} - E_{{{\text{MM1,q}}_{{1}} { = 0}}}^{{{\text{HL}}}}$$
(2)

where \(E_{{\text{QM1+ptch23}}}^{{{\text{HL}}}}\) is the QM energy of the QM system truncated by HL atoms and embedded in the set of point-charge modelling systems 2 and 3 (but excluding the self-energy of the point charges). \(E_{{{\text{MM1,q}}_{{1}} = 0}}^{{{\text{HL}}}}\) is the MM energy of the QM system, still truncated by HL atoms, but without any electrostatic interactions. Finally, \(E_{{{\text{MM123,q}}_{{1}} { = 0}}}^{{{\text{CL}}}}\) is the classical energy of all atoms in the system with CL atoms and with the charges of the QM region set to zero (to avoid double-counting of the electrostatic interactions). Thus, ComQum employs a subtractive scheme with electrostatic embedding and van der Waals link-atom corrections [71]. No cutoff is used for any of the interactions in the three energy terms in Eq. 2.

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. For all structures, the QM/MM geometry optimisations were performed using both TPSS-D3 and B3LYP-D3 methods [47, 57] with the def2-SV(P) [48] basis set. Single-point QM/MM energies were calculated also at the TPSS-D3/def2-TZVPD level.

Result and discussion

In this paper, we study the binding of N2H2 (we use N2H2 as a common name of both HNNH and NNH2) to the FeMo cofactor in nitrogenase with QM/MM methods. All structures were optimised with QM/MM, using the TPSS-D3/def2-SV(P) method. All calculations were performed for the quartet spin state, which is the observed spin state for E0 [3]. For the successful optimisations (i.e. leading to the desired structure), we calculated also single-point TPSS-D3 energies with the large def2-TZVPD basis set and reoptimised the structures with the B3LYP-D3/def2-SV(P) approach, because we have found very large differences between the results of pure and hybrid functionals both in this and previous studies [22, 23, 25]. We also repeated the TPSS QM/MM geometry optimisations with the MM system free to relax.

In a recent article, we examined the problem of the BS states and designed a procedure to deal with the BS states [22]. In this study, we do a similar investigation: all calculations were first studied in the BS7-235 state, which is the most stable state for the resting E0 state [22] and also for some of the protonated and reduced states [23], but in some cases, the BS state changed during the geometry optimisation. For all N2H2 binding modes within 20 kJ/mol of the best one (and also several other interesting states), we studied all possible BS states. For simplicity, we discuss only the results obtained with the best BS state, whereas in the tables, results are given also for the BS7-235 state, if available.

Below, we describe the obtained structures for end-on NNH2, side-on HNNH, as well as end-on HNNH in separate sections. We discuss first the results obtained with the TPSS functional and a fixed protein surrounding outside the QM system. Finally, we discuss how the results change when the surrounding protein is allowed to relax or when the functional is changed to B3LYP.

End-on binding of NNH2

To begin with, we studied end-on binding modes of NNH2. We tested both terminal binding to a single Fe ion, bridging between pairs of Fe ions, or binding to four Fe ions. However, we ignored binding to the Fe1 and Mo atoms, because such structures were high in energy in our previous work [23] (we checked some cases in this project to confirm that it is true also for N2H2) and Hoffman and coworkers have argued against N2 binding to Mo [3, 72]. For the terminal binding, we tested only one conformation, viz. trans to the central carbide ion, because such structures were most stable for hydride binding [23] and because the NNH2 molecule is quite large, so there is not much room for alternative binding modes in the protein structure. When NNH2 binds to two Fe ions belonging to the same subcluster (e.g. Fe2–Fe4 or Fe5–Fe7), only one conformation is possible, giving rise to 3 + 3 = 6 possible structures, as shown in Table 1. Each structure is named by giving the number of the metal ions to which the N atoms bind, e.g. Fe2/4, indicating that the N atom binds to both Fe2 and Fe4 (the numbering of the Fe ions is shown in Fig. 1b). However, when NNH2 bridges Fe ions on different sides of the cluster, Fe2/6, Fe3/7 and Fe4/5, two conformations are possible, depending on which side of the belt (µ2) sulfide ions they are located. These two conformations are named after which belt sulfide atom they are directed towards, e.g. Fe2/6(3) or Fe2/6(5), indicating that NNH2 is on the side of the S2B ion that is directed towards S3A or S5A, respectively (again the names of the sulfide ions are taken from the crystal structure [6] and are shown in Fig. 1b). Finally, for the binding to four Fe ions, there is only a single conformation for each face of the cluster, which we call Fe2/3/6/7, Fe3/4/5/7 and Fe2/4/5/6, indicating the four Fe ions involved. This gives 21 possible structures, shown in Table 1.

Table 1 Binding modes of N2H2 tested
Full size table

In many cases, the optimisation failed to give the desired structures, e.g. because the structure reorganised to another structure or the protein structure is too crowded. If so, we tried to obtain it using one or two N–Fe distance restraints and if the geometry optimisation was successful, the restraints were removed and the structures were reoptimised (thus, all presented structures were obtained without any restraints). If this did not give the desired structure, we made no further attempts to get it.

The results of the optimisations are shown in Tables 2 (energies) and 3 (structures). It can be seen that we managed to obtain 19 of the structures. The most favourable binding mode was to Fe6. We found four different structures with NNH2 binding end-on to Fe6. One structure was stabilised by a hydrogen bond between the N atom that binds to Fe6 and the alcoholic proton of homocitrate (HCA; 1.56 Å; Fig. 2a). Consequently, we call this structure Fe6(HCA). Moreover, one of the two protons on NNH2 forms another hydrogen bond to the acetate group of homocitrate (which normally forms an internal hydrogen bond to the alcohol group; 1.60 Å). It had Fe6–N and N–N bond lengths of 1.85 and 1.26 Å. The latter is somewhat longer than in free NNH2, 1.21 Å, calculated at the same level of theory. The Mulliken spin population on Fe6, 1.6, is appreciably smaller than those on other iron ions, 2.5–3.2 (in absolute terms; shown in Table S1). For this structure, we studied the relative stabilities of all possible BS states and found that BS10-135 was most stable (Table S2), 8–23 kJ/mol more stable than BS7-235 with both functionals and basis sets. The Fe6–N bond length increased by 0.03 Å when going from BS7-235 to BS10-135 and the spin population on Fe6 increased by 0.5.

Table 2 Relative energies of the various structures obtained
Full size table
Table 3 Geometries of the various structures obtained
Full size table
Fig. 2
figure 2

The most stable structures for end-on binding of NNH2, obtained at the TPSS/def2-SV(P) level of theory: a Fe6(HCA), b Fe6 and c Fe6(HNNH2). The Fe, Mo, S, C, N, O and H atoms are orange, cyan, yellow, grey, blue, red and white, respectively. Hydrogen bonds to homocitrate and His-195 are indicated with broken lines

Full size image

If the hydrogen bonds between NNH2 and the homocitrate ligand are replaced with the normal internal hydrogen bond within homocitrate (we call that structure simply Fe6, shown in Fig. 2b), the structure was 18 kJ/mol less stable (16 kJ/mol with the large def2-TZVPD basis set). It turned out to be more stable in the BS10-147 state, which was 51 kJ/mol more stable than the BS7-235 state.

However, if the proton is transferred from the homocitrate alcohol group to the substrate, giving HNNH2, still bound end-on to Fe6 [therefore called Fe6(HNNH2)], the structure (shown in Fig. 2c) is actually stabilised by 25–29 kJ/mol, therefore constituting the most stable structure of this type. It had Fe6–N and N–N bond lengths of 1.84 and 1.28 Å (the latter is 0.05 Å longer than in free HNNH2). It forms the same two hydrogen bonds with homocitrate as the Fe6(HCA) structure (although with the opposite polarity for one of them), with H–O distances of 1.69 and 1.60 Å. This structure was also found to be more stable in the BS10-147 state. The spin population was still lowest on Fe6, although it was somewhat higher, 1.8, compared to 2.5–3.2 for the other Fe ions.

Finally, we found a fourth structure, in which there is a close interaction between the N atom that binds to Fe6 and S2B (N–S distance of 1.70 Å). This corresponds to N–S bonds, reflecting that the QM calculations allow for chemical reactions. Such bonds are observed in nine of the structures in Table 3 and they are also characterised by a significantly elongated N–N bond (1.38 Å). They are not among the most stable structures [the present structure, called Fe6(S2B), is the best, 57–79 kJ/mol less stable than Fe6(HNNH2)], but their energies are comparable to many of the other structures without N–S bonds.

Two structures with other binding modes were more stable than the latter structure. The Fe4/5(5) binding mode was 69 kJ/mol less stable than Fe6(HNNH2) (81 kJ/mol with the def2-TZVPD basis set). In this structure, the N atom bridges two Fe ions from different subclusters, but it is much closer to Fe4 (1.81 Å) than to Fe5 (2.18 Å). The N–N bond length is the same as in the Fe6(HNNH2) binding mode, 1.28 Å. The Fe3/4 structure has a similar stability, being 69–72 kJ/mol less stable than Fe6(HNNH2) with the two basis sets. In this structure, the nitrogen atom bridges between Fe3 and Fe4 (i.e. Fe ions from the same subcluster) with distances of 1.83 and 1.76 Å. The N–N bond length is 1.31 Å, i.e. 0.03 Å longer than in Fe6(HNNH2). Interestingly, S4A has dissociated from Fe4 in this structure (3.6 Å distance, cf. Figure S1 in the supplementary material). Both structures were most stable in the BS7-235 state.

The other structures of this type are at least 93 kJ/mol less stable than Fe6(HNNH2). It can be seen from Table 2 that there is no clear relation between the type of structures and their energies, indicating that the stability is mainly determined by whether NNH2 can fit into the structure without clashing with the surrounding protein. In general, structures with NNH2 binding to two Fe ions give longer N–N bond lengths (1.28–1.35 Å) than those binding to one Fe ion (1.25–1.28 Å, disregarding those with close N–S interactions), and those with NNH2 binding to four Fe ions give even longer N–N bonds (1.35–1.39 Å).

Side-on binding of HNNH

Next, we performed a similar investigation for complexes with side-on binding. In this case, we assumed that the substrate binds in the form of cis-HNNH, where the hydrogen atoms are on the same side of the N–N bond and therefore do not interfere with the binding. We tested side-on binding to the same Fe ion, to two or four different metal ions. Each structure is named by giving the number of the metal ions to which each of the N atoms bind, e.g. Fe2Fe2 (same metal), Fe2Fe3 (bridging two metal ions) and Fe2/3Fe6/7 (indicating that the first N atom binds to Fe2 and Fe3, whereas the second binds to Fe6 and Fe7). This way, it can be directly seen if we discuss a side-on structure (two “Fe”) or an end-on structure (one “Fe”). It can also be seen to which metal ions the various N atoms bind. All side-on binding modes to two Fe ions had the Fe–Fe and N–N vectors parallel; the FeMo cluster turned out to be too crowded to allow for any transverse side-on binding modes.

For most side-on structures, we studied only a single conformation (with the two H atoms on HNNH pointing away from the cluster. However, when HNNH bridges the closest Fe ions between the two subclusters, (Fe2Fe6, Fe3Fe7 and Fe4Fe5), two conformations are possible, depending on which side of the belt sulfide ions they are located. Again, the name reflects the direction of the group, e.g. Fe2Fe6(3) or Fe2Fe6(5), indicating that HNNH is on the same side as S3A or S5A, respectively. When HNNH bind to four Fe ions on the three faces of the FeMo cluster, two conformations are possible, depending on whether the N–N bond is parallel with or perpendicular to the approximate C3 axis of the cluster (through the Mo, Fe1 and C atoms). However, these already have different names, e.g. Fe2/3Fe6/7 and Fe2/6Fe3/7, because the N atoms bind to different Fe ions. Finally, HNNH can also bridge more distant Fe ions, diagonally over the cluster face, e.g. Fe2Fe5. In total, we tested 30 possibilities, as it is shown in Table 1.

The results of this investigation are also included in Tables 2 and 3. It can be seen that we found only 20 of the tested structures. In particular, all structures involving binding to two Fe ions within the same subcluster reorganised to other structures. Likewise, only two structures with HNNH binding diagonally between the two subclusters were found.

At the TPSS level, the best binding mode is Fe3/7Fe4/5. As can be seen in Fig. 3, HNNH bridges between the Fe3/7 and Fe4/5 pairs, with the N–N bond perpendicular to the Fe1–Mo axis. All Fe–N bonds are of a similar length, 1.94–1.97 Å. The N–N bond length is 1.40 Å, i.e. appreciably longer than in free cis-HNNH (1.24 Å). This structure is 27–31 kJ/mol less stable than the Fe6(HNNH2) end-on binding mode with the two basis sets. An investigation of all 35 BS states showed that BS2-234 was lowest in energy, 19 kJ/mol more stable than BS7-235 (33 kJ/mol with the big basis set). It had low spin populations on all four Fe ions binding HNNH, 0.8–1.8 (in absolute terms). The spin population on the other three Fe ions was higher, but still rather low, 2.4–3.0.

Fig. 3
figure 3

The most stable structure for side-on binding of HNNH, obtained at the TPSS/def2-SV(P) level of theory, Fe3/7Fe4/5. The water molecule is HOH-525C in the crystal structure [6]

Full size image

The second-best side-on binding mode was similar, viz. the Fe2/6Fe3/7tr structure, where HNNH also binds to four Fe ions, although on another face of the cluster. However, this structure is peculiar in that one of the protons of HNNH has moved to S5A, whereas the ligand has taken a proton from Arg-96, forming a short hydrogen bond to this group (1.57 Å; cf. Figure S2). This structure is 30 kJ/mol less stable than the Fe3/7Fe4/5 structure (the corresponding structure without this transfer, Fe2/6Fe3/7, is 115 kJ/mol less stable). Five structures [Fe2Fe2, Fe3Fe7(3), Fe2/Fe6(3), Fe3/7Fe3(2) and Fe4Fe5(5)] were 67–80 kJ/mol less stable than Fe3/7Fe4/5, whereas the other structures were 110–236 kJ/mol less stable than Fe3/7Fe4/5. Again, structures with HNNH binding to four Fe ions gave longer N–N bond lengths (1.35–1.40 Å) than those binding to two Fe ions (1.27–1.33 Å), with the exception of the Fe4Fe4 complex (1.37 Å).

End-on binding of HNNH

Several of the side-on structures ended up in structures in which HNNH instead bound terminally to a Fe ion. In principle, this is not unexpected, because for free N2H2, cis-HNNH is 60–73 kJ/mol more stable than NNH2, depending on the functional, basis set and whether the calculations are performed in vacuum or in a COSMO continuum solvent with a dielectric constant of 80 (water). Moreover, cis-HNNH is 21–27 kJ/mol less stable than trans-HNNH in vacuum, but only 11 kJ/mol less stable in the water-like continuum solvent. Therefore, we decided to study also such complexes systematically, looking for complexes with cis- or trans-HNNH terminally bound to any of the seven Fe ions and also Mo.

For trans-HNNH, we found six complexes. Interestingly, the complex with trans-HNNH terminally bound to Fe2 [we will call it Fe2(trans) in the following to discern it to the corresponding Fe2(NNH2) and Fe2(cis-HNNH) complexes] turned out to be the most stable complex in this study. It has a N–N bond length of 1.26 Å, i.e. only slightly longer than free trans-HNNH (1.25 Å). The Fe2–N bond is 1.94 Å, which is longer than for all terminal NNH2 complexes (1.79–1.92 Å; the shortest bond if more than one). As can be seen in Fig. 4a, the HNNH group is stacked between His-195 and Ser-278, but it is not stabilised by any hydrogen bonds. It was found to be 9 kJ/mol more stable than the Fe6(HNNH2) complex (3 kJ/mol with the large basis set) and 40 kJ/mol more stable than the Fe3/7Fe4/5 complex (31 kJ/mol with the large basis set). We performed a full investigation of the BS states and the BS10-147 state turned out to be the lowest in energy, 32 kJ/mol more stable than the BS7-235 state. All states are shown in Table S2 and it can be seen that three additional states are within 3–8 kJ/mol, BS7-346, BS10-135 and BS6-157. The BS10-147 state had the lowest spin population on Fe2 (2.1), slightly less than for the other Fe ions, 2.4–3.2 (in absolute terms).

Fig. 4
figure 4

The most stable structures for end-on binding of trans-HNNH, obtained at the TPSS/def2-SV(P) level of theory, a Fe2(trans), b Fe6(trans) and c Fe4(trans)

Full size image

We found five additional complexes with trans-HNNH. The one binding to Fe6 was also quite low in energy, 10 kJ/mol less stable than the Fe2(trans) complex (19 kJ/mol with the large basis set), which is 1–16 kJ/mol higher than the Fe6(HNNH2) complex. As can be seen in Fig. 4b, HNNH forms a hydrogen bond to the alcohol O atom of homocitrate (1.86 Å), but this atom is also involved in an internal hydrogen bond with its own acetate group (1.50 Å) and the other proton of the substrate does not form any hydrogen bond. It turned out to be most stable in the BS10-135 state, but the BS10-147 was only 6 kJ/mol less stable. The spin population on Fe6, – 1.3, was much lower (in absolute terms) than for the other six Fe ions, 2.1–3.2.

The Fe4(trans) complex was also quite stable, 29 kJ/mol less stable than the Fe2(trans) structure (28 kJ/mol with the large basis set). It had a shorter Fe–N bond length of 1.88 Å, whereas the N–N bond is 1.26 Å. As can be seen from Fig. 4c, the HNNH group forms a weak hydrogen bond to S4A (2.31 Å). As expected, it has a low spin population on Fe4 (1.2, compared to 2.4–3.1 for the other Fe ions). The other trans-HNNH complexes were appreciably higher in energy, 69–158 kJ/mol above Fe2(trans).

We studied also the corresponding cis-HNNH complexes. We found seven such complexes, as can be seen in Tables 2 and 3. Again, the complex involving binding to Fe2 [Fe2(cis)] turned out to be most stable, only 17 kJ/mol less stable than the Fe2(trans) complex (12 kJ/mol with the large basis set). The structure (Fig. 5a) and the spin population were very similar to the Fe2(trans) complex. HNNH is still stacked between His-195 and Ser-278, without forming any hydrogen bonds. We performed a full investigation of the BS states, and the best turned out to be BS10-147, 34 kJ/mol more stable than the BS7-235 state. Again, BS7-346, BS10-135 and BS6-157 were close in energy (1–6 kJ/mol).

Fig. 5
figure 5

The most stable structures for end-on binding of cis-HNNH, obtained at the TPSS/def2-SV(P) level of theory, a Fe2(cis) and b Fe6(cis)

Full size image

The second-best structure was Fe6(cis), 48 kJ/mol less stable than Fe2(trans) (32 kJ/mol with the large basis set). It had Fe6–N and N–N bond lengths of 1.90 Å and 1.25 Å and the distal NH group forms a hydrogen bond to the acetate group of homocitrate (1.96 Å). The other structures turned out to be appreciably less stable, 67–187 kJ/mol less stable than Fe2(cis).

Considering the Fe Mulliken spin populations for all types of complexes (shown in Table S1), some general trends can be observed. First, the highest spin (throughout this paragraph, we discuss only the absolute values of the TPSS spin populations) is nearly always found on Fe1, 3.0–3.5. The only exceptions are when N2H2 binds to Fe1 and for the Fe2/4 and Fe3/4 complexes. The latter two complexes also differ that the largest spin population (on Fe3 or Fe2) is lower, 2.8–2.9, i.e. more similar to those of the second and third largest population for the other complexes. As mentioned above, the spin population is typically lowest on the Fe ions binding N2H2, but there are several exceptions, especially when N2H2 binds to several Fe ions. The lowest spin population varies from 0.03 for Fe4 in the Fe4Fe5(5) structure to 2.2 for Fe6 in the end-on Fe4 structure (in which Fe4 has a population of 2.3). The average is 1.2, illustrating that the spin is appreciably higher than for the N2 complexes studied by Dance, who got a spin population below 0.5 in 65% of the studied structures (and below 0.1 in 39% of the structures) [14]. There is no correlation between the lowest Fe spin population and the shortest Fe–N bond (R = –0.03). Fe6 and Fe7 often have relatively low spin populations (averages of 1.8 and 2.0 over all complexes), whereas Fe2–Fe5 have larger populations (averages of 2.5). Mo typically has a low spin population, 0–0.5, with an average of 0.2. However, in the end-on Fe7 complex, it is 0.9 and in the end-on Fe6Mo complex, it is 0.6. In about 75% of the complexes, the spin on Mo is negative.

Calculations with a relaxed protein

In the calculations discussed so far, protein residues outside the QM system in Fig. 1 were kept fixed at their crystal positions. However, we observed a rather large variation in the MM energy of the surrounding protein (i.e. in the \({E}_{\text{MM123,}{\text{q}}_{1}= \text{0} }^{\text{CL}}-{E}_{\text{MM1,}{\text{q}}_{1}= \text{0} }^{\text{HL}}\) term in Eq. 2), up to 77 kJ/mol. This is much larger than in our previous studies (less than 10 kJ/mol) [22, 23], reflecting that the N2H2 ligand is rather large and may clash into the surrounding protein, so that its binding to the FeMo cluster may require significant reorganisation of the surrounding protein residues. Therefore, for all structures in Table 2, we also performed another QM/MM optimisation in which all MM residues with at least one atom within 6 Å of the QM system (always whole residues) were allowed to relax by a MM optimisation in each QM/MM geometry iteration.

The results of these calculations are also included in Tables 2and3. It can be seen that the relative energies change quite extensively (by up to 70 kJ/mol, both increasing or decreasing the relative stability), but the correlation between the two data sets is quite good, R = 0.89. In particular, Fe2(trans) remains the most stable complex and it is stabilised compared to the other low-energy complex, 23, 29, 86, 63, 85, 34 and 91 kJ/mol more stable than Fe6(HNNH2), Fe6(trans), Fe2(cis), Fe4(trans), Fe3/7Fe4/5, Fe6(cis) and Fe6 complexes. Thus, the Fe6(cis) complex becomes the fourth best structure. The Fe6(HCA) complex spontaneously reorganised to the Fe6(HNNH2) structure. In general, the Fe–N and N–N distances do not change much when the surroundings are allowed to relax as can be seen in Table 3. Fe5(cis) is the only structure that changes qualitatively, in that the HNNH ligand dissociates (it was bound weakly already in the structure with fixed surroundings, Fe5–N = 2.40 Å). The N–N distances change by no more than 0.03 Å. Likewise, the shortest Fe–N distance does not change by more than 0.05 Å.

B3LYP results

All results described so far were obtained with the TPSS functional. However, for all structures, we also performed a B3LYP/def2-SV(P) geometry optimisation. The results of these are also collected in Tables 2 and 3. It can be seen that the two DFT functionals give quite similar structures for all complexes. In general, B3LYP gives a slightly shorter N–N bond, by 0.02 Å on average (0.01 Å for the isolated molecules). For the Fe–N distances, the variation is much larger and more varying (with differences from –0.5 to 0.3 Å), but on average, B3LYP gives slightly longer bonds (0.03 Å, median 0.02 Å). As for the relaxed structures, HNNH dissociated from Fe5(cis) with B3LYP, which is the only qualitative difference. However, in many cases, the relative lengths of Fe–N bonds to the same N atom change significantly.

The relative energies show larger and more systematic differences (the correlation coefficient to the TPSS energies is 0.84). B3LYP strongly disfavours all side-on structures, so that there are no such structures within 164 kJ/mol of the best structure. The end-on NNH2 structures (especially those involve more than one Fe ion) are also mostly disfavoured, but to a smaller extent. The Fe2(trans) complex is the most stable structure also with B3LYP. However, the second-best structure is Fe2(cis), which is only 18 kJ/mol less stable. The Fe6(HNNH2) and Fe6(trans) structures are still low in energy, 29 and 35 kJ/mol less stable than Fe2(trans). Next come the Fe6(cis), Fe6(HCA) and Fe4(trans) complexes, 64–75 kJ/mol less stable than Fe2(trans).

For the best structures, we also run B3LYP optimisations with the surrounding protein free to relax. These results are shown in the last column in Table 2. It can be seen that the Fe2(trans) structure is still best, being 34 kJ/mol more stable than the Fe4(trans) structure.

B3LYP and TPSS sometimes give differences in the most stable BS states. For example, the Fe2(trans) and Fe2(cis) structure are most stable in the BS8-236 state with B3LYP, but in the BS10-147 state with TPSS. Likewise, Fe6(trans) is most stable in the BS2-234 state with B3LYP, but BS10-135 with TPSS. This gives some significant differences in the geometry: for the Fe2(trans) complex, the Fe2–N bond length is 0.19 Å longer and the N–N bond is 0.03 Å shorter with B3LYP than with TPSS in the best BS states.

As usual, B3LYP gives larger and more similar Fe spin populations than TPSS (cf. Table S1). For the best Fe2(trans) structure, the spin on Fe2 with TPSS is 2.1, whereas the other Fe ions have spin populations of 2.4–3.2 (in absolute terms). However, with B3LYP, all Fe ions have 3.5–3.7 (3.6 on Fe2).

Conclusions

We have studied the binding of the substrate to the FeMo cluster in nitrogenase. To reduce the very large number of possible structures, we have studied the binding of N2H2 to the E0 state of the cluster, rather than N2 to the E4 state. We have systematically studied side-on binding of HNNH and end-on binding of both NNH2 and HNNH to one, two or four Fe ions. We concentrated the study to the inner Fe ions (Fe2–Fe7), because both in this and previous studies [23], we have found that binding to the terminal Mo and Fe1 ions is unfavourable and Hoffman and coworkers have argued that N2 binding to Mo is unlikely [3, 72].

Interestingly, our results show that the binding of N2H2 is primarily determined by interactions with the surrounding protein and not by the intrinsic stability of the various binding modes. Thus, among the best structures, we find structures representing all four tested binding modes: side-on binding of HNNH, end-on binding of NNH2 and end-on binding of either cis- or trans-HNNH. This shows that it is absolutely necessary to study N2 binding with QM/MM methods; otherwise the results will be strongly biased by the selection of the QM model, unless all residues in the second coordination sphere of any atom in the FeMo cluster are included in the model. On the other hand, this also makes it harder to automatically set up all possible binding structures: in practice, we had to build up each binding mode by hand, carefully considering all surrounding residues. This should reduce the risk that some binding modes have been disfavoured by the use of a poor starting structure.

We find that the Fe2(trans) binding mode is most favourable, i.e. with trans-HNNH binding terminally to Fe2 (Fig. 4a). This mode is stabilised by HNNH being stacked between His-195 and Ser-278. The binding energy of HNNH to this structure (compared to the resting E0 state and trans-HNNH in a water-like continuum solvent) is favourable by 33 kJ/mol. The second most stable structure with TPSS has HNNH2 bound to Fe6. It is formed by a proton transfer from the homocitrate ligand and it is stabilised by two hydrogen bonds to the homocitrate ligand (Fig. 2c). The proton transfer stabilises the structure by 25–29 kJ/mol and the hydrogen bonds by 16–18 kJ/mol. This structure is 3–9 kJ/mol less stable than the Fe2(trans) structure at the TPSS level, but by 23 kJ/mol if the surroundings are allowed to relax. This structure is attractive, because it may explain why homocitrate is a compulsory ligand of the FeMo cluster and cannot be replaced by any related group, except erythro-1-fluorohomocitrate [73].

With B3LYP, instead the structure with cis-HNNH binding to Fe2 is second best, 18 or 46 kJ/mol less stable than Fe2(trans), depending on whether the surroundings are relaxed or not. This structure is fourth best with TPSS, 17 kJ/mol less stable than Fe2(trans). Two structures with trans-HNNH bound to Fe6 or Fe4 are also rather low in energy.

With TPSS, a structure with cis-HNNH bound side-on to the Fe3–Fe4–Fe5–Fe7 face of the FeMo cluster is also rather low in energy, 31–40 kJ/mol less stable than Fe2(trans). However, this structure is strongly disfavoured by B3LYP [183 kJ/mol less stable than Fe2(trans)] and also destabilised if the surroundings are relaxed [85 kJ/mol less stable than Fe2(trans) with TPSS].

Several other groups have studied the binding of N2 to the FeMo cluster. Dance studied the binding of N2 to the Fe2 or (mainly) Fe6 ions for a minimal cluster model with BLYP and the numerical DNP basis set [14, 74]. He compared 54 structures, differing in the number (0–4) and positions of protons bound to the cluster, the N2 binding mode (end-on or side-on) and direction, as well as the spin state, whereas only one BS state was considered, BS7-247. He found that N2 never bridges two metal ions, that end-on binding is more stable than side-on binding and that it is most favourable if N2 binds trans to the central carbide, which typically leads to cleavage of the Fe–C bond (but the carbide ion always retains at least five Fe–C bonds). The most stable structure had three protons (on S2B, Fe2 and Fe6) and N2 bound to Fe6, trans to C in the triplet state. However, in his suggested reaction mechanism, instead a structure with another proton on S3B and N2 binding side-on to Fe6 was selected, although it was 105 kJ/mol less stable. Moreover, as soon as it starts to be protonated, bridging structures were preferred. Therefore, for the HNNH state, he suggested a structure with HNNH bridging Fe2 and Fe6 asymmetrically (one N bridges Fe2 and Fe6, whereas the other binds only to Fe6. This is clearly not the most stable structure in our study, but it differs from our structures in that it contains three additional protons (still on S2B, Fe2 and Fe6). We do not observe any cleavage of the trans Fe–C bonds in our N2H2-bond structures, only a slight elongation (0.05 Å for the best Fe2(trans) structure, 0.13 Å for the Fe6(HNNH2) and Fe6(HCA) structures, but a shortening by ~ 0.02 Å for the Fe4(cis) and Fe3/7Fe4/5 structures, compared to the crystal structure of the resting state [6]). This is probably an effect of the inclusion of the surrounding protein with its steric and electrostatic stabilisation.

Hallmen and Kästner also studied the binding of N2 to the FeMo cluster with the PBE functional and plane-wave basis set for a minimal cluster model of the E2 state [33]. They considered 12 structures, differing in the binding position of N2 (Fe2, Fe3, Fe6, Fe7 or Mo), whereas the two protons were kept on S2B and S5A. In contrast to Dance, they reported several bridging structures, but the structure with N2 binding end-on to Fe7 (not trans to C) and with a cleaved Fe7–S5A bond was best, 10 kJ/mol more stable than a structure with N2 binding end-on to Mo. In an earlier study of the full reaction mechanism (with a central N3– ion), Kästner and Blöchl suggested that HNNH bridges Fe3 and Fe7 first symmetrically and then asymmetrically [32]. Neither of these structures are supported by our calculations, showing that systematic searches of all binding possibilities are needed, as well as a detailed account of the surrounding protein with its sterical and hydrogen-bonding interactions.

The electronic structure of the FeMo cluster is an important ingredient of DFT studies of nitrogenase. Most groups argue that it does not significantly affect the results and that it is enough to restrict the study to a few BS states [14, 16]. Unfortunately, the results in Table S2 indicate that this an oversimplification. The most stable BS can be either BS10-147, BS10-135, BS7-235, BS2-234 or (only B3LYP) BS8-236. Thus, these belong to four out of Noodleman’s ten types of BS states, without any clear logical connection. Moreover, in total, 18 BS states (i.e. more than half of the 35 possible BS states) are within 10 kJ/mol of the best state for at least one of the 18 investigated structures (including all Noodleman states, except BS1 and BS9) and all 35 BS states are within 50 kJ/mol of the best state for at least one structure. This makes it hard to predict beforehand which state will be most favourable. It also indicates that the automatic procedure employed by Kästner and Blöchl [32], indicating that the best states always belong to BS6 or BS7, does not work properly. Even worse, if we select to always use the BS10-147 state, which is best for six structures and within 10 kJ/mol of the best BS state for four additional structures, it is up to 120 kJ/mol worse than the best BS state for other structures, which could lead to erroneous predictions. BS7-235, which is best for six structures and within 10 kJ/mol for additional two, seems to be a slightly better choice, but it can still give errors of up to 89 kJ/mol (or BS2-234, which may give errors of 82 kJ/mol). Therefore, we recommend a complete BS investigation for the best structures, including full geometry optimisations.

Finally, we note that it is satisfying that our all six best structures involve N2H2 binding to Fe2 and Fe6, because experimental observations have suggested the Fe2–Fe3–Fe6–Fe7 face as the reactive side of the cluster [3, 72, 75]. This shows that our QM/MM approach works properly and is accurate enough to find the most reactive sites without restricting the search by experimental information. This makes our structures the most likely candidates for the N2-bound structure of nitrogenase. In future investigations, we will study what implications this finding has for the interpretation of the E4 state and how the N2-bound states may continue to react and form the NH3 product. We will also investigate the possible dissociation of S2B [7, 76], which bridges the same Fe2 and Fe6 ions, forming a natural binding site between these two ions. In several of our studied structures, this ligand dissociates from one of the two Fe ions when the substrate binds (but never from both).