Computational Approach for Probing Redox Potential for Iron-Sulfur Clusters in Photosystem I

Simple Summary Many biological systems contain iron–sulfur clusters, which are typically found as components of electron transport proteins. Continuum electrostatic calculations were used to investigate the effect of protein environment on the redox properties of the three iron–sulfur clusters in the cyanobacterial photosystem I. Our results show a good correlation between the estimated and the measured reduction potential. Moreover, the results indicate that the low potential of FX is shown to be due to the interactions with the surrounding residues and ligating sulfurs. Our results will help in understanding the electron transfer reaction in photosystem I. Abstract Photosystem I is a light-driven electron transfer device. Available X-ray crystal structure from Thermosynechococcus elongatus showed that electron transfer pathways consist of two nearly symmetric branches of cofactors converging at the first iron–sulfur cluster FX, which is followed by two terminal iron–sulfur clusters FA and FB. Experiments have shown that FX has lower oxidation potential than FA and FB, which facilitates the electron transfer reaction. Here, we use density functional theory and Multi-Conformer Continuum Electrostatics to explain the differences in the midpoint Em potentials of the FX, FA and FB clusters. Our calculations show that FX has the lowest oxidation potential compared to FA and FB due to strong pairwise electrostatic interactions with surrounding residues. These interactions are shown to be dominated by the bridging sulfurs and cysteine ligands, which may be attributed to the shorter average bond distances between the oxidized Fe ion and ligating sulfurs for FX compared to FA and FB. Moreover, the electrostatic repulsion between the 4Fe-4S clusters and the positive potential of the backbone atoms is lowest for FX compared to both FA and FB. These results agree with the experimental measurements from the redox titrations of low-temperature EPR signals and of room temperature recombination kinetics.


Introduction
The photosynthesis process is the process that guarantees the existence of our life. In photosynthesis, solar energy is harvested by pigments associated with the photosynthetic machinery and stored as energy-rich compounds [1]. Initial energy conversion reactions take place in special protein complexes known as Type I and Type II reaction centers [2], which are classified according to the type of terminal electron acceptor used, iron-sulfur clusters (Fe-S) and mobile quinine for type I and type II, respectively [2][3][4][5][6][7]. Photosystem I (PS I) is the Type I reaction center found in the thylakoid membranes of chloroplasts and I (PS I) is the Type I reaction center found in the thylakoid membranes of chloroplasts and cyanobacteria [6,8]. PS I is a very interesting electron transfer machine that converts solar energy to a reducing power with a quantum yield close to 1 [9][10][11]. It, mainly, mediates the transfer of electrons from either cytochrome c6 or plastocyanin to the terminal electron acceptor at its stromal side through a series of redox reactions along electron transfer chains. The crystal structure of a trimeric cyanobacterial PS I is resolved at an atomic resolution of 2.5 Å [12], where each monomer consists of about 12 polypeptide chains (PSaA-PsaX) (Figure 1a,b). There are three highly conserved chains in PS I: PsaA, PsaB and PsaC [13]. The first two chains form the heterodimeric core, which noncovalently bonds most of the antenna pigments, redox cofactors employed in the electron transfer chains (ETCs) and the inter-polypeptide iron-sulfur cluster FX [14,15]. PsaC comprises two iron-sulfur clusters FA and FB, and it forms, with PsaE and PsaD, the stromal hump providing a docking site for soluble protein ferredoxin [16,17] (Figure 1a). Cofactors employed in the ETCs are a chlorophyll (a) dimer P700, two pair of chlorophyll a molecules A/A0 and two phylloquinones A1. These cofactors are arranged in two nearly symmetric branches A and B, from P700 at the luminal side to FX at the PsaA and PsaB interface followed by the two terminal iron-sulfur clusters FA and FB ( Figure 1c) [8,18,19].  [12]: Twelve protein subunits in the polypeptide structure of cyanobacterial PS I monomer viewed perpendicular to the plane of the thylakoid membranes. (a) Front view; (b) back view; (c) electron transfer chains (ETCs) in PS I, including primary electron donor P700 (Chl a dimer), primary electron acceptors A/A0 (Chl a molecules), secondary electron acceptor A1 (phylloquinone molecule PQN), tertiary electron acceptor X (FX) and terminal electron acceptor A (FA) and B (FB) [8].
Upon photoexcitation of a primary electron donor P700, an electron will transfer to the primary electron acceptor A/A0, within ~100 fs [20], followed by an electron transfer to the phylloquinone molecule within 20-50 ps [19]. Then, the electron is transferred, sequentially, to the three iron-sulfur clusters, FX, FA and FB, within ~1.2 μs [19]. It was Figure 1. PDB code 1JB0 [12]: Twelve protein subunits in the polypeptide structure of cyanobacterial PS I monomer viewed perpendicular to the plane of the thylakoid membranes. (a) Front view; (b) back view; (c) electron transfer chains (ETCs) in PS I, including primary electron donor P700 (Chl a dimer), primary electron acceptors A/A 0 (Chl a molecules), secondary electron acceptor A 1 (phylloquinone molecule PQN), tertiary electron acceptor X (F X ) and terminal electron acceptor A (F A ) and B (F B ) [8].
Upon photoexcitation of a primary electron donor P700, an electron will transfer to the primary electron acceptor A/A 0 , within~100 fs [20], followed by an electron transfer to the phylloquinone molecule within 20-50 ps [19]. Then, the electron is transferred, sequentially, to the three iron-sulfur clusters, F X , F A and F B , within~1.2 µs [19]. It was shown that the reduced F B will directly reduce a soluble protein ferredoxin (Fd), which in turn will reduce the NADP+ to NADPH in the ferredoxin-NADP+ reductase complex (FNR) [3][4][5][6][7][21][22][23]. Knowing the redox potentials of these cofactors is crucial for understanding the primary photosynthetic processes. However, the complexity of the PS I protein complex and the electrostatic nature of interactions between charged groups and among redox centers make it difficult to assign the measured signals to a specific redox-active center. Thus, computational methods could be a complementary technique for the characterization of redox reactions. The three iron-sulfur clusters in PS I are 4Fe-4S clusters; 4Fe-4S is a distorted cube of four iron atoms linked by four bridging sulfur atoms and ligated by four cysteine ligands [8]. The PsaC polypeptide chain provides the cysteine ligands to both clusters F A and F B : C53, C50, C20 and C47 for F A and C10, C57, C13 and C16 for F B . The F X cluster is ligated by four cysteines: two from the PsaA chain (C578 and C587) and two from the PsaB chain (C565 and C574). They are mainly distinguished by their low-temperature EPR spectrum [24,25]. In PS I, F X , F A and F B are known as low-potential [4Fe-4S] clusters that employ the 2 + /1 + redox couple [26][27][28]. In its oxidized state, a lowpotential [4Fe-4S] cluster has two ferric and two ferrous Fe atoms and possesses a total spin S = 0. In its reduced state, it has one ferric and three ferrous Fe atoms with total spin S = 1/2. This is due to the paramagnetic pairing between an equal-valence pair Fe +2 −Fe +2 and a mixed-valence pair Fe +2.5 −Fe +2.5 [8]. In PS I, the redox potentials of 4Fe-4S clusters vary in a wide range from −730 to −440 mV [19]. Low-temperature electron paramagnetic resonance (EPR) spectroscopy studies showed that the midpoint potentials are −705 ± 15, −530 and −580 mV for F X , F A and F B clusters, respectively [8,19]. However, other studies suggested that the midpoint potentials of these clusters would be positively shifted at room temperature [29][30][31][32]. Redox potentials of iron-sulfur clusters in PS I were calculated previously by Torres et al. [33]. They reported the E m values for F X , F A and F B to be −980, −510, and −710 mV, respectively, where the E sol m was obtained by correcting the ionization potential calculated by gas-phase DFT with the solvation effects and referencing the calculated potential to the standard hydrogen electrode (∆SHE = −4.5 eV). Torres et al. employed a model with three dielectric regions, the continuum solvent (ε wat = 80), the protein (ε prot = 4) and ε = 1 for the redox site to reflect the little screening effect of protein due to hydrogen bonding in the vicinity of the clusters. In their paper, Ptushenko et al. [34] argued the implausibility of the proposed three dielectric regions model by Torres due to the overestimation of the amide field in the vicinity of clusters, which lead to a negatively deviated midpoint potential from experimental values by 275 to 330 mV for F X and 130 to 245 mV for F B . In the work of Ptushenko and coworkers [34], they calculated midpoint potentials for all redox cofactors in PS I, including the three [4Fe-4S] clusters. Their reported values are −654, −481 and −585 mV for F X , F A and F B, respectively. In their calculations, they employed the semicontinuum electrostatic model, where two dielectric constants for proteins were used, the optical dielectric constant (ε o = 2.5) for pre-existing permanent charges and a static dielectric constant (ε s = 4) for charges formed due to the formation of ions in protein upon ionization reaction. In their calculations, Torres and Ptushenko included all protein subunits and other prosthetic groups in the PS I complex.
Here, we report the calculated relative midpoint potential of [4Fe-4S] clusters F X , F A and F B , using Multi-Conformer Continuum Electrostatics (MCCE) [35][36][37]. Although this model is based on classical electrostatics and does not consider the quantum effects such as spin-spin and spin-orbit interactions, it successfully produced the experimental pK a and E m values for oxo-manganese complexes and Mn and Fe superoxide dismutase [38,39]. Moreover, we provide an insight into the conformational changes and the interactions that induce the differences in the redox potential of the three [4Fe-4S] clusters from the classical electrostatics perspective and their implication on the electron transfer reaction.

Structural Model
PS I crystal structure of Thermosynechococcus elongatus (PDB code: 1JB0 [12]), at resolution 2.5 Å, was used as the basis for all calculations. All crystallographic water Biology 2022, 11, 362 4 of 14 molecules and lipids were removed prior to our calculations. There are 49 cofactors and 35 amino acids with missing atoms. Moreover, atomic coordinates of~91 residues were not reported during experiments. Prior to MC calculations, PDBfixer (Python-based software) [40] was used to add missed inter-polypeptide residues only and to add all other missed atoms in residues and cofactors.
The fixed structure was solvated by water molecules (more than 200,000) to fill a cubic box (~19 × 19 × 19 Å). Sodium counterions were added by replacing water molecules to provide a neutral simulation box. The resulting complex comprised approximately 700,000 atoms. Atomic interactions for the standard amino acids and cofactors (except 4Fe-4S clusters) in PS I were based on the developed AMBER molecular potentials for PS II [41]. The parameterization of the iron-sulfur clusters was based on the work of Carvalho et al. [42]. Energy minimization calculation with the steepest descent method followed by conjugate gradient energy minimization was conducted using the Gromacs program. The minimized structure was further equilibrated for 100 ps in the NVT ensemble followed by another NPT equilibration for 100 ps. Then, a short MD simulation (10 ns) was performed, where trajectories were collected every 10 ps. The final coordinates retrieved from MD simulation were used for further calculations. On the other hand, the geometry for [4Fe-4S] clusters F X , F A and F B , surrounded by~10 Å nearby residues, as shown in Figure 2, were extracted from the crystal structure and were optimized at the DFT/B3LYP level of theory, with LANL2DZ basis sets [43] for Fe metal centers and 6-31G* basis set for other atoms, using the Gaussian09 package [44]. The [4Fe-4S] core was set to the reduced state with total spin S = 1 2 using the broken symmetry wavefunction [45]. Each DFT optimized structure was docked into the complete PS I structure retrieved from MD simulations such that the coordinates of side chains that are not included in the DFT clusters and all backbone atoms remain the same as in MD-obtained structures. This final complex was used as an initial state for MC simulations (structures are available in Files S1-S3).

Structural Model
PS I crystal structure of Thermosynechococcus elongatus (PDB code: 1JB0 resolution 2.5 Å, was used as the basis for all calculations. All crystallographic wa ecules and lipids were removed prior to our calculations. There are 49 cofactors amino acids with missing atoms. Moreover, atomic coordinates of ~91 residues w reported during experiments. Prior to MC calculations, PDBfixer (Python-based so [40] was used to add missed inter-polypeptide residues only and to add all other atoms in residues and cofactors. The fixed structure was solvated by water molecules (more than 200,000) to f bic box (~19 × 19 × 19 Å). Sodium counterions were added by replacing water mo to provide a neutral simulation box. The resulting complex comprised approx 700,000 atoms. Atomic interactions for the standard amino acids and cofactors (exc 4S clusters) in PS I were based on the developed AMBER molecular potentials f [41]. The parameterization of the iron-sulfur clusters was based on the work of C et al. [42]. Energy minimization calculation with the steepest descent method follo conjugate gradient energy minimization was conducted using the Gromacs progra minimized structure was further equilibrated for 100 ps in the NVT ensemble follo another NPT equilibration for 100 ps. Then, a short MD simulation (10 ns) was per where trajectories were collected every 10 ps. The final coordinates retrieved fr simulation were used for further calculations. On the other hand, the geometry f 4S] clusters FX, FA and FB, surrounded by ~10 Å nearby residues, as shown in F were extracted from the crystal structure and were optimized at the DFT/B3LYP theory, with LANL2DZ basis sets [43] for Fe metal centers and 6-31G* basis set fo atoms, using the Gaussian09 package [44]. The [4Fe-4S] core was set to the reduc with total spin S = ½ using the broken symmetry wavefunction [45]. Each DFT op structure was docked into the complete PS I structure retrieved from MD simulatio that the coordinates of side chains that are not included in the DFT clusters and a bone atoms remain the same as in MD-obtained structures. This final complex w as an initial state for MC simulations (structures are available in File S1, File S2 a S3).

Multi-Conformer Continuum Electrostatics (MCCE) Calculations
MCCE generates different conformers for all amino acid residues and co These conformers undergo a preselection process, which discards conformers tha rience vdW clashes [36]. All crystallographic water molecules and solvated i stripped off and replaced with a continuum dielectric medium. The electrostatic p

Multi-Conformer Continuum Electrostatics (MCCE) Calculations
MCCE generates different conformers for all amino acid residues and cofactors. These conformers undergo a preselection process, which discards conformers that experience vdW clashes [36]. All crystallographic water molecules and solvated ions are stripped off and replaced with a continuum dielectric medium. The electrostatic potential of the protein is calculated by solving Poisson-Boltzmann equation [46] using DelPhi [47]. In this calculation, the surrounding solvent (water) was assigned a dielectric constant of ε = 80, and ε = 4 was assigned for protein [48]. The partial charges and radii used for amino acids in MCCE calculations are taken from the PARSE charges [49]. The probe radius for placing water is 1.4 Å and 0.15 M salt concentration is used. For 4Fe-4S clusters, each Fe ion, bridging sulfurs S and each ligand as separate fragments with an integer charge interact with each other through electrostatic and Lennard-Jones potentials [39].
The Fe atoms have formal charges of +2 or +3, while each bridging sulfur atom has a charge of −2.
For each conformer i, DelPhi calculates different energy terms, the polar interaction energy (∆G pol,i ), desolvation energy ∆∆G rxn,i and pairwise electrostatic and Lennard-Jones interactions with other conformers j (∆G ij ). For M conformers, the ∆∆G rxn,i and ∆G pol,i energy terms will be collected into two matrices with M rows while the ∆G ij energy term will be collected into the M × M matrix [37]. A single protein microstate x is defined by choosing one conformer for each cofactor and residue. Therefore, the number of possible microstates of the system is very high. As a final step, MCCE uses Monte Carlo sampling to compute the probability of occurrence for each conformer in the Boltzmann distribution for a given pH and electron concentration (E h ) [37,50].
The total energy of each microstate G x with M conformers is the sum of electrostatic and non-electrostatic energies, and it is computed according to the following equation [51,52]: where δ x,i is equal to 0 if microstate x lacks conformer i and 1 otherwise; m i takes the values 0, 1 and −1 for neutral, base and acid conformers, respectively; n i is the number of electrons transferred during redox reactions; pK a,sol,i and E m,sol,i are the reference pK a and E m for i-th group in the reference dielectric medium (e.g., water); F is the Faraday constant; K b is the Boltzmann constant; T is temperature (298 K in our calculations); ∆∆G rxn,i is the desolvation energy of moving conformer i from solution to its position in the protein; ∆G ij is the pairwise interaction between different conformers i and j; and ∆G pol,i is the pairwise interaction of conformer i with other groups with zero conformational degrees of freedom (e.g., backbone atoms).
The reference solution E m,sol for Fe ions is obtained according to the thermodynamic cycles shown in Figure 3. The solution redox potential, E m,sol is fixed at −170 mV, a value that reproduced the redox chemistry of F X (−715 mV vs. SHE) [39].  is the midpoint potential in reference medium and is the midpoint potential in situ calculated by MCCE.

Mean-Field Energy (MFE) Analysis
MCCE determines the in situ midpoint potential of the redox centers as shifted by the protein environment. This shift is due to the loss in the reaction field energy ∆∆ and other electrostatic interactions. Mean-field energy analysis (MFE) allows decomposition of these energetic terms to determine what factors yield the reported midpoint potentials in protein, Equation (2)

Mean-Field Energy (MFE) Analysis
MCCE determines the in situ midpoint potential E m of the redox centers as shifted by the protein environment. This shift is due to the loss in the reaction field energy ∆∆G rxn and other electrostatic interactions. Mean-field energy analysis (MFE) allows decomposition of  (2) where ∆G bkbn is the electrostatic and non-electrostatic interactions of the redox cofactor with the backbone atoms of protein and ∆G MFE res is the mean-field electrostatic interaction between the redox cofactor and the average occupancy of conformers of all other residues in the protein in the Boltzmann distribution at each E h [53]. Other terms are the same as shown in Equation (1).

Results and Discussion
Molecular structures for [4Fe-4S] clusters in PS I were investigated by extended X-ray absorption fine structure (EXAFS), which revealed two peaks at~2.27 and~2.7 Å, which are attributed to the backscattering from sulfur and iron atoms, respectively [54][55][56][57]. The results of geometry optimization of three extracted structures with total spin S = 1 2 and with [4Fe-4S] in their reduced state are reported in Table 1a. Our calculated Fe-S (bridging sulfur atoms), Fe-SG (organic sulfur atoms) and Fe-Fe bond distances are shown, generally, to be longer than the XRD [12] and EXAFS reported distances (Table 1).  Bold values are the distances between the second oxidized Fe ion and the four sulfur ligands (three bridging sulfurs and one sulfur from cysteine), while Avg. is the average over these distances for each 4Fe-4S cluster (see Figure 4).

The Midpoint Potentials (E m ) of F X , F B and F A
In our calculations, we considered the oxidation potential of the second oxidized Fe ion as the oxidation potential of the cluster from [4Fe-4S] +1 to [4Fe-4S] +2 . The measured E m values of F X , F A and F B are reported in Table 2. For F X , E m is −715 mV, which is 0-65 mV more negative than experimental values [8,58,59]. The E m values calculated for F A and F B are −355 and −346 mV, respectively, and are positively shifted by more than 80 mV compared to experimentally determined values [19,60,61]. Our results are shown to agree with the experimental values within the error range of the method [37,62,63]. To better understand the effect of ligands and other residues in the model structures on the calculated E m , mean-field energy (MFE) analysis was performed for each 4Fe-4S cluster at its calculated E m to determine the different factors contributing to the stabilization of the ionization state of clusters in the protein (Equation (2)). Results from the MFE analysis are reported in Table 3. The calculated E m values are shown to deviate from the reference value by about 13 kcal/mol for F X , while for both F A and F B it is shown that the values are only shifted by about 4 kcal/mol, which may be attributed to the degree of burial of each 4Fe-4S cluster. Our calculations demonstrated that the desolvation energy ∆∆G rxn was always positive and unfavorable. It destabilizes the ionization state by~80 kcal/mol for F X and by~76 kcal/mol for F A and F B . This unfavorable interaction is, nearly, compensated by the electrostatic interactions with the surrounding residues ∆G resd for F X , F A and F B to be −104, −97 and −92 kcal/mol, respectively. Table 2. Calculated midpoint potential for redox couples +2/+1 (in units of mV).

Cal. E m
Exp. E m The bold is the E m value used for determining the solution reference redox potential; i [24,64], j [61], k [58], l [31], m [65]. Moreover, interactions with the backbone ∆G bkbn disfavor the oxidized form of the Fe. For F A and F B, this effect is shown to be~1.13-fold more than that in F X . By further breaking down the contribution from different residues and ligands (see Table 4), it is shown that stabilization of ionization state of 4Fe-4S clusters is mainly controlled by the classical electrostatic interactions between Fe ions and both bridging sulfurs and cysteine ligands. Oppositely, the electrostatic interaction with positively charged residues and other Fe ions is shown to destabilize the oxidized form of Fe ion.  Furthermore, to investigate the role of different interactions in altering the redox potential of 4Fe-4S clusters, interaction energies between bridging sulfurs and the average occupancy of conformers of all other residues in the protein complex were calculated at the E m of each cluster. The total interaction energy of S3 in the F B cluster with other residues Biology 2022, 11, 362 9 of 14 and cofactors in PS I protein is shown to be the strongest interaction compared to both S3 in F X and F A and compared to other sulfur atoms in the same cluster (Figure 5b). By further breaking down the contribution from different residues, it is shown that interaction energies are favorable between inorganic sulfurs and the charged and polar amino acids in the vicinity of the clusters. The electrostatic interaction between Lys-C51 and S2 of F A clusters is shown to be the strongest interaction compared to other bridging sulfurs in the same cluster or in the other two clusters with interaction energy~17 kcal/mol (see Table 5 and Figure 5a). The bridging sulfurs in the inter-polypeptide cluster F X are hydrogen-bonded to the charged residues, mainly Lys and Arg residues, in the vicinity of the cluster with a total energy of~−152.5 kcal/mol. However, the stromal clusters exhibited less interaction between bridging sulfurs and polar and charged residues, where the total interaction energies with inorganic sulfurs in F A and F B are −126 and −30 kcal/mol, respectively. Bold values are the distances between the second oxidized Fe ion and the four sulfur ligands (three bridging sulfurs and one sulfur from cysteine), while Avg. is the average over these distances for each 4Fe-4S cluster (see Figure 4).
* total interaction energy of bridging sulfurs with each residue in the vicinity of the cluster.

Conclusions
Based on our calculations, the low potential of F X is shown to be due to the backbone and residue sidechain contributions. Moreover, the distances between ligating sulfurs and the second oxidized Fe ion were found to be, on average,~2.3 Å for F X and~2.5 Å for F A and F B . This could explain the stronger effect of sulfurs in F X for shifting the redox potential. Moreover, the bridging sulfurs in the inter-polypeptide cluster F X interact with the charged residues in the cluster's proximity, namely Lys and Arg residues, with a total interaction energy of −152.5 kcal/mol. The stromal clusters, on the other hand, showed reduced contact between bridging sulfurs and polar and charged residues, with total interaction energies with inorganic sulfurs of −126 and −30 kcal/mol in F A and F B , respectively. When compared to S3 in F X and F A , as well as other sulfur atoms in the same cluster, the overall interaction energy of S3 in the F B cluster with other residues and cofactors in PS I protein is demonstrated to be the greatest. Finally, when compared to other bridging sulfurs in the same cluster or in the other two clusters, the electrostatic interaction between Lys-C51 and S2 of F A clusters is revealed to be the strongest interaction, with an interaction energy of 17 kcal/mol. In this study, our calculations demonstrated how each 4Fe-4S cluster is exposed to a different protein electrostatic environment, which tunes the redox properties of each 4Fe-4S cluster.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biology11030362/s1, File S1: The optimized structure of F X and surrounding residues within~10 Å; File S2: The optimized structure of F A and surrounding residues within~10 Å; File S3: The optimized structure of F B and surrounding residues within~10 Å

Data Availability Statement:
The structural models that were used are available in the SI online.