Electronic and Structural Properties of the Double Cubane Iron-Sulfur Cluster

: The double-cubane cluster (DCC) refers to an [Fe 8 S 9 ] iron-sulfur complex that is otherwise only known to exist in nitrogenases. Containing a bridging µ 2 -S ligand, the DCC in the DCC-containing protein (DCCP) is covalently linked to the protein scaffold via six coordinating cysteine residues. In this study, the nature of spin coupling and the effect of spin states on the cluster’s geometry are investigated computationally. Using density functional theory (DFT) and a broken symmetry (BS) approach to study the electronic ground state of the system, we computed the exchange interaction between the spin-coupled spins of the four FeFe dimers contained in the DCC. This treatment yields results that are in excellent agreement with both computed and experimentally determined exchange parameters for analogously coupled di-iron complexes. Hybrid quantum mechanical (QM)/molecular mechanical (MM) geometry optimizations show that cubane cluster A closest to charged amino acid side chains (Arg312, Glu140, Lys146) is less compact than cluster B, indicating that electrons of the same spin in a charged environment seek maximum separation. Overall, this study provides the community with a fundamental reference for subsequent studies of DCCP, as well as for investigations of other [Fe 8 S 9 ]-containing enzymes.


Introduction
Enzymes catalyzing ATP-driven redox reactions usually consist of two components. The first component is an ATPase, which contains an FeS cluster that provides the electron that is to be transferred. The second component is another metalloenzyme, the electron acceptor, which is ready to accept a highly energetic electron since it contains one or more metal centers that cannot easily be reduced at physiological reduction potentials [1]. Recently, the two components of a system of this type were characterized and the electronaccepting component was found to be a novel double-cubane [Fe8S9]-cluster (DCC) [1] ( Figure 1). This cluster catalyzes reductive reactions otherwise associated only with the complex iron-sulfur clusters of nitrogenases [2]. Although double-cubane [Fe8S9]-clusters have been reported only once for a biological system, they have already been synthesized and studied [3][4][5][6].  [1], shows the A-cluster coordinated to sidechains of Cys75, Cys308, and Cys340 and the B-cluster coordinated to sidechains of Cys113, Cys143, and Cys373. Acluster is linked to B-cluster via the μ2-S ligand (S9). Each cubane cluster contains an exchange interaction between the mixed-valence Fe dimers within the cubane (Jcubane) (dimer high spin "up" shown in blue; dimer high spin "down" shown in red). An exchange energy also arises from two mixed-valence Fe dimers interacting via the μ2-S ligand (Jcubane-cubane).
Stack et al. reported the first example of a sulfide-bridged double-cubane cluster in 1989, hypothesizing that such a Fe-S complex may exist in proteins, as it is constructed of discrete cubane units [3,7]. Several variants of the inorganic DCC with different ligands have been studied, showing the lability of the bridging μ2-S [3,4,7,8] and very weak magnetic coupling between the two subclusters. However, the electron transfer is strongly coupled and two distinct redox reactions occur [3]. If such a FeS cluster is embedded in a so-called "double-cubane cluster" protein (DCCP), diverse environments of the two subclusters may lead to differing reduction potentials [1].
Notably, the corresponding cluster in Carboxydothermus hydrogenoformans, DCCPCh, exhibits similarities regarding structure and reactivity with the FeMoCo iron-molybdenum cofactor of nitrogenases [9]. Both enzymes contain two Rossmann-fold domains with similar orientation, thereby providing the six coordinating Cys residues [1]. Furthermore, both enzymes are capable of an ATP-dependent two electron reduction of acetylene and hydrazine [1]. Such a two-electron reduction, so far only observed in nitrogenases, is believed to start with acetylene binding to the metal cluster and transiently acting as a bridging ligand between two adjacent Fe-ions, which are 3.9 Å apart from each other [10]. A similar distance (3.74 Å) exists between Fe1 and Fe8 in DCCPCh, bridged by μ2-S, and facing the potential substrate channel [1]. Hence, this analogy with nitrogenases led to speculation that during catalysis the bridging sulfide site may dissociate leaving space for acetylene to bind Fe1 and Fe8, activating catalysis [1]. Further evidence that DCCP may function analogously to nitrogenase comes from inhibition experiments. As in nitrogenase, DCCPCh reduction of acetylene is reversibly inhibited by CO; this inhibition is speculated to occur when CO replaces the labile ligand [1].
For metal clusters, in which multiple iron centers are joined by sulfur atoms with a variety of geometric arrangements, an accurate description of the electronic structure is required. The electronic state of the transition metal cluster may have a fundamental influence, not only on the geometry, but also on the spectroscopic properties of the cluster. Different well-established experimental techniques exist for probing the various states of the catalytic cycle, as well as the physicochemical properties of FeS clusters. Among those, electron paramagnetic resonance (EPR) [10][11][12], infrared IR [13][14][15][16], resonance Raman (RR) [13][14][15][16][17], and Mössbauer [18] spectroscopy are powerful complementary techniques allowing the description of structural and electronic properties of the metal clusters [19]. EPR, based on microwave radiation, is well-suited to studying systems containing unpaired electrons, common for high-spin Fe centers [20]. Mössbauer spectroscopy relies on the resonant and recoil-free emission and absorption of gamma radiation by a chemical species. The gamma rays are accelerated through a range of velocities and the detected intensity is plotted as a function of the source velocity. When gamma rays are absorbed by the sample, a drop in measured intensity occurs. The Mössbauer spectrum thus reveals information about the chemical environment of the absorbing nuclei ( 57 Fe) [18].
The most common theoretical method for describing metal clusters containing multiple iron centers is the broken symmetry (BS) approach [21][22][23]. The need for a BS treatment is that electronic exchange interactions among multiple Fe centers are difficult to treat with density functional theory (DFT), which is limited to single-determinant wave functions. The BS methodology attempts to capture the multiple-determinant groundstate wave function by using a weighted average of pure spin states arising from a specific distribution of localized electrons. Alternatively, more recently, a type of multireference configuration interaction (CI) approach has been applied to treat coupled Fe spin states, relying on a finite complete active space (CAS) of interacting iron and sulfur orbitals [24,25]. To treat the [Fe2S2] 2+ cluster correctly, the minimum CAS of all 3d orbitals of iron atoms and 3p orbitals of the bridging sulfur atoms would lead to 16 [26]. Others have since applied similar methodology to study excited states of iron-sulfur clusters in various oxidation and spin states [27,28].
Therefore, for two Fe atoms, and analogously for four Fe atoms, the total S = 9/2 dimer spins are antiferromagnetically coupled to produce the resulting S = 0 (MS = 1) ground state ( Figure 2). Exchange interactions within the DCC cluster may be described by two predominant types of couplings: (1) within a single cubane, the two [Fe 2+ Fe 3+ ] dimers may couple (Jcubane); (2) between the two cubane clusters, two [Fe 2+ Fe 3+ ] dimers may couple via the bridging between μ2-S atom (Jcubane-cubane) ( Figure 1). This scheme is based on the assumption that the double-exchange interactions within the cubane are so strong that the mixed-valence [Fe 2+ Fe 3+ ] pair can be treated as a single spin, an assumption that is supported by spectroscopic studies [29,30]. A similar approach was successfully applied by Fiedler and Brunold in a "three-spin" model to study the exchange interactions between the cubane cluster [Fe4S4] and di-iron subunit [2Fe]H within the Fe-only hydrogenase [31]. This model assumes that the interaction between distantly separated [Fe 2+ Fe 3+ ] dimers is negligible [31].
According to the Heisenberg-Dirac-Van Vleck (HDVV) Hamiltonian, the exchange coupling constant J can be related to the energy through the Hamiltonian H = J(SA•SB) where SA and SB are the spin operators for centers A and B, respectively [21]. The Hamiltonian can be extended to account for double exchange via H = J(SA•SB) ± B(Stotal + 1/2), where B is the resonance delocalization parameter [23]. If the mixed-valence [Fe 2+ Fe 3+ ] pair is treated as a single spin, the delocalization parameter B is not treated explicitly [31].
Given these assumptions, and following the "three-spin" model of Fielder and Brunold [31], a "four-spin" model accounting for the coupling of four spin systems can be constructed for the DCC system ( Figure 3). The spin coupling constant, J, which describes the exchange interaction between two spins (SA, SB) can therefore be computed from the energies of the high spin (HS) and BS states via the following equation: Here, only the HS state is an eigenfunction of S 2 and the BS state, as it is not a pure spin state, is an eigenfunction of MS.
Multiple BS states can be conceived for eight Fe ions. Here, in addition to two high spin states (αααα-αααα and αααα-ββββ), 22 BS states were modeled in order to determine the electronic ground state of DCCPCh and the magnitude of Jcubane and Jcubane-cubane ( Figure 4). These 24 states are constructed to probe the wavefunction stability as a function of spin distribution within the DCC. Please note here that our spin notation (xxxx-xxxx) corresponds to Fe atoms (1234-5678), where Fe1 and Fe8 are the two bridging Fe atoms (Figure 1). Therefore, states with αxxx-xxxα have bridging Fe atoms of the same spin whereas αxxx-xxxβ have bridging Fe atoms of opposite spin. We have loosely adapted the BS numeric label syntax of each broken symmetry state, e.g., BS34(ααβ-αββα) from the literature on analogous [Fe4S4] clusters [32]. Our numeric labels are extended to two FeS clusters, where the first number refers to the spin states of Fe atoms 1234 and the second number refers to the spin states of Fe atoms 5678. The largest overall energy splitting is expected from the difference in energy between the high spin state E(HS = αααα-αααα) and a S = 0 broken symmetry state, e.g., E(BS = αβαβ-αβαβ) (Figure 3). Within a single cubane cluster, the splitting between the HS (αααα) and BS (e.g., βαβα) states is given by: EHS (Smax = 9/2 + 9/2 = 18/2) − EBS (MBS = 9/2 − 9/2 = 0) = 2 J (9/2)(9/2) = (81/2)J. (2) The energy coupling within a single cubane, Jcubane, of DCC can be obtained through an average over all EHS corresponding to one high-spin cubane (A-cluster with αααα) and the second cubane in all possible high and low spin BS states ( Figure 3).

Exchange Parameters Jcubane and Jcubane−cubane
The energies of all 24 quantum mechanically (QM) optimized spin states are listed in Table S1  cubane is found to be ten times smaller [31,35,36]. Upon ligand binding, J increases due to a change in the electronic distribution in the 2-Fe cluster. Specifically, prior to ligand binding at the distal Fe site, the unpaired electron is located at the distal Fe; following binding, the unpaired electron is located at the proximal Fe [35]. Similar arguments can be invoked to explain the magnitude of Jcubane−cubane associated with four coupled Fe-Fe dimers in DCC. Next, we turn our attention to hybrid quantum mechanical/molecular mechanical (QM/MM) geometry optimizations of DCCP.

QM/MM Optimized Geometries
QM/MM optimized geometries from the nine (S = 0) broken symmetry states were analyzed to deduce stabilizing structural features. The nine initial spin configurations were maintained throughout the whole geometry optimization, so that each final wave function describes each state's geometry and spin density properties uniquely. For the lowest energy state, BS25, the corresponding spin densities after optimization are reported in Table 1 and are shown in Figure 5.  The active sites of all nine QM/MM broken symmetry states show small root mean square deviation (RMSD) (0.1 Å) from the active site geometry in the crystal structure (PDB: 6ENO) [1]. The cubane clusters of all QM/MM optimized states preserve the starting crystal geometry (RMSD changes ranging between 0.10 Å and 0.17 Å) ( Table 2). In other studies of BS configurations of Fe4S4 atoms, similarly minimal structural differences were observed [31]. Small structural rearrangements in surrounding protein atoms are observed in all nine states: (1) a shift (0.7 Å) of one of the oxygen atoms from a nearby water molecule (w1) to bring it closer to the peptide backbone of Cys113; (2) a twist of the sidechain of Tyr376 away from its initial position, leading to a 0.7 Å movement of the hydroxyl oxygen atom to bring it closer to the positively charged sidechain of Arg312 (Arg312-Cz-OTyr376 is 3.31Å in the crystal structure and ~2.89 Å after QM/MM geometry optimization); (3) a small shift of Glu140 sidechain ( Figure 6, shown for lowest energy state, BS25).  The relative energies of each of the nine BS states indicate that the arrangement of unpaired electrons on the Fe atoms gives rise to energy differences (<5 kcal/mol) that are related to the electron configuration, in agreement with QM/MM calculations reported previously for H-cluster in Fe-only hydrogenase [37].
In each BS state, the electrons belonging to Fe ions bridging the μ2-S atom are spinpaired (i.e., Fe1 β and Fe8 α). In general, each symmetry state gives rise to a mixture of four short (type s = d < 2.25 Å) and four long (type l = d > 2.25 Å) Fe-S bonds (pattern shown for BS25 A-cluster in Table S2). Aside from BS35, the average bond lengths within the A-cluster are larger than those of the B-cluster ( Figure S2). The average cluster bond length ranges from 2.265 Å to 2.270 Å ( Figure S2).
Lower energy geometries can be linked to the spin densities' distributions on the different sub-clusters. The three lowest energy states, BS25(αβα-βαβα), BS15(ββααβαβα), and BS35(βααβ-βαβα) have in common the βαβα pattern of electron spins on Fe atoms 5,6,7,8 on the A-cluster, nearest to the residues Arg312, Glu140, Lys146 with chargebearing side chains. These three lowest energy states demonstrate a nearly identical pattern of bond lengths on the A-cluster ( Figure 7A), suggesting that this FeS cluster is highly sensitive to the charged protein environment and is the driving factor for overall DCC stability inside the protein scaffold.  In general, if two electrons of the same spin are located on iron atoms on the face of the cubane cluster, a lengthening of the sides of the face is observed. This phenomenon is further enhanced in the presence of an additional electrostatic environment, e.g., if amino acids with charge-bearing sidechains are nearby. Compared to the other BS states, the three lowest energy states, BS25, BS15, BS35, demonstrate a lengthened Fe6(α)-Fe8(α) distance (2.80 Å) and a shortened Fe5(β)-Fe7(β) distance (2.62 Å) ( Figure 8B). In other words, the iron atoms on the A-cluster face that is nearest to the charges localized on the Glu140 and Lys146 sidechains (-1 and +1, respectively) are farther apart than the iron atoms on the face toward the single charged sidechain of Arg312 (+1). The longest Fe-Fe separation (~2.85 Å) is observed for Fe5-Fe8 of states BS14, BS24, BS34, which have in common the (αββα) configuration on A-cluster iron atoms Fe5, Fe6, Fe7, Fe8. Iron atoms Fe5 and Fe8 carrying α-spin are located on the face pointing toward the B-cluster. In the absence of the protein environment, the Fe-Fe separations are largest (2.95 Å) for BS15 and BS16 ( Figure 7C), whose Fe1-Fe2 atoms both carry β spin and are located vis-à-vis the Acluster face containing Fe5 and Fe8, which, in these two BS states, also carry an unpaired electron of β spin. The iron atoms react to this high spin environment by seeking greater separation. Compared to vacuum ( Figure 8D-F), the protein environment reduces the scattering in Fe-S bond lengths in the DCC, and more systematic geometrical patterns emerge ( Figure 8A-C). Particularly for the BS24, BS25, BS26 states, Fe-S bond lengths are conserved in the B-cluster (circled purple squares in Figure 8A-C), reflecting the βαβαspin pattern that is common to these three states (Figure 4). Interestingly, the Fe8-S9 (Acluster Fe8 with bridging μ2-S atom) is consistently between 2.15-2.18 Å in both protein and vacuum environments. The Fe1-S9 bond, on the other hand, is lengthened considerably in the protein scaffold, reaching values as high as 2.2 Å (Figure 8A-C).
Taken together, these geometry analyses indicate that electrons of the same spin in a charged environment seek maximum separation. For the BS states investigated with the DCC system, these observations explain the overall more compact nature of the B-cluster, which is not in the direct vicinity of amino acids with charge-bearing sidechains. In general, cluster reduction by the addition of an electron would increase S-S atom repulsion and lead to overall cluster expansion [38]. For the DCC system, the energetic cost of reduction may therefore be smaller if an electron is added at the B-cluster as this cluster is not surrounded by amino acids carrying charged sidechains.

Broken Symmetry States
In addition to the HS1 state (αααα-αααα), following the BS-DFT approach to multicenter transition metal clusters, various parallel (ferromagnetic) or antiparallel (antiferromagnetic) alignments of the Fe iron site spins were considered, resulting in a total number of 22 different BS states defined following the numbering of the eight Fe atoms in Figure 1, i.e., 1234-5678 ( Figure 4). The semi-high spin HS2(αααα-ββββ) was also considered (see Introduction for a detailed description). In Figure 4, the grey shaded area includes BS states in which one cluster is high spin and the other cluster is low spin, either αxxx-xxxα or αxxx-xxxβ across bridging 2-S. The state BS52 (Figure 4) was examined only in the context of QM spin analyses. It is important to point out that, in the asymmetric protein matrix, each of these spin states is unique due to surrounding amino acids. In vacuum, six of the BS states (44,45,46,54,55,56) are expected to collapse into the six of corresponding patterns (37, 27, 17, 38, 28, 18, respectively), due to the symmetry of the DCC.
In each calculation, each of the sub-clusters is considered to be in the [Fe4S4] 2+ oxidized state. Hence, the formal number of ferrous (Fe 2+ , d 6 , S = 2) and ferric (Fe 3+ , d 5 , S = 5/2) iron sites is two and two, respectively, for a total Ms = 1 (Figure 2). For a single cubane cluster in an [Fe4S4] 2+ oxidized state with two each ferrous (Fe 2+ , d 6 , S = 2) and ferric (Fe 3+ , d 5 , S = 5/2) iron sites, two electron exchange phenomena should be considered [21][22][23]: (1) double exchange leads to ferromagnetic coupling of spins for a pair of Fe centers to yield two Fe-Fe dimers with S = 9/2 each; (2) according to Heisenberg exchange, the two S = 9/2 dimer spins are coupled antiferromagnetically to produce S = 0 ground state for a total Ms = 1 (Figure 2). Only the nine S = 0 BS states are considered in the discussion of DCC structural features (Figure 4, black text without gray shading).

Model for QM Geometry Optimization
Geometry optimizations of each of the spin states were carried out on an active site containing 65 atoms in vacuum. For these calculations the models considered only the atoms comprising the DCC, including four iron and four sulfur atoms for each sub-cluster, one μ2-S bridging the two cubane clusters and three cysteine thiolates binding each subsystem (Cys340, Cys308, Cys75 in the A-cluster; Cys113, Cys143, Cys373 in the Bcluster). The C-C covalent bonds connecting the C and the backbone of the six coordinating cysteine sidechains were cut and the valency of C atoms was satisfied by adding a proton to each of them. During geometry optimization of the DCC in vacuum, the C atoms were kept fixed to maintain their positions in the protein. A figure showing the 65 atoms included in the QM geometry optimization is included in the Supplementary Materials ( Figure S1). The geometry optimizations, carried out for each of the nine different spin configurations, were performed using the DFT functional BP86 [39,40], def2-TZVP basis set for Fe and S and 6-31g* for non-metals, respectively, combined with the resolution-of-identity (RI) technique as implemented in TURBOMOLE [41,42].

Model for QM/MM Geometry Optimization
In a second set of calculations, a hybrid QM/MM approach was used to investigate the effect of protein matrix on the DCC active site. For the QM/MM geometry optimizations, the system is divided into three computational layers (Figure 9). The first layer, treated at the QM level, consists of the 157 DCC atoms coordinated to subsystems (Cys340, Cys308, Cys75 in sub-cluster 1; Cys373, Cys113, Cys143 in sub-cluster 2) and atoms of sidechains (up to and including C) of the following amino acids: Cys94, Glu140, Lys146, Arg312, and Tyr376 ( Figure 9A). The side chains of Lys146 and Arg312 were considered to carry a positive charge, while Glu140 was modeled with a negative charge. Three water molecules (w1, w2, w3) were also included in this QM region. All atoms in the first computational layer were treated with DFT using the functional BP86 [39,40], def2-TZVP basis set for Fe and S and 6-31g* for non-metals, respectively, combined with the resolution-of -identity (RI) technique as implemented in TURBOMOLE [41,42]. The second layer of the model includes all protein and water atoms within a radius of 15 Å from the iron atom Fe8 located on the A-cluster. All atoms within this second layer were treated explicitly using the charmm22/CMAP force field [43] and allowed to relax during the geometry optimization cycles. The third layer includes all remaining protein, solvent, and salt ions extending 40 Å from Fe8 ( Figure 9C); the atoms in this third layer were held fixed during the optimization procedure.
The protein model for DCC was constructed using the 1.63 Å structure (6ENO.pdb) obtained through X-ray diffraction [1]. All enzyme atoms, including hydrogen atoms, were modeled with CHARMM [43], using the charmm22/CMAP force field [43,44]. Point charges of the DCC cluster atoms with coordinating Cys residues were computed from the electrostatic potential (ESP), as implemented in Gaussian16 [45] using DFT with the functional BP86 40], def2-TZVP basis set for Fe and S and 6-31g* for non-metals. The van der Waals parameters for the Fe atoms were assigned values reported in the literature for Fe in similar enzymatic systems [46].
The enzyme was next solvated in a box of dimension 105 Å × 119 Å × 121 Å with 42,393 explicit TIP3 water molecules [47], 32 sodium ions, and eight chlorine ions to neutralize charge. The positions of protein hydrogen atoms and solvent water atoms were optimized using a conjugate gradient energy minimization routine, as implemented in NAMD [48].
The charmm22/CMAP force field was used to describe all bonding and non-bonding interactions. Covalent bonds traversing the QM-MM border (between the Cα and C) were cut and hydrogen link-atoms were inserted to satisfy valency. According to established QM/MM protocols, the energetic coupling between the first and second layers (QM and mobile MM) was treated using electrostatic embedding with a charge-shift scheme as implemented in ChemShell [49,50]. For each of the nine BS states considered here, a geometry optimization was carried out.

Conclusions
Using a single-determinant DFT approach, the ground electronic state of spincoupled systems cannot be described, as the state is defined by a combination of multiple determinants. Within the framework of the BS formalism, however, these spin interactions can be modeled by localizing opposite spins of the wave function in different parts of the molecule. For enzymes with multi-metal centers, such as the novel multicenter transition metal cluster DCC, this approach provides us access to the lowest energy state as well as the exchange parameters between coupled spins. Here, 23 spin states were studied, and the calculated Jcubane and Jcubane−cubane show excellent agreement with values determined computationally and experimentally for analogous FeS systems. QM and QM/MM geometry optimizations were carried out for nine BS states. The protein scaffold was shown to play an essential role in determining the lowest energy state.
The catalytic activity of DCCPCh protein containing the DCC is comparable to the activity of enzymes such as the nitrogenases, which are able to reduce small molecules such as acetylene. In addition to this similar catalytic activity, the geometric symmetry of the P-cluster of nitrogenase suggests that the 8Fe cluster is formed from two closelyspaced 4Fe clusters [51]. Moreover, the precursor to the P-cluster is catalytically active and catalyzes the reduction of CO and CN to alkanes and alkenes [52]. These striking analogies suggest that DCCP and nitrogenase may have common ancestors. Here, the focus of the investigation was to provide the community with (1) an initial characterization of the DCC spin states, and (2) QM/MM optimized geometries of the DCCP active site based on a BS wavefunction. With this foundation, ongoing investigations are aimed at modeling substrate and inhibitor molecules, including acetylene and CO. Also, the protein scaffold's effect on the reduction potential in DCC will be characterized. DCCP may therefore serve as a model for a class of DCC-containing enzymes with unexplored reactivity and enzymatic potential.