Theoretical Study on Redox Potential Control of Iron-Sulfur Cluster by Hydrogen Bonds: A Possibility of Redox Potential Programming

The effect of hydrogen bonds around the active site of Anabaena [2Fe-2S] ferredoxin (Fd) on a vertical ionization potential of the reduced state (IP(red)) is examined based on the density functional theory (DFT) calculations. The results indicate that a single hydrogen bond increases the relative stability of the reduced state, and shifts IP(red) to a reductive side by 0.31–0.33 eV, regardless of the attached sulfur atoms. In addition, the IP(red) value can be changed by the number of hydrogen bonds around the active site. The results also suggest that the redox potential of [2Fe-2S] Fd is controlled by the number of hydrogen bonds because IP(red) is considered to be a major factor in the redox potential. Furthermore, there is a possibility that the redox potentials of artificial iron-sulfur clusters can be finely controlled by the number of the hydrogen bonds attached to the sulfur atoms of the cluster.


Introduction
In biosystems, the energy is sometimes transferred by electrons with various potential energies in a stepwise manner [1][2][3]. For this reason, organisms usually have a variety of electron transfer proteins. Iron-sulfur clusters are often found in the active sites of those electron transfer proteins in a wide range of organisms from bacterium to higher organisms [2]. It has been reported that the electron transfer proteins containing such ironsulfur clusters show a broad range of redox potentials although they have the same clusters in the active sites [3]. For example, [4Fe-4S] ferredoxin (Fd) and high potential iron-sulfur protein (HiPIP) are both the electron transfer proteins that have the same 4Fe-4S clusters in the active sites; nevertheless, their redox potentials are −400 ([Fe 4 S 4 (Cys) 4 ] 2−/3− ) and +360 ([Fe 4 S 4 (Cys) 4 ] 1−/2− ) mV, respectively [3].
A mechanism of how they control the redox potentials using the same 4Fe-4S clusters has attracted much attention in both bioscience and material science. Up to now, experimental and theoretical studies have proposed many hypotheses for the mechanism [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] such as a distortion of the cubane-type cluster [4], a hydrophobic atmosphere around the active site [6], a protein environment [7][8][9], and hydrogen bonds with solvent waters [10]. On the other hand, a development of X-ray crystallographic analyses has clarified a relationship between the number of hydrogen bonds around the active site and the redox potential [3,11,12]. There have also been some reports that support the importance of the hydrogen bonds by the theoretical calculations. For example, Gámiz-Hernández and co-workers reported that there is a relationship between the hydrogen bonds and the redox potential using optimized [Fe(SCH 2 CH 3 ) 4 ] structures under the solvent condition [13]. They clarified the importance of the number of hydrogen bonds for the redox potential, however, the effect of each hydrogen bond and its mechanism in the real active sites are still unknown. Usually the iron-sulfur clusters have several [N-H . . . S] hydrogen bonds between the peptide N-H bonds and sulfur atoms of the cluster as illustrated in Scheme 1. For instance, [4Fe-4S] Fd shows the lower redox potential has eight or nine hydrogen bonds [12], while HiPIP shows higher redox potential has five [6]. In general, the ironsulfur clusters are negatively charged if coordinating cysteines are included, therefore, a weak positive charge on the hydrogen atom by the charge polarization (-N δ− -H δ+ ) is expected to stabilize the active site. This hypothesis gives important insight to connect the protein structure and redox potential, however, there has not been direct evidence. Our group has studied the electronic structures of iron-sulfur complexes by theoretical calculations [17][18][19]. From the result of the 4Fe-4S cluster of HiPIP, we found that the hydrogen bonds around the cluster scarcely change orbital energy gaps, shapes of Kohn-Sham molecular orbitals, and magnetic interactions between Fe ions, although they stabilize orbital energies [19]. In addition, the calculated vertical ionization potential of the reduced state (IP(red)) suggests that the hydrogen bonds around the active site stabilize the reduced state by making a positive electrostatic environment. Because the ionization potential is one of the dominant factors of the redox potential of the protein [7,20], the result strongly suggests an importance of the hydrogen bonds for control of the redox potential. A detail of the relationship between the hydrogen bonds and the redox potential, however, has not been made clear yet. The iron-sulfur clusters are often utilized for various catalysts [21][22][23], so that the relationship contributes to a design guideline of such catalysts.
In this study, therefore, we aim to clarify the relationship between the hydrogen bonds around the active site in the iron-sulfur protein and the ionization potential by theoretical calculations in detail. As the first step, we focus on [2Fe-2S] ferredoxin ([2Fe-2S] Fd), which has a simple structure of the active site. In living organisms, one can find various proteins containing the 2Fe-2S clusters that work for electron transfer. For example, the plant-type and adrenodoxin-type Fds that are found in photosynthesis systems and in the metabolic system, respectively, are the well-known 2Fe-2S electron transfer proteins. The [2Fe-2S] Fds show various redox potentials from lower (ca. −450mV [24]) to higher (ca. −270 mV [25]) values, but nevertheless have the same cluster structure. In order to elucidate how the proteins control the redox potentials of the active site, the effect of the hydrogen bonds on the IP(red) values is examined by DFT calculations using the model structures of the [2Fe-2S] Fd.

Calculated Model Structures
A model structure of the [2Fe-2S] Fd is constructed from a result of the X-ray crystallographic analysis for an oxidized form of the plant-type Anabaena PCC7119 (PDB ID: 1QT9) [26]. Here, the hydrogen bonds around the 2Fe-2S active site are decided by two definitions as follows: (i) the distance between the hydrogen atom (H) and sulfur (S) (r H-S ) is shorter than 3.2Å; (ii) the angle among proton donor (D), H, and S (∠D-H· · · S) is larger than 139 • . Around the active site of Anabaena PCC7119, there are 11 hydrogen bonds including nine [N-H· · · S] hydrogen bonds from the peptide N-H bonds as well as two [O-H . . . S] bonds from side chains of Thr48 and Thr78 as illustrated in Figure 1A and as listed in Table 1. In order to focus on those 11 hydrogen bonds and to reduce a computational cost, we construct the minimal model including 13 amino acid residues (Phe39-Cys49, Thr78 and Cys79) that contain the 11 hydrogen bonds as illustrated in Figure 1A. All those amino acid residues are included in the model as a charge-neutral state. The N and C terminals of the amino acid residues in the model are capped by the hydrogen atoms. The Cartesian coordinates of all hydrogen atoms in the model including the hydrogen bonds are optimized by the DFT calculations, and other atoms are fixed to the X-ray structure.  Table 1. Structural parameters about the 11 hydrogen bonds around the active site included in the model structure (1 H ). The proton donor (D) is the amino acid residue that provides an N-H (or O-H) bond. The D-S distance is the distance between N and S (or O and S). The geometry of hydrogen atoms is optimized but the position of other atoms is fixed to the X-ray structure. Groups of the hydrogen bonds are classified by attaching sulfur atoms as shown in Figure 2A.  In order to clarify the effect of the hydrogen bonds, a naked cluster model in which the four cysteines attached to the Fe ions are replaced by methanethiolate (CH 3 S − ) ligands is also constructed as illustrated in Figure 1B. Here the model structures with/without the 11 hydrogen bonds are abbreviated as 1 H and 1 N , respectively. Furthermore, the effect of each hydrogen bond is also examined by changing amino acid residues in the model structure as explained below. The Cartesian coordinates of 1 H are summarized in Supplementary Materials Table S1.

Type of Hydrogen Bond Proton Donor (D) Sulfur (S) D-S Distance
The charge states of the active site for the reduced and oxidized states are [Fe(II,III) 2 S 2 (Cys) 4 ] 3− and [Fe(III,III) 2 S 2 (Cys) 4 ] 2− , respectively, including the coordinating Cysresidues [26]. As mentioned above, the attached amino acid residues (oligomer) are the charge neutral states, therefore, net charges of both 1 H and 1 N are the same. All Fe ions are a high spin state (Fe(II) s = 2, Fe(III) s = 5/2), and an antiferromagnetic coupling state is assumed as reported in [17] and [18].

Calculation of the Vertical Ionization Potential of the Reduced State (IP(red))
In experiments, a relative stability among several charge states is often discussed by using the redox potential (E redox ) that is defined by the difference in the free energies between the oxidized (ox) and reduced (red) states. In recent years, the redox potentials were calculated by the DFT calculations effectively [15,16,27,28]. As discussed later, the quantitative estimation of the redox potential is one of the important issues in the quantum chemistry. Some papers reported the good agreement of the computational results with the experimental ones [27,28], however, it is still a difficult task especially for large systems surrounded by the proteins. To avoid such difficulties, Noodleman and co-workers proposed a simplified method [7], where ∆E solv is a difference in solvation-free energies between the ox and red states. ∆SHE is a correction constant for the standard hydrogen electrode (SHE). We note that ∆SHE is −4.43 V instead of the 4.43 V [7]. IP(red) is the vertical ionization potential values of the red state calculated by where E(ox) and E(red) are calculated total energies of the ox and red states by the DFT method in the gas phase. It is reported that the IP(red) values are a major factor for the E redox in the 4Fe-4S systems, and dominant in some cases such as between [Fe 4 S 4 (SCH 3 ) 4 ] 0 and [Fe 4 S 4 (SCH 3 ) 4 ] 1− [7]. For this reason, we discussed the effect of the hydrogen bonds on the redox potential by using IP(red) under the assumption that the value is also key factor for the redox potential of the 2Fe-2S systems. The effect of solvation is discussed below.

Electronic Structures of the Oxidized and Reduced States
The calculated electronic structures for the ox and red states are confirmed by the atomic spin densities on Fe and S atoms as summarized in Table 2. Positive (up) and negative (down) spins on Fe ions indicate the antiferromagnetic spin coupling state. In the case of the ox state i.e., [Fe(III,III) 2 S 2 (Cys) 4 ] 2− , the values of the spin densities on Fe1 and Fe2 ions are smaller than the value expected from the high spin Fe(III) (five spins), indicating a small orbital overlap between Fe1 and Fe2 ions as reported in previous studies [17,18]. In addition, the beta spin densities found on 41Sγ and 46Sγ (or the alpha spin densities on 49Sγ and 79Sγ) suggest a partial alpha (or beta) electron transfer from the cysteine residues to the Fe(III) ions, which also decreases the spin densities on Fe(III) ion. In the case of the red state i.e., [Fe(II,III) 2 S 2 (Cys) 4 ] 3− , on the other hand, the spin density on Fe1 significantly decreases in comparison with the value of Fe2, indicating that Fe1 is reduced to Fe(II). In the red state, the negligible spin densities of 41Sγ and 46Sγ suggest that the partial electron transfer only occurs from cysteines to Fe(III) but not to Fe(II). In addition, the quantitatively equivalent spin density distributions between 1 H and 1 N indicate that the two models are the same electronic configuration.  Figure 1A.

Calculated Vertical Ionization Potential (IP(red)) of 1 H and 1 N
The calculated total energies of the ox and red states, and IP(red) values of 1 H and 1 N are summarized in Table 3. The calculated IP(red) values are negative, meaning that the ox state is more stable in comparison with the red state on both models. The absolute value of 1 H (1.65 eV) is, however, significantly smaller than the value of 1 N (5.06 eV). The difference in IP(red), which is about 3.4 eV, indicates that the red state is drastically stabilized by the amino acid residues around the active site compared to the ox state. In other words, the existence of the surrounding amino acids increases the relative stability of the red state. As mentioned above, the attached amino acid residues (oligomer) in 1 H are neutral charge, so that the net charge of the system is equivalent between 1 H and 1 N . In addition, the electronic configurations of 1 H and 1 N are also the same as shown in Table 2, therefore, the difference in IP(red) does not originate in the change of the charge or spin states, but in the effect of the surrounding amino acid residues. The results suggest two possibilities of mechanisms for change in IP(red) by the amino acid residues: (1) a local interaction by the hydrogen bonds, and (2) a whole electrostatic potential by residues. If the effect originates in the hydrogen bonds, IP(red) must be changed by single or several amino acids. In order to clarify the effect of the hydrogen bonds on IP(red) in detail, therefore, the effect of each hydrogen bond on IP(red) is separately examined.   Table 1, those hydrogen bonds are classified into three groups (Group 1-3) based on the attached sulfur atoms (Figure 2A). The hydrogen bonds in the group 1 and group 3 are bound to sulfur atoms of cysteines coordinating to Fe1 and Fe2, respectively, while those in group 2 are bound to bridging (inorganic) sulfur atoms. For each [N-H . . . S] hydrogen bond, we construct simple models, namely the peptide model, that consists of 1 N and a peptide bond moiety as illustrated in Figure 2B are no special hydrogen bonds that dominantly contribute to IP(red). In other words, all hydrogen bonds seem to contribute to IP(red) equivalently. In order to confirm a synergistic effect among hydrogen bonds on IP(red), a combination of hydrogen bonds within each group are examined as summarized in Table 4 Table 4. With an increase in the amino acid residues of the computational models, i.e., an increase of the number of hydrogen bonds, the absolute values of the calculated IP(red) become smaller. In Figure 3A, those calculated IP(red) values are plotted against the number of the hydrogen bonds. The result clearly indicates that the IP(red) values are in direct proportion to the number of the hydrogen bonds. It shows that the number of the hydrogen bonds contributes to IP(red), but the synergistic effect is not found. The slope of the plot is 0.33, indicating that a single hydrogen bond increases the relative stability of the red state and shifts IP(red) to the reductive side by 0.33 eV.

Effect of [O-H . . . S] Hydrogen Bond from the Side Chains on IP(red)
The effect of the [O-H . . . S] hydrogen bonds is also confirmed by using the models including the side chains. The models including the side chains, namely the amino acid model, are constructed by using 1 N and peptide bond moiety with the side chains as illustrated in Figure 2C. In the model, both N and C terminals are capped by hydrogen atoms. The geometry of the hydrogen atoms is optimized, but other atoms are fixed to 1 H . Some examples of the models are depicted in Supplementary Materials Figure S1. According to the strategy used for the peptide model, the IP(red) values are calculated by using the combination of the amino acid residues summarized in Table 5 Figure S1.
The calculated total energies of the ox and red states and IP(red) values are summarized in Table 5. The calculated IP(red) values are weakly perturbed by the side chains in comparison with the peptide models. The plot in Figure 3B, however, keeps a proportional relationship between the number of the hydrogen bonds and the IP(red) values even with the side chains. The slope (0.31) expresses that a single hydrogen bond shifts IP(red) to the reductive side by 0.31 eV, which is almost equivalent to the value of the peptide model. In contrast to the peptide model, the amino acid model shows the larger dispersion. The difference in the IP(red) values between the amino acid and peptide models is significant in Ala43, Thr48, Ser40, and Arg42. Those amino acids involve side chains that contribute the hydrogen bonds except for Ala43. In addition to the difference in a strength of the hydrogen bond between [N-H . . . S] and [O-H . . . S] [13], a part of the dispersion of the IP(red) values is considered to originate in the inter-amino acid interaction, although further analysis must be needed. On the other hand, a distribution of the H-S distance (r) and N-H-S angle (θ) does not affect the dispersion of IP(red) so much as shown in Supplementary Materials Figure S2. Table 5. Calculated total energies of the ox and red states and the IP(red) values of the amino acid models. The result under the solvent (water) effect for 1 H and 1 N is written in square parentheses.

Possible Mechanism for the Effect of Hydrogen Bond on IP(red)
In order to explain why the hydrogen bonds further stabilize the red state and shift IP(red) to the reductive side, an electron density difference map (∆ρ) between the ox and red states (ρ(ox)-ρ(red)) is depicted in Figure 4A, where ρ is the electron density distribution. The blue (negative) area is largely distributed especially on the active site because the red state has one more electron than the ox state. On the contrary, the hydrogen atoms concerning the hydrogen bonds around the active site are found to be colored by purple (positive). Figure 4B that focuses on the positive distribution around active sites clearly indicates that the hydrogen atoms concerning the hydrogen bonds become positive, indicating that the electron densities at those hydrogen positions are decreased in the reduced state. This suggests that the hydrogen atom concerned with the hydrogen bonds shows a stronger charge polarization (N δ− -H δ+ ) in the reduced state. As a result, the local positive electrostatic field by the hydrogen bonds is enhanced especially for the red state. The hydrogen bonds, therefore, are considered to shift IP(red) to the reductive side.

Materials and Methods
All calculations are performed by using the spin-unrestricted hybrid DFT method with the B3LYP functional set on the Gaussian program package [29]. Two different levels of basis sets i.e., BS-I (Fe, C, H: 6-31G, N, O: 6-31G* and S: 6-31+G*) and BS-II (Fe, C, N, O, H: 6-31G* and S: 6-31+G*) are used for the geometry optimization of hydrogen atoms, and following single-point energy calculations, respectively. Those calculations are carried out in the gas phase condition, unless otherwise specifically noted.

Conclusions
In this paper, the effect of the hydrogen bonds around the active site of Anabaena In addition, the value can be changed by the number of hydrogen bonds. Those results strongly suggest that the difference in IP(red) between the 1 N and 1 H (3.4 eV) originates in a summation of the local effects of the 11 hydrogen bonds included in the 1 H model rather than the effect by the whole amino acid residues (oligomer) included in the model. The results also suggest that the redox potential of [2Fe-2S] Fd can be accurately controlled by the number of hydrogen bonds. The information gives us some useful insights for mechanisms of the electron transfer as well as the control of the redox potential. For example, the X-ray structural analysis has revealed that the reduced structure of this protein (PDB ID: 1CZP) has one additional [N-H . . . S] hydrogen bond from Cys46 around the active site (12 hydrogen bonds in total) [26]. In other words, there is some possibility of controlling the electron transfer (giving and receiving) by the change in the number of the hydrogen bonds around the active site. Furthermore, the Fe1 site, which has more hydrogen bonds, is expected to be easily reduced, as supported by the calculated result in Table 1.
The importance of the hydrogen bonds around the active site is explained in this study, however, the obtained value (−1.65 eV is considered to be −1.65 V in the case of the one-electron transfer) is different from the experimental values even if ∆SHE is considered. The result suggests the strong contribution of environmental effect from the outer region i.e., the surrounding protein or water molecules. In fact, it is reported that ∆E solv becomes more important especially for more negatively charged [Fe 4 S 4 (SCH 3 ) 4 ] model systems [7]. As mentioned above, in addition, Gámiz-Hernández et al. estimated a change in the redox potential by the number of hydrogen bonds using the optimized [Fe(SCH 2 CH 3 ) 4 ] model clusters under the solvent condition. Our result is a few times larger than their conclusion (a change is 0.064-0.095 eV per hydrogen bond) [13]. In order to evaluate the importance of the surrounding environment, we calculate IP(red) of 1 H and 1 N under a solvent (water) condition with the polarizable continuum model (PCM) using the integral equation formalism variant (IEFPCM). The calculated total energies and IP(red) are summarized in Table 5. The calculated IP(red) values (3.32 and 2.89 eV for 1 H and 1 N , respectively) is drastically changed by the solvent water. It is reported that the solvent effect considered by the IEFPCM sometimes overestimates the electrostatic effect [19], however the result strongly suggests the importance of the environmental effect from the outer region on the redox potential, too. In addition, the differences between 1 H and 1 N (0.43 eV) are decreased a lot by the solvent effect, indicating that the dielectric field of the polar solvent stabilized the red state as a substitute for the positive electrostatic field from the hydrogen bonds in 1 N . Conversely, the effect of hydrogen bonds around the active sites is expected to stand out especially in a hydrophobic atmosphere condition such as HiPIP. In fact, the difference between 1 H and 1 N (3.4 eV) can be also said to be almost comparable to the difference between the gas phase and solvent of 1 H (4.97 eV), suggesting the strong contribution of the local hydrogen bonds around the active site.
Considering the functionality of the protein, in addition, the electron transfer (giving and receiving) must be controlled accurately. From the above results, a local rearrangement of the hydrogen bonds is considered to be suitable for the accurate control of the redox for the electron transfer, while the surrounding protein and solvents are suited to shift the redox potential significantly. In order to clarify roles of the hydrogen bonds and protein, ONIOM calculations including a whole protein and water molecules is necessary, but we focus on the effect of the hydrogen bonds and do not discuss further about the effect from the outer region here. In addition to the consideration of the surrounding protein, the quantitative estimation of the redox potential is also required for the further elucidation of the mechanism. There are many reports about the theoretical calculation of the redox potentials [15,16,27,28,30]. Those approaches must also be considered for the quantitative analyses and further discussions.
There have been, on the other hand, many reports about the synthesized iron-sulfur cluster complexes that aim to add catalytic functionality and so on [22,23]. The above results suggest that the redox potentials of such artificial iron-sulfur clusters can be finely controlled by the number of the hydrogen bonds attached to the sulfur atoms of the clusters. The tune or control of the redox potential by the number of hydrogen bonds can be an effective designing guideline for the iron-sulfur clusters, suggesting a possibility of potential redox programming by the surrounding ligands.