Quantum Mechanics/Molecular Mechanics Study of Resting-State Vanadium Nitrogenase: Molecular and Electronic Structure of the Iron-Vanadium Cofactor.

The nitrogenase enzymes are responsible for all biological nitrogen reduction. How this is accomplished at the atomic level, however, has still not been established. The molybdenum-dependent nitrogenase has been extensively studied and is the most active catalyst for dinitrogen reduction of the nitrogenase enzymes. The vanadium-dependent form, on the other hand, displays different reactivity, being capable of CO and CO2 reduction to hydrocarbons. Only recently did a crystal structure of the VFe protein of vanadium nitrogenase become available, paving the way for detailed theoretical studies of the iron-vanadium cofactor (FeVco) within the protein matrix. The crystal structure revealed a bridging 4-atom ligand between two Fe atoms, proposed to be either a CO32- or NO3- ligand. Using a quantum mechanics/molecular mechanics model of the VFe protein, starting from the 1.35 Å crystal structure, we have systematically explored multiple computational models for FeVco, considering either a CO32- or NO3- ligand, three different redox states, and multiple broken-symmetry states. We find that only a [VFe7S8C(CO3)]2- model for FeVco reproduces the crystal structure of FeVco well, as seen in a comparison of the Fe-Fe and V-Fe distances in the computed models. Furthermore, a broken-symmetry solution with Fe2, Fe3, and Fe5 spin-down (BS7-235) is energetically preferred. The electronic structure of the [VFe7S8C(CO3)]2- BS7-235 model is compared to our [MoFe7S9C]- BS7-235 model of FeMoco via localized orbital analysis and is discussed in terms of local oxidation states and different degrees of delocalization. As previously found from Fe X-ray absorption spectroscopy studies, the Fe part of FeVco is reduced compared to FeMoco, and the calculations reveal Fe5 as locally ferrous. This suggests resting-state FeVco to be analogous to an unprotonated E1 state of FeMoco. Furthermore, V-Fe interactions in FeVco are not as strong compared to Mo-Fe interactions in FeMoco. These clear differences in the electronic structures of otherwise similar cofactors suggest an explanation for distinct differences in reactivity.

Introduction

Nitrogenases are nature’s solution to the difficult problem of converting atmospheric nitrogen into a bioavailable form.[1,2] These enzymes catalyze the reaction of dinitrogen to two molecules of ammonia in an adenosine triphosphate (ATP)-dependent process according to the stoichiometryas recently established for all three types of nitrogenase.[3] For each molecule of reacted dinitrogen, a molecule of dihydrogen is produced. This obligatory dihydrogen formation was hypothesized early on[4,5] as being related to a catalytically active state that included hydrides. Now it is well established that the reductive elimination of two hydrides to form dihydrogen is a part of the mechanism.[6]

Three different types of nitrogenase enzymes exist: molybdenum nitrogenase (Mo-nitrogenase),[7,8] vanadium nitrogenase (V-nitrogenase)[9−11] and iron-only nitrogenase (Fe-nitrogenase),[12,13] which differ in their catalytic activity. Mo-nitrogenase is the most active at ambient temperature[14] and is preferably expressed by organisms. It is also the nitrogenase that has been most extensively studied. V-nitrogenase, on the other hand, has received less attention because of difficulties in expressing it, as it is only produced in Mo-deficient conditions.[10,15] While V-nitrogenase is less active under ambient conditions than Mo-nitrogenase, at cooler temperatures, V-nitrogenase is more active than Mo-nitrogenase.[15] In recent years, V-nitrogenase has gained more attention because of its capability of binding CO in its resting state,[16] reducing CO to hydrocarbons,[17] and even reducing CO2 to CO.[18] Little is known about the third type of nitrogenase, the iron-only nitrogenase. It is even less active than V-nitrogenase and seems to be only expressed in Mo- and V-deficient conditions.[19,20]

Mo-nitrogenase consists of the molybdenum–iron protein (MoFe protein), an α2β2 heterotetramer, and the iron protein (Fe protein), which is a δ2 homodimer. The MoFe protein contains a catalytically active site, the iron–molybdenum cofactor (FeMoco) a [MoFe7S9C-homocitrate] cluster, and the P-cluster an [Fe8S7] cluster.[21] In the dithionite reduced state (resting state), FeMoco exhibits an S = 3/2 spin state, as revealed via electron paramagnetic resonance (EPR) spectroscopy, whereas the P-cluster is EPR-silent (in an integer spin state).[1,22−24] The Fe protein is a reductase containing an [Fe4S4] cluster, which shuttles electrons to FeMoco through the P-cluster in an event driven by the hydrolysis of ATP.[25]

Like Mo-nitrogenase, V-nitrogenase is a heteromultimeric protein complex comprised of the vanadium–iron protein (VFe protein) and the Fe protein. The Fe protein of V-nitrogenase shows high similarity to the Fe protein of Mo-nitrogenase, as recently revealed by crystallography.[26] However, the VFe protein, as it appears in the 1.35 Å resolution X-ray structure (PDB 5NY6),[27] is an α2β2γ2 heterohexamer, containing two more peptide chains than the MoFe protein. The role of these additional chains is currently not clear. The VFe protein contains an iron–vanadium cofactor (FeVco) and a P-cluster, with the latter being more or less structurally identical with its MoFe counterpart but with different redox behavior.[28] While the P-cluster of the MoFe protein is EPR-silent in a dithionite solution, there have been reports on S = 1/2[29] and 5/2[30] EPR signals in similar VFe dithionite preparations, which have usually been attributed to a singly oxidized VFe P-cluster, suggesting that a VFe P-cluster is more easily oxidized than a MoFe one.[28,30]

FeMoco of the MoFe protein is now well characterized via high-resolution crystallography and spectroscopy. The cofactor resembles two fused iron–sulfur cubanes ([MoFe3S3] and [Fe4S3]) but with seven Fe ions, nine sulfides, a Mo ion, and a central carbide. It contains a Mo-bound homocitrate ligand and is bound to the protein via a cysteine residue (to one of its Fe atoms) and a histidine residue (to the Mo atom). The overall oxidation state of FeMoco, and the local oxidation state of individual metals, has been a subject of debate. In recent years, the literature has converged on a [MoFe7S9C]− resting state via analysis of the 57Fe Mössbauer isomer shifts,[31,32] Spatially resolved anomalous dispersion (SpReAD) refinement of the MoFe protein[33] and a computational study by us[34] [we found that quantum mechanics/molecular mechanics (QM/MM)-optimized structures only favored this charge state when compared to the 1.0 Å crystal structure of the MoFe protein (PDB 3U7Q)].[35] A Mo(III) oxidation state was discovered in FeMoco via Mo X-ray absorption spectroscopy (XAS) and X-ray magnetic circular dichroism experiments and theoretical calculations.[36−38] A 3Fe(II)4Fe(III) oxidation state is suggested by SpReAD[33] and Se XAS[39] experiments, while theoretical calculations suggest more delocalization of electrons.[32,34]

Through multiple spectroscopic studies, FeVco was found to be structurally similar to FeMoco, with a proposed VFe7 metal core[40−43] and also featuring a central carbide, as evidenced by Fe X-ray emission spectroscopy.[44] While the first crystal structure of the VFe protein confirmed the overall similar nature of FeVco compared to FeMoco, the crystal structure surprisingly revealed a 4-atom bridging ligand between two Fe atoms (Fe4 and Fe5) instead of a bridging sulfide. The electron density map suggested that either a nitrate (NO3–) or a carbonate (CO32–) is present (Figure ),[27] and this same ligand also appears in another crystal structure of the VFe protein.[45] While spectroscopic confirmation of the carbonate/nitrate ligand is lacking for the VFe protein in solution, it will be assumed for the purposes of this study that the carbonate/nitrate ligand, as revealed by the X-ray structure, is an integral part of FeVco that is always present.

Figure 1

Left: FeVco within the protein matrix as it appears in the PBD5N6Y X-ray structure.[27] Right: FeVco with atoms labeled. The 4-atom bridging ligand (shown as XO3) that is proposed to be either carbonate or nitrate.[27]

FeVco, like FeMoco, appears to have a noninteger spin of S = 3/2 in a dithionite solution (in its resting state),[29,46] analogous to synthetic [VFe3S4]2+ cubanes synthesized by Holm and co-workers.[47,48] This assignment is complicated, however, by the occurrence of multiple S = 3/2 EPR signals with different rhombicities as well as additional S = 1/2 and 5/2 signals. Although these are plausibly attributed to one-electron-oxidized P-clusters, some debate remains about the origin of these signals.[28]

In the synthetic V-cubanes, the V atom is in a V(III) oxidation state and the Fe atoms in Fe(2.5) and Fe(II) oxidation states according to 57Fe Mössbauer studies.[48] From V XAS, comparing FeVco and a [VFe3S4]2+ cubane, the V ion of FeVco was also found to be in a V(III) oxidation state.[40,41] A V(III) oxidation state results in a d2 configuration, in contrast to the d3 configuration for Mo(III) in FeMoco and [MoFe3S4]3+ cubanes.[36,49] A joint Fe high-energy-resolution fluorescence detection (HERFD)-XAS and density functional theory (DFT) study comparing [MoFe3S4]3+/[VFe3S4]2+ cubanes and MoFe/VFe proteins suggests a more reduced Fe part in FeVco compared to FeMoco.[50]

In our previous theoretical studies[34,51] of FeMoco in the MoFe protein, we demonstrated that the geometric structure of FeMoco is well described by our computational protocol by comparison to the 1.0 Å X-ray structure of the MoFe protein. Our protocol accounts for the protein environment via a systematically improvable QM/MM model and describes the electronic structure via broken-symmetry DFT (BS-DFT) calculations using the TPSSh exchange-correlation functional, which we have found to describe the complex electronic structure of the cofactor better than other functionals.[51] Furthermore, we have shown that the calculated structures are highly sensitive to the redox state of the cofactor and that the charge state of FeMoco could be unambiguously determined by the structural comparison. The analysis furthermore indicated a specific electronic state (BS determinant) to be in better agreement with the experimental structure than the other low-lying states.[34]

In comparison to FeMoco, few computational studies have focused on FeVco,[44,50,52−55] and to the best of our knowledge, only a single study[55] from our own group has presented calculations on FeVco that included the new 4-atom ligand (carbonate or nitrate). Our previous study presented calculations on a nonresting ligand-bound state and focused primarily on the assignment of a light-atom ligand (NH or OH) bridging Fe2 and Fe6 of FeVco in a recent crystal structure.[45] Our QM/MM calculations found a structure with an OH ligand to be more consistent with the crystal structure than an NH ligand. The results were not as clear-cut, however, for the identity of the 4-atom ligand bridging Fe4 and Fe5, with calculated structures slightly favoring carbonate over nitrate.

In this study, we present QM/MM calculations of FeVco in the VFe protein in its resting state. The aim is to characterize the electronic structure of the cofactor and compare it to the resting state FeMoco in the MoFe protein. Thus, we seek to critically examine, via comparison to the recent crystal structure,[27] whether the calculations are more consistent with a carbonate or a nitrate ligand and whether the Fe redox state of the cofactor is more reduced or oxidized than that in FeMoco. We also set out to understand the differences in the electronic structures between both cofactors that are known to give rise to different reactivities.

Computational Details

MM Model Preparation

The VFe protein was first modeled classically using the 1.35 Å X-ray structure of the VFe protein from Azotobacter vinelandii as a starting point.[27] The protocol is similar to that of our previous model for the MoFe protein.[34] All molecules and residues present in the crystal structure were included, and no attempt was made to model missing residues. GROMACS, version 5.1.4,[56−58] was used to prepare the MM model and add missing H atoms. The VFe protein is an α2β2γ2 heterohexamer (it contains two VnfDKG trimers), and the protonation state of titrable residues was determined in a single αβγ trimer by visual inspection of the hydrogen-bonding patterns [see the Supporting Information (SI) for specific assignments]. It is assumed here that the protonation state of the titrable residues is the same in both αβγ subunits. The CHARMM36 protein force field[59] was used in all MM and QM/MM calculations (see the SI for details on force-field parameters for inorganic residues). The protein was solvated in a cubic periodic box and as the overall charge of the protein was 62 -, the charge was neutralized by adding 62 Na ions ions, giving a total system size of 531080 atoms. The MM model was then simulated in the NVT ensemble at 300 K for 5 ns. Further details on the MM setup are given in the SI.

QM/MM Preparation and Calculations

A spherical QM/MM model (Figure ) was generated from a snapshot from the MM molecular dynamics trajectory). All residues from a single αβγ trimer (VnfDKG) and all residues from the second β monomer (VnfK) of the VFe protein are in the QM/MM model (chains A–C and E, as labeled in PDB 5N6Y),[27] i.e., an αβ2γ heterotetramer. Additionally, all water molecules and other crystallographically determined ions within ∼42 Å of the carbide of a single FeVco unit, as well as 35 Na ions, are included to keep the system charge-neutral. See the SI for further information. The QM/MM model size is 32562 atoms. An alternative QM/MM model consisting of the full protein was also prepared, as detailed in the SI and discussed in the Results and Discussion section.

Figure 2

Left: 531080-atom MM model. Right: 32562-atom QM/MM model.

Chemshell, version 3.7,[60,61] was used for all QM/MM calculations using ORCA, version 4.1.0[62,63] (unless otherwise stated), as the QM code within an electrostatic embedding QM/MM coupling scheme. Broken-symmetry solutions of FeVco were found with the spin-flipping procedure, as implemented in ORCA from a high-spin ferromagnetic state (MS = 35/2, 34/2, or 33/2, depending on the redox state).

Link atoms were used to terminate the QM/MM border using a charge shift scheme, as implemented in Chemshell.[60] For the QM calculations, we used the same methodology as that in previous studies,[34,51,55] the TPSSh hybrid density functional[64,65] with a relativistically recontracted def2 Ahlrichs basis set[66,67] with a triple-ζ ZORA-def2-TZVP for V, Fe, S, carbide, and the XO3 ligand, whereas a double-ζ ZORA-def2-SVP was used for other atoms.[67] The RIJCOSX approximation[68,69] was used with the default grid setting to speed up the Coulomb and exchange integrals. The D3BJ dispersion correction[70,71] and ZORA relativistic approximation[72,73] were also used. The MM calculations used the CHARMM36 force field within the DL_POLY program,[74] as implemented in Chemshell, and QM/MM geometry optimizations were performed with DL-FIND.[75]

The QM/MM optimizations used an active region of 1038 atoms (except where otherwise stated), with all other atoms of the model frozen. Three QM regions were used, with the smallest region being 57 atoms, the intermediate region being 83 atoms, and the largest region being 181 atoms (link atoms are not included in the atom count); see Figures S1–S3. Single-point QM/MM calculations on the crystal structure geometry were calculated using ORCA with the MM point charges included. The VMD program[76] was used to render the figures of molecular structures in this Article and to calculate RMSDs of the structures. Polarized QM energies (i.e., QM energies polarized by the MM point charges) are primarily discussed instead of total QM/MM energies.

Results and Discussion

A computational model of the resting state FeVco depends on the identity of the 4-atom ligand, redox state, total spin state (MS value in BS-DFT), and actual spin-coupling configuration (an electronic or a broken-symmetry state in BS-DFT). We are describing the spin coupling of this complicated system via BS-DFT, which unfortunately is not capable of describing pure spin states. Nonetheless, for the related FeMoco system, BS-DFT has been rather successful for describing its various geometric and spectroscopic properties,[32,44,50,77] and the landscape of broken-symmetry solutions has been extensively explored in multiple studies.[53,78,79] This, however, is not the case for FeVco since the discovery of the 4-atom ligand. We will thus start our discussion by exploring the broken-symmetry landscape of FeVco via single-point calculations on the X-ray crystal structure for various models (section 1), before moving on to QM/MM geometry optimizations for the most plausible broken-symmetry states compared to the crystal structure (section 2). We note that the X-ray structure of the VFe protein has a resolution of 1.35 Å[27] compared to the X-ray structure of the MoFe protein with 1.0 Å resolution.[35] As will be shown, the lower resolution is still sufficient for distinguishing between different cofactor models. The electronic structure of the final FeVco model is then analyzed in detail and compared to the more studied FeMoco (section 3). Finally, we discuss the problem of redox stability of FeVco in these calculations (section 4)

Broken-Symmetry Solutions of FeVco (Using the X-ray Structure)

Noodleman and co-workers originally proposed 10 different BS solutions for FeMoco,[79] assuming 3-fold symmetry of the cofactor. For FeMoco and FeVco, this symmetry is broken, both by the ligation on Mo/V and by the protein environment. In recent studies that account for the protein environment,[51,80−82] it has become common to consider all 35 BS ways of flipping the Fe ions in FeMoco, and as shown in our QM/MM study on FeMoco,[34] spin-isomeric solutions (i.e., solutions equivalent under 3-fold symmetry) give important geometric differences. In the case of FeVco, the new 4-atom bridging ligand further breaks the symmetry, and hence it becomes necessary to consider all possible broken-symmetry solutions. We label the broken-symmetry solutions according to the Noodleman classification[79] but also according to which Fe ions are spin-down, e.g., “BS7-235” (BS7 class according to Noodleman), and where Fe2, Fe3, and Fe5 (crystal-structure numbering) are spin-down. Because the V ion is part of the spin coupling in the cluster, it could also be considered part of the spin-flipping problem. In practice, we have found (similar to FeMoco) that the V ion will always flip to its lowest configuration on its own (see the SI for more information).

Because the charge, ligand, and metal oxidation states of FeVco in the X-ray structure are not completely clear, we will consider three different redox states for each carbonate and nitrate ligand model: [VFe7S8C(CO3)]0,1–,2– and [VFe7S8C(NO3)]1–,0,1+. We note that an Fe XAS study showed that FeVco (in the VFe protein) in its S = 3/2 resting state likely contains a more reduced Fe part compared to FeMoco.[50] While the protein crystals in the X-ray crystallography study were not subject to an EPR investigation, the dithionite conditions employed in crystallization of the protein should result in crystallized resting state (i.e., the S = 3/2 state). We considered two redox states compatible with S = 3/2 ([VFe7S8C(CO3)]0,2– and [VFe7S8C(NO3)]1–,1+) and one integer-spin redox state ([VFe7S8C(CO3)]− and [VFe7S8C(NO3)]0) for comparison. All models with different redox states, different ligands, and different total spins (MS in BS-DFT) and their energies are tabulated in Table S1.

The electronic energies of the 35 broken-symmetry solutions (single-point QM/MM energies on the X-ray structure; QM region of 57 atoms) for each model are shown in Figure . For simplicity, only the lowest-energy total spin states (MS), where the lowest-energy BS solution is found, are shown (see Figures S4–S18 and Table S1 for data on other spin states). In the following discussion, we use the abbreviation [V-XO3], where V indicates the VFe7S8C part of the cofactor, X is either C or N in the 4-atom ligand, and n indicates the total charge.

Figure 3

Relative energies (polarized QM energies) of all 35 broken-symmetry solutions of FeVco. All calculations were performed on the crystal structure geometry of the cofactor, with MM point charges included.[27] In all cases, the BS7-235 state was found to be favored, and all energies are therefore relative to the BS7-235 solution for each redox state and ligand model. (a) Relative energies for carbonate models. (b) Relative energies if nitrate is the 4-atom ligand. The QM-region size in the QM/MM model is 57 atoms.

In the case of [V-CO3]0, the MS = 3/2 BS7-235 solution is favored, while MS = 1/2 BS7-235 (+3.27 kcal/mol), BS10-147 (+4.95 kcal/mol), and BS10-146 (+5.06 kcal/mol) are the second, third, and fourth lowest in energy, respectively. The BS solution that is the second lowest for MS = 3/2 is BS6-167 (+6.09 kcal/mol), whereas the BS7-346 (+7.63 kcal/mol) and BS7-247 (+8.30 kcal/mol) with MS = 3/2 are even higher in energy.

For the integer-spin redox state [V-CO3]−, BS7-235 with MS = 1 is favored, followed closely by MS = 2 BS7-235 (+1.51 kcal/mol; Figure S5), MS = 2 BS7-346 (+8.67 kcal/mol; Figure S5), and MS = 2 BS7-247 (+9.37 kcal/mol; Figure S5). For [V-CO3]2–, BS7-235 with MS = 3/2 is favored, followed by BS8-237 (+7.58 kcal/mol), BS7-247 (+7.83 kcal/mol), and BS7-346 (+8.18 kcal/mol). The MS = 1/2 BS8-236 (+11.08 kcal/mol) is the lowest MS = 1/2 solution.

Exchanging the CO3 ligand for a NO3 ligand results in some changes to the energy landscape. The BS7-235 solution with MS = 3/2 is still preferred for [V-NO3]+, followed by MS = 1/2 solutions BS7-235 (+2.51 kcal/mol), BS10-146 (+4.22 kcal/mol), and BS10-147 (+4.70 kcal/mol), whereas the BS6-167 (+6.50 kcal/mol) solution with MS = 3/2 is the fifth lowest in energy.

For [V-NO3]0, the MS = 1 BS7-235 is favored, with MS = 1 BS10-147 (+6.10 kcal/mol) being the second lowest in energy, whereas MS = 1 BS10-146 (+6.38 kcal/mol) and MS = 1 BS6-167 (8.10 kcal/mol) are the third and fourth lowest in energy, respectively.

For [V-NO3]−, MS = 3/2 BS7-235 is the lowest in energy, with BS7-346 (+3.30 kcal/mol) being second lowest, whereas BS7-247 (+4.09 kcal/mol) and BS8-237 (+4.41 kcal/mol) are the third and fourth lowest in energy, respectively.

Overall, the BS7-235 broken-symmetry solution is always favored, regardless of whether a CO3 or a NO3 ligand is considered or what the redox state is. Unlike FeMoco, the other spin-isomeric BS7 solutions (BS7-247 and BS7-346) are usually higher in energy for FeVco than other non-BS7 solutions. For FeMoco ([MoFe7S9C]−; MS = 3/2), the BS7 class of solutions is the lowest in energy[34,79] and energetically indistinguishable but does result in distinct geometric differences. As discussed in our previous study,[34] the FeMoco BS7-235 solution yields a calculated structure in very good agreement with the geometry of the crystal structure (and captures trends in the metal–metal distances not seen for the other BS solutions) but is not quite the lowest in energy (being 0.7 kcal/mol higher in energy than BS7-346).

The functional dependence of these single-point energy calculations was explored by redoing the calculations on [V-CO3]2– with the nonhybrid TPSS functional and the 20% Hartree–Fock (HF) exchange hybrid B3LYP functional (Figure S19); BS7-235 was still always preferred. Calculations were also performed without an explicit protein environment using a continuum solvation model (CPCM with a dielectric constant of 4) instead of QM/MM (Figure S20), but the preference for BS7-235 was retained.

Mulliken spin populations of the different [VFe7S8CXO3] models are tabulated in Tables S4–S17. We note that if the self-consistent field did not converge to the specific BS solution attempted but rather another one, then no attempt was made to force convergence.

QM/MM-Optimized Models

As discussed, the BS7-235 solution is overall favored by a few kilocalories per mole for all redox states and ligands considered. When considering the other BS7 solutions, this preference for BS7-235, compared to BS7-346 and BS7-247, seems to be an effect exerted by both V and XO3, which stabilize this particular BS solution (see “the BS7 solutions” chapter in the SI). Because of this strong preference, we will only discuss the BS7-235 solution from now on. QM/MM geometry optimizations were performed for the BS7-235 solution for each ligand and each redox state considered: [V-CO3]0,1–,2– and [V-NO3]1+,0,1–. For the noninteger spin redox states, we only considered a MS = 3/2 spin state due to the strong preference for this spin state (and due to the experimental spin state being S = 3/2), as previously discussed, while for [V-CO3]−, both MS = 1 and 2 spin states were considered (because of the small energy gap between these states).

The root-mean-square deviations (RMSDs) of the optimized cofactors (the [VFe7S8C(XO3)] part) with respect to both cofactors in the crystal structure are shown in Table for two different QM regions. This RMSD definition is used throughout the Article.

Table 1

RMSDsb (in Å) of Different QM/MM-Optimized Geometries Compared to the Experimentally Determined X-ray Structurea,[27]

structure and spin	57-atom QM region	181-atom QM region
[V-CO3]2–, MS = 3/2	0.082	0.079
[V-CO3]−, MS = 1	0.085	0.079
[V-CO3]−, MS = 2	0.111	0.104
[V-CO3], MS = 3/2	0.104	0.091
[V-NO3]−, MS = 3/2	0.095	0.085
[V-NO3], MS = 1	0.095	0.082
[V-NO3]+, MS = 3/2	0.093	0.088

The BS7-235 solution is used for all calculations.

The RMSD is defined as the deviation of the [VFe7S8C(XO3)] part with respect to both cofactors in the X-ray structure.

For a small 57-atom QM region, the optimized [V-CO3]2– model yields a structure that is in best agreement with the X-ray structure, based on the RMSDs. The MS = 1 [V-CO3]− model, however, gives an RMSD value of only 0.003 Å higher. Other models give larger deviations and appear at first glance less likely. The MS = 2 [V-CO3]− model is interestingly in worst agreement with the crystal structure, despite the previous single-point calculations showing the MS = 1 and 2 states as close in energy (1.51 kcal/mol in favor of MS = 1). Once optimized, the energy difference between the two MS states increases to 8.60 kcal/mol for the 57-atom QM-region QM/MM model and further increases to 13.08 kcal/mol for the 181-atom QM-region QM/MM model.

As the QM-region size is increased from 57 to 181 atoms, the spread in the RMSD values for the models becomes smaller. [V-CO3]2– and [V-CO3]− now yield equally good structures, and [V-NO3]0 is only 0.003 Å worse. The simple RMSD metric thus appears to not be beneficial in distinguishing between different FeVco models. However, if instead we focus on comparing individual metal–metal distances of the optimized structures to the X-ray structure (presented in Figure and Table ), a clear distinction between models becomes apparent.

Figure 4

Metal–metal distances of six QM/MM models (QM region: 181 atoms) calculated using the BS7-235 solution compared to the X-ray structure. The green and blue ellipsoids highlight the Fe2–Fe3 and Fe6–Fe7 distances that are strongly affected by the redox state. The red ellipsoid highlights the change in the Fe4–Fe5 distance with the bridging XO3 ligand. The metal–metal distances of the crystal structure are the average of the two cofactors found in the 5N6Y crystal structure. All distances are in angstroms.

Table 2

Various Bond Lengths (Å), Atom–Atom Distances (Å), and a Dihedral Angle (deg) of the Two Instances of FeVco in the PBD 5N6Y Crystal Structure and Values from Relaxed QM/MM (181 QM-Region Atoms) Structures

model	crystal A	crystal B	[V-CO3]2–	[V-CO3]−	[V-CO3]	[V-NO3]−	[V-NO3]	[V-NO3]+
spin (MS)	N/A	N/A	3/2	1	3/2	3/2	1	3/2
BS	N/A	N/A	BS7-235	BS7-235	BS7-235	BS7-235	BS7-235	BS7-235
Fe1–Fe2	2.71	2.68	2.69	2.71	2.71	2.69	2.71	2.70
Fe1–Fe3	2.66	2.65	2.63	2.65	2.65	2.62	2.65	2.65
Fe1–Fe4	2.58	2.57	2.58	2.57	2.60	2.56	2.55	2.56
V–Fe5	2.70	2.69	2.67	2.67	2.68	2.67	2.68	2.67
V–Fe6	2.77	2.75	2.81	2.78	2.82	2.81	2.79	2.81
V–Fe7	2.73	2.74	2.69	2.68	2.72	2.72	2.72	2.75
Fe2–Fe3	2.63	2.65	2.60	2.81	2.77	2.62	2.83	2.78
Fe2–Fe4	2.62	2.63	2.64	2.65	2.68	2.61	2.63	2.66
Fe3–Fe4	2.63	2.64	2.65	2.67	2.68	2.64	2.64	2.65
Fe5–Fe6	2.60	2.60	2.64	2.62	2.65	2.61	2.60	2.62
Fe5–Fe7	2.63	2.64	2.65	2.65	2.65	2.61	2.60	2.60
Fe6–Fe7	2.57	2.57	2.56	2.51	2.71	2.59	2.54	2.73
V–Fe1	7.11	7.09	7.05	7.00	7.00	7.04	7.01	6.99
C–Fe1	3.48	3.49	3.47	3.39	3.44	3.45	3.38	3.42
C–V	3.63	3.60	3.59	3.61	3.57	3.59	3.63	3.57
Fe2–Fe6	2.61	2.63	2.62	2.60	2.60	2.63	2.60	2.60
Fe3–Fe7	2.60	2.60	2.58	2.58	2.57	2.59	2.58	2.57
Fe4–Fe5	2.76	2.78	2.70	2.67	2.63	2.64	2.63	2.59
Fe4–O1XO3	1.97	1.93	1.96	1.93	1.90	2.05	2.04	2.01
Fe5–O2XO3	1.93	1.94	1.97	1.94	1.91	2.05	2.02	2.00
∠Fe4–O1XO3-O2XO3-Fe5	–11.30	–10.93	–13.85	–14.42	–15.09	–11.15	–10.76	–10.74
V–Oalcohol	2.17	2.17	2.16	2.14	2.12	2.15	2.12	2.07
V–Ocarboxy	2.12	2.10	2.08	2.05	2.03	2.06	2.03	2.02
V–NHis423	2.30	2.31	2.21	2.21	2.20	2.21	2.21	2.21
Fe1–SCys257	2.29	2.32	2.29	2.27	2.25	2.28	2.26	2.25
Fe2–S2B	2.21	2.23	2.21	2.17	2.17	2.19	2.16	2.16
Fe6–S2B	2.17	2.21	2.19	2.17	2.15	2.17	2.17	2.14
Fe3–S5A	2.23	2.23	2.20	2.17	2.16	2.19	2.16	2.15
Fe7–S5A	2.27	2.25	2.21	2.19	2.17	2.20	2.18	2.15

We note in this context that there is little variation in the metal–metal distances of FeMoco in X-ray structures of the MoFe protein of differing resolution, as is shown in Table S26. In the comparison of a 1.0-Å-resolution X-ray structure (PDB 3U7Q) to a 1.43-Å-resolution X-ray structure (PDB 4TKU),[84] there is a maximum of 0.04 Å difference in the metal–metal distance, whereas the mean absolute deviation is only 0.01 Å. We therefore consider the 1.35 Å resolution of the VFe protein sufficient for our comparison of computed models to X-ray geometries.

Considering first the X-ray structure, a comparison of distances involving Fe1 and Fe2, Fe3, and Fe4, respectively (in the Fe4S3C cubane part of FeVco), reveals Fe1–Fe4 as the shortest (2.57 Å), Fe1–Fe3 in the middle (2.66 Å), and Fe1–Fe2 as the longest (2.69 Å). All computed models also show this trend, revealing an asymmetry in the Fe–Fe interactions. Turning our attention to the Fe2–Fe3 distance next, for [V-CO3]− and [V-NO3]0, the Fe2–Fe3 distance is significantly overestimated compared to the experiment by 0.17 and 0.19 Å, respectively. The same applies for the [V-CO3]0 and [V-NO3]+ models, where the same distance is overestimated by 0.13 and 0.14 Å, respectively, whereas the more reduced [V-CO3]2– and [V-NO3]− models underestimate it by only 0.04 and 0.02 Å, respectively. Hence, a clear sensitivity of the cofactor geometry with respect to the redox state is revealed.

In the VFe3S3C cubane, the X-ray structure reveals the trend of Fe5–Fe7 being longer (2.64 Å) than the Fe5–Fe6 (2.60 Å) and Fe6–Fe7 (2.57 Å) distances. This trend is not captured at all for the more oxidized models, [V-CO3]0 and [V-NO3]+, which instead show Fe6–Fe7 as the longest distance. The carbonate models [V-CO3]− and [V-CO3]2– capture the trend reasonably well and the nitrate models less well.

These redox-dependent structural changes for the Fe2–Fe3 and Fe6–Fe7 pairs are particularly significant and can be understood. The [V-CO3]2– and [V-NO3]− models are in good agreement with the X-ray structure for these distances. For Fe2–Fe3, the deviations are only 0.04 and 0.02 Å, respectively, but when an electron is removed, the distance increases by 0.21 Å for both the [V-CO3]− and the [V-NO3]0 models.

Similarly, the Fe6–Fe7 distance is affected as the cluster is further oxidized. Upon going from the [V-CO3]− and [V-NO3]0 models to the more oxidized models [V-CO3]0 and [V-NO3]+, this distance is increased by ∼0.2 Å. Clearly, these redox-dependent structural changes correspond to the removal of electrons at specific parts of FeVco and, as discussed later, correlate well with the presence or absence of the mixed-valence delocalized electrons of FeVco.

As for the V–Fe5, V–Fe6, and V–Fe7 distances in the crystal structure, the V–Fe5 distance is the shortest (2.70 Å), V–Fe7 is in the middle (2.73 Å) and V–Fe6 is the longest (2.76 Å). All models capture this trend, but the models do not capture the absolute values well.

Overall, only the [V-CO3]2– model reproduces both the metal–metal distances well enough as well as captures all of the observed trends within the cofactor. On the basis of the specific redox-dependent structural changes, the more oxidized [V-CO3]1–,0 and [V-NO3]0,1+ models can all be safely ruled out. Additionally, the 0.13 Å shorter Fe4–Fe5 distance in [V-NO3]− makes this model an unlikely candidate. Therefore, we conclude that the model that best reproduces the crystal structure is [V-CO3]2–.

Electronic Structure of the E0 State of FeVco in Comparison to the E0 State of FeMoco

Until now, we have not discussed the electronic structure of FeVco in any detail. The electronic structure of iron–sulfur cofactors is complicated[85−87] with FeMoco and FeVco probably being the most complicated examples. Because of the many unpaired electrons, spin coupling, strong covalency, mixed-valence spin delocalization, unusual ligand environment (interstitial carbide), and metal–metal bonding interactions, the electronic structure of these clusters is truly exotic, presenting a true challenge to current quantum chemistry approaches. While the ability of the single-determinant BS-DFT strategy to describe these complex systems should be continuously scrutinized, the good agreement previously obtained for the structure of FeMoco[34,55] and now for the [V-CO3]2– model of FeVco compared to the high-resolution X-ray structure, especially in reproducing trends in metal–metal distances, suggests that we are describing at least aspects of the complicated electronic structure correctly. A simple analysis of the spin density via atomic spin populations of the cofactors is typically not very informative. However, as our previous work on FeMoco[32,34,36] has revealed, utilizing localized orbital analysis of the broken-symmetry determinant allows one to gain detailed insight into the cofactor, and approximate oxidation state assignments are possible despite extensive spin delocalization. As discussed in our QM/MM study of FeMoco, the approximate oxidation state assignments from localized orbital analysis could be correlated well with both the experimental and calculated metal–metal distances. We now present a localized orbital analysis of FeVco and compare it to a localized orbital analysis of FeMoco utilizing the intrinsic atomic orbital–intrinsic bond orbital (IAOIBO) method.[88] While previous work by us used Pipek–Mezey localization,[89] the use of better-behaved IAO charges should reduce possible artifacts in the localization procedure arising from maximizing the atomic charge (with the basis set dependence of Mulliken charges being the major issue) of each orbital.

Figure shows our interpretation of the BS7-235 electronic state of [V-CO3]2– (181-atom QM/MM model) and the BS7-235 electronic state of FeMoco (247-atom QM/MM model from a previous study by us[34]) based on assignment of the IAOIBO localized orbitals to specific atoms. Isosurfaces of the localized orbitals themselves are shown in Figures S21–S23.

Figure 5

Proposed electronic structure as interpreted from IAOIBO localized orbitals of the BS7-235 broken-symmetry determinant for both FeVco and FeMoco computational models. Left: [V-CO3]2– 181-atom-QM-region-optimized QM/MM model. Right: 254-atom-QM-region-optimized QM/MM model of FeMoco (from a previous study).[34] The values in the boxes give the population of the localized orbitals (see Figures S21–23 for the localized orbitals) that indicate the magnitude of delocalization of the minority spin electrons between atoms as well as delocalization between the heterometal (Mo/V) and Fe ions.

The Fe ions in FeMoco/FeVco are all high-spin Fe(III) or Fe(II). Assigning the five majority-spin electrons (either α or β depending on the spin coupling) on each Fe ion is unambiguous because they are always well-localized, and the only difficulty is in assigning the minority-spin electrons.

The minority-spin electrons can either localize on a specific Fe ion [making it localized Fe(II)] or alternatively delocalize between two Fe ions (if the majority spins of both ions are parallel) to create a mixed-valence Fe(2.5)–Fe(2.5) pair. However, depending on competing interactions in the cofactor or the effect of the protein environment, one also observes partial delocalization, which is, furthermore, sensitive to the theory level.

Starting with Fe1 in FeVco and FeMoco, five α electrons are easily assigned via the localized orbital analysis. Additionally, there is a β electron delocalized between Fe1 and Fe4 for both cofactors. In FeMoco, the electron is only partially delocalized (Fe1–Fe4: 0.60/0.36), whereas in FeVco, the analogous β electron is more delocalized with even slightly more Fe4 character than Fe1 (Fe1–Fe4: 0.47/0.51). For FeMoco, the Fe1–Fe2, Fe1–Fe3, and Fe1–Fe4 distances are more or less equal (∼2.66 Å) in the PBD 3U7Q X-ray structure,[35] whereas for FeVco, we observe distinctly different distances. Fe1–Fe4 (X-ray: 2.57 Å) is the shortest distance, Fe1–Fe3 (X-ray: 2.66 Å) is in the middle, and Fe1–Fe2 (X-ray: 2.69 Å) is the longest. This difference between FeMoco and FeVco can be rationalized as being due to the Fe1–Fe4 β electron being more delocalized in FeVco than in FeMoco, resulting in stronger Fe–Fe bonding and a shorter Fe1–Fe4 distance in FeVco than in FeMoco. We speculate that the carbonate ligand, bound to Fe4, is behind this increased delocalization of the β electron, possibly by increased back-bonding between Fe4 and the carbonate.

The minority-spin β electron between Fe6 and Fe7 is almost equally delocalized in FeMoco (Fe6–Fe7: 0.49/0.41) as in FeVco (Fe6–Fe7: 0.50/0.43). When the experimental Fe6–Fe7 distances in both cofactors are compared, they are very similar, but FeVco has a 0.03 Å shorter Fe6–Fe7 distance (X-ray: 2.57 Å) compared to FeMoco (X-ray: 2.60 Å).

The minority-spin α electron between Fe2 and Fe3 is somewhat more localized on Fe3 in FeMoco (Fe2–Fe3: 0.34/0.57) than in FeVco (Fe2–Fe3: 0.38–0.54), but classifying Fe2–Fe3 as a Fe(2.5)–Fe(2.5) delocalized pair still seems appropriate. However, the distance between Fe2 and Fe3 in both FeVco (X-ray: 2.64 Å) and FeMoco (X-ray: 2.65 Å) is longer or about equal to those in Fe2–Fe4 and Fe3–Fe4 in the Fe2–Fe3–Fe4 triangle, which is seemingly in disagreement with the assignment of Fe2–Fe3 being a mixed-valent pair but Fe2–Fe4 and Fe3–Fe4 not. We note, however, that calculations of both cofactors are still in good agreement with the respective X-ray structures.

The localized orbital analysis allows us to rationalize the redox-dependent structural changes shown in Figure for the [V-CO3]2–, [V-CO3]−, and [V-CO3]0 models. As electrons are removed from [V-CO3]2–, we get first a 0.21 Å elongation of the Fe2–Fe3 distance in [V-CO3]−, and this is followed by a 0.20 Å elongation of Fe6–Fe7 in [V-CO3]0. These changes obviously correspond to removal of the minority-spin delocalized electrons between these two pairs of Fe ions. A localized orbital analysis of the [V-CO3]− and [V-CO3]0 models confirms this interpretation (Figures S22 and S23). Interestingly, the minority-spin on Fe5 remains in these models. This suggests that the delocalized electrons in FeVco occupy higher energy levels than the localized electron on Fe5.

Analyzing the heterometal orbitals, we see first that the oxidation state assignment of V in FeVco appears to be unambiguous based on the presence of two strongly localized β-spin electrons (and no appreciable V–Fe delocalization) to give a V(III) 3d2 assignment. This is in contrast to Mo in FeMoco, where an unusual Mo(III) 4d3 configuration is observed with an appreciable Mo–Fe delocalization. As previously discussed, the Mo(III) 4d3 configuration is a highly unusual one, breaking Hund’s rule, and has been referred to as a “non-Hund” configuration.[34,36,50] Furthermore, the 4d electrons of Mo are not as localized compared to the 3d electrons of V, with slightly more pronounced Fe character of the β orbital; Mo–Fe6, 0.72/0.22; Mo–Fe7, 0.68/0.27. The α electron on Mo in FeMoco is even more delocalized (Mo–Fe5: 0.50/0.46), which is suggestive of Mo–Fe bonding. The electronic structure of Mo in FeMoco is thus more ambiguous from the localized orbital analysis than that of V in FeVco. We note that V K-edge XAS of the VFe protein supports the V(III) assignment in FeVco,[41] and Mo K- and L-edge XAS supports the Mo(III) assignment.[36,37] More recently, experimental support for the unusual non-Hund configuration at Mo comes from Mo L-edge X-ray magnetic circular dichroism spectroscopy of a related [MoFe3S4] model cubane that shares the unusual electronic structure of FeMoco.[38]

In the X-ray structure of FeMoco (PDB 3U7Q), the Mo–Fe5 distance (X-ray: 2.73 Å) is longer than the Mo–Fe6 (X-ray: 2.67 Å) and Mo–Fe7 (X-ray: 2.68 Å) distances. In our 247-atom QM-region QM/MM model of FeMoco, the calculated Mo–Fe5 distance is also longer (2.71 Å) than the Mo–Fe6 and Mo–Fe7 distances (2.65 and 2.62 Å, respectively). Although the α electron of Mo is somewhat more delocalized than the β electrons of Mo, in our computed model, the Mo–Fe5 distance ends up being the longest of the three. This may be related to the higher effective charge at Fe5 [being assigned as Fe(III) rather than Fe6/Fe7 (assigned as Fe(2.5)–Fe(2.5)] or possibly the overall ferromagnetic alignment of Fe5 and Mo.

While V in FeVco lacks a third d electron, an additional localized α electron is instead present on Fe5 (Fe5–V: 0.82/0.16), making it of primarily localized Fe(II) character. In the PBD 5N6Y crystal structure, the V–Fe5 distance (X-ray: 2.70 Å) is shorter than the V–Fe6 (X-ray: 2.76 Å) and V–Fe7 (X-ray: 2.73 Å) distances, in contrast to the analogous Mo–Fe distances in FeMoco. In our 181-atom QM-region QM/MM model of FeVco, the V–Fe5 distance is also shorter (2.67 Å) than the V–Fe6 and V–Fe7 distances (2.81 and 2.69 Å, respectively).

The Mayer bond order (MBO)[90] is a useful metric for the magnitude of bonding between two atoms. The absolute values of MBO between the metals in FeVco and FeMoco are tabulated in Table S20. In the 247-atom QM-region QM/MM model of FeMoco and 181-atom QM-region QM/MM model of FeVco, we observe similar and relatively low MBO values between Fe atoms that do not have a delocalized minority-spin electron (values ranging from 0.21 to 0.31).

Fe atoms that are ferromagnetically aligned and have a delocalized minority-spin electron, on the other hand, show higher MBO values (ranging from 0.41 to 0.47), indicating more bonding character between these Fe atoms compared to Fe interactions with no delocalized electrons.

The MBOs for the heterometal and Fe interaction in FeVco and FeMoco paint a starkly different picture between the two cofactors. The MBOs for V–Fe5 (0.37), V–Fe6 (0.18), and V–Fe7 (0.30) are much lower than the MBOs for Mo–Fe5 (0.47), Mo–Fe6 (0.45), and Mo–Fe7 (0.49), indicating that Mo in FeMoco interacts more strongly with the Fe atoms of FeMoco than V in FeVco.

This difference in interaction between the heterometal and Fe ions in the cofactor is intriguing because the average Mo–Fe distance in FeMoco is shorter than the average V–Fe distance in FeVco. This difference between the average Mo–Fe and V–Fe distances is seemingly in contradiction because one would expect Mo, a larger ion than V, to have longer Mo–Fe distances than the V–Fe distances. In the 1.0-Å-resolution X-ray structure of the MoFe protein (PDB 3U7Q),[35] the average Mo–Fe distance is 2.69 Å, whereas it is 2.70 Å according to extended X-ray absorption fine structure (EXAFS).[91] In comparison, the average V–Fe distance in the 1.35-Å-resolution X-ray structure of the VFe protein (PDB 5N6Y) is 2.73 Å, whereas it is 2.76 Å as interpreted from EXAFS.[40] In our models of FeMoco and FeVco, we underestimate these distances. The average Mo–Fe distance is 2.66 Å (underestimation of 0.03 or 0.04 Å) and the average V–Fe distance is 2.72 Å (underestimation of 0.01 or 0.04 Å). However, the average M–Fe (M = Mo, V) distance agrees with the EXAFS data, where the increase in the average M–Fe (M = Mo, V) distance between Mo–Fe and V–Fe is 0.06 Å, which we capture exactly in our computed models. This difference between EXAFS and the crystal structure for V–Fe may be attributed to the lower-resolution crystal structure (1.35 Å for VFe).

The delocalization or bonding interactions that we observe between the Mo ion and the Fe ions but less so between the V ion and the Fe ions in our computed models are a likely explanation for the difference between the Mo–Fe and V–Fe distances. It seems likely that this effect is at least partially due to more diffuse 4d orbitals of Mo, resulting in larger overlap between Mo and Fe orbitals than V and Fe orbitals.

These differences in the electronic structure between the two cofactors, which are subsequently reflected in the molecular structure, are likely one of the contributing factors to the difference of the catalytic activity between the two cofactors. Aside from differences in the catalytic activity, another important difference between the two cofactors is that FeVco is capable of binding CO in its resting state, whereas FeMoco can only do so under turnover conditions.[92] One hypothesis is that the more ferrous iron environment in FeVco aids CO binding to an Fe ion in the resting state, something not possible in FeMoco. In fact, the resting state FeVco can be interpreted as being analogous to an unprotonated E1 state in FeMoco, and it has been proposed in a recent joint EXAFS–QM/MM study that the MoFe3S3 subcubane of FeMoco receives the added electron when reduced to the E1 state,[93] which is analogous to the reduced VFe3S3 subcubane of our resting state FeVco model.

Stability of the [VFe7S8(CO3)]2– Redox State

The previously discussed results, in our view, clearly establish the resting state structure of FeVco as containing a bound CO3 ligand and a redox state consistent with [VFe7S8(CO3)]2– and a formal oxidation state distribution as V(III)Fe(II)4Fe(III)3. Furthermore, the calculated electronic state according to a BS-DFT description is MS = 3/2 (which is consistent with the experimental S = 3/2 spin state) and having a lowest-energy BS7-235 solution with spin delocalization, as indicated in Figure . Only this specific ligand (CO3), this specific redox state ([V-CO3]2–), and this specific alignment of spins (BS7-235) can account for the specific geometric properties of FeVco seen in the experimental X-ray structure.

Until now, we have, however, not discussed a problem that arises from this more negative charge of the resting state FeVco ([VFe7S8(CO3)]2–) compared to the resting state FeMoco ([MoFe7S9]−), namely, that of unbound electrons. The QM/MM model discussed (calculated using a 181-atom QM region and the TPSSh functional) results in 12 α and 11 β electrons residing in orbitals with positive orbital energies (Figures S25–S27). Positive orbital energies are unphysical in this context because they suggest that spontaneous detachment of the electrons from the cofactor should occur. While unbound electrons are sometimes an indication of an inflexible basis set, because of the large polarized triple-ζ basis set used here, this is unlikely to be the case. We considered three possible reasons for the presence of unbound electrons in the [V-CO3]2– computational model: (i) the cofactor is, unlike FeMoco, protonated in the resting state, (ii) the unbound electrons arise because of the DFT self-interaction error (SIE), and (iii) the cofactor is insufficiently stabilized by the QM/MM model.

Protonated Cofactor

A QM/MM model (58-atom QM region) with an additional proton on the longer carboxylate arm of the homocitrate ligand was found to result in a reduction in the number of unbound electrons to 2 α and 3 β electrons (with a subsequent increase in the RMSD to 0.095 Å compared to 0.080 Å for an unprotonated model). Models with protonated belt sulfides S2B or S5A (Table S25) were found to have no unbound electrons. However, the protonation of a belt sulfide was accompanied by a 0.10–0.13 Å elongation of the Fe–S2B/S5A distances, a difference that would likely be resolved in the crystal structure if it were present. Additionally, the RMSD was considerably worse than the unprotonated model (0.106 and 0.121 Å for S2B and S5A protonation, respectively, compared to 0.080 Å for an unprotonated model). Thus, we consider a protonated resting state FeVco model unlikely.

Self-Interaction Error

The unbound electrons could arise because of the approximate nature of the density functional used. Approximate density functionals do not show the correct long-range asymptotic behavior of the potential and suffer from SIEs.[94−97] While hybrid density functionals reduce the SIE via the inclusion of exact HF exchange and range-separated hybrids can enforce the correct long-range behavior, the low amount employed in most popular hybrid functionals used for transition-metal chemistry (typically 0–25%) is usually not sufficient to resolve the issues. As revealed in Figure S28, while increasing the global HF exchange in the standard hybrid DFT form results in a decrease of unbound electrons of FeVco, this is not a realistic strategy for describing FeVco because the quality of the geometry deteriorates significantly, as revealed by the strong increase in RMSD with respect to the crystal structure. The 50–70% HF exchange required to remove unbound electrons results in an unacceptable RMSD (>0.20 Å), suggesting that the electronic structure is very poorly described with those functionals. This problem of high HF exchange functionals deteriorating the electronic structure of nitrogenase cofactors has been noted previously by us and others for FeMoco.[51,98]

Environmental Effects

Because the spherical QM/MM model does not account for mutual polarization between the QM and MM regions and does not contain the whole VFe protein, we considered the possibility of protein polarization or longer-range bulk electrostatics playing a role in stabilizing the cofactor.

The protein and solvent environment clearly play a role in describing the cofactor, as can be seen by considering a simple 59-atom cluster model of the cofactor in a vacuum; this results in a very unstable cofactor with 195 unbound electrons. Such a vacuum model, furthermore, gives a poor geometric structure (RMSD of 0.166 Å). Using a simple continuum solvation model, COSMO[99] with a dielectric constant of 4, the electronic structure is stabilized considerably, yet 13 unbound electrons remain. A different continuum model, CPCM (using a Gaussian charge scheme and a scaled vdW cavity[100,101]), however, stabilizes the electronic structure, resulting in no unbound electrons. These results thus clearly reveal a sensitivity of the computational model to bulk electrostatics and polarization effects. These continuum model calculations, however, clearly result in cofactor structures inferior to the QM/MM models, as can be judged by the relatively high RMSDs of the cofactor geometries (Table ), being not much better than the vacuum-optimized geometry.

Table 3

Unbound Electrons and RMSD (with Respect to the X-ray Geometry) in Cluster and QM/MM Models of FeVco

	system size (atoms)	QM-region size (atoms)	charge (QM region)	unbound electrons	RMSDc (Å)
vacuumb	59	59	–6	195	0.166
COSMOb	59	59	–6	13	0.140
CPCMa,b	59	59	–6	0	0.141
QM/MM	32562	57	–6	29	0.082
QM/MM	32562	83	–3	32	0.085
QM/MM	32562	181	–3	23	0.079
QM/MM	47516	57	–6	62	0.084
QM/MM	531080	57	–6	0	0.085
QM/MM	47516	181	–3	70	0.080
QM/MM	531080	181	–3	0	0.079

Using a Gaussian charge scheme and a scaled van der Waals surface.[101]

An additional two protons are added to cap residues 423His and 257Cys.

The RMSD is defined as the deviation of the [VFe7S8C(XO3)] part with respect to both cofactors in the crystal structure.

The lack of protein polarization by the MM environment in the QM/MM models could only be partially explored in this study via a systematic increase of the QM-region size in the calculations. Upon going from a 57-atom QM-region model to a 181-atom model, the number of unbound electrons is only reduced from 29 to 23. Longer-range polarization effects affecting the electronic structure stability can presently not be ruled out and could be explored using even larger QM regions or via polarizable QM/MM methods.

The effect of long-range bulk electrostatics was explored in more detail by the creation of an additional QM/MM setup. Instead of a spherical model (Figure , right), centered on the cofactor of half of the hexameric protein, a new QM/MM model consisting of the whole VFe protein was cut from the periodic MM-box model (Figure , left) of the whole solvated VFe protein. A 3 Å solvation shell and all Na+ ions were additionally included in the model, resulting in an overall neutral model of 47516 atoms. This model (Figure S29), consisting of the whole protein, is slightly larger than the spherical model of half of the protein. Interestingly, however, this larger model that avoids cutting protein chains and includes the whole protein and all Na+ ions (as MM point charges) actually results in an increase in unbound electrons to 62 (for a 57-atom QM region). This unfavorable result, nonetheless, suggests a strong sensitivity to the bulk electrostatic environment through the presence and/or distribution of MM point charges. When the size of the MM point-charge environment was increased even further, a model where all 531080 atoms (additional water molecules) of the solvated protein box (Figure , left) were included in the electrostatic embedding was calculated. Remarkably, this results in complete stabilization of all electrons. QM/MM geometry optimizations at this level of theory further revealed almost negligible structural changes compared to the 47516-atom model, demonstrating that these unphysical unbound electrons do not affect the local structure of the cofactor. Calculations with an even larger 181-atom QM region with the extended MM region resulted in a very similar geometry and no unbound electrons.

The problem of unbound electrons in computational models of FeVco thus appears primarily to be a question of insufficient solvent bulk electrostatic effects. Fortunately, insufficient account of this bulk effect has little consequence for the quality of the cofactor geometry (Table ) or the spin-density distribution (see the Mulliken spin populations in Table S26). Clearly, however, a proper account of solvent bulk electrostatics needs to be carefully considered in future calculations of the redox properties of these cofactors.

Conclusion

We have performed, to our knowledge, the first extensive computational study on the resting state FeVco of V-nitrogenase from A. vinelandii using a QM/MM model starting from the 1.35-Å-resolution X-ray structure (PDB 5N6Y). We systematically explored 35 different broken-symmetry solutions for all plausible oxidation states of FeVco, using either nitrate or carbonate as a possible 4-atom bridging ligand, and found that the BS7-235 solution (with Fe atoms labeled 2, 3, and 5 being spin-down) is energetically favored in all cases. Through a careful comparison of the QM/MM-optimized structures to the X-ray structure by analysis of the Fe–Fe and V–Fe distances, we find that only the [VFe7S8C(CO3)]2– model is in good agreement with the X-ray structure. Analysis of the electronic structure of [VFe7S8C(CO3)]2– via the localized orbitals reveals a reduced VFe3 subcubane compared to the MoFe3 subcubane in FeMoco. This is in good agreement with a previous Fe K-edge XAS study of MoFe and VFe proteins[50] and, interestingly, suggests that the Fe electronic structure of the resting state FeVco is analogous to the E1 state of FeMoco, recently characterized by a joint Fe/Mo EXAFS and QM/MM study.[93] The additional electron in the VFe3 subcubane is localized on Fe5, making it locally ferrous, in contrast to the ferric Fe5 in the analogous computational model of FeMoco (as interpreted from the BS7-235 solution). The oxidation state of V is found to be V(III), in agreement with experimental V K-edge XAS studies,[41] and the V(III) ion is found to have less bonding interactions with Fe ions than Mo(III) in FeMoco.

These electronic structure differences between FeVco and FeMoco are proposed as a likely reason for the experimentally observed differences in the reactivity and catalytic activity.

Finally, we discussed the stability of the electronic structure of the FeVco model with respect to unbound electrons. It is found that the more reduced FeVco model (compared to FeMoco) is highly sensitive to solvent bulk electrostatic effects and that a large number of explicit solvent molecules are required to stabilize the highest electron energy levels of the cofactor.