Redox potentials along the redox-active low-barrier H-bonds in electron transfer pathways

Abstract

Low-barrier H-bonds form when the pK_a values of the H-bond donor and acceptor moieties are nearly equal. Here, we report redox potential (E_m) values along two redox-active low-barrier H-bonds in the water-oxidizing enzyme photosystem II (PSII), using a quantum mechanical/molecular mechanical approach. The low-barrier H-bond between D1-Tyr161 (TyrZ) and D1-His190 is located in the middle of the electron transfer pathway. When the proton is at D1-His190, E_m(TyrZ) is the lowest and can serve as an electron donor to the oxidized chlorophyll P_D1˙⁺. E_m(TyrZ) and E_m(D1-His190) are equal, and the TyrZ⋯D1-His190 pair serves as an electron acceptor to Mn₄CaO₅ when the proton is at TyrZ. In the low-barrier H-bond between D1-His215 and plastoquinone Q_B, located at the terminus of the electron transfer pathway, the driving force of electron transfer and electronic coupling between Q_A and Q_B are maximized when the proton arrives at Q_B. It seems likely that local proton transfer along redox-active low-barrier H-bonds can alter the driving force or electronic coupling for electron transfer.

---

The driving force of water oxidation in photosystem II (PSII) is provided by light-induced charge separation in the reaction center. In PSII, the reaction center has a pair of chlorophylls (P_D1/P_D2), accessory chlorophylls (Chl_D1/Chl_D2), pheophytins (Pheo_D1/Pheo_D2), and plastoquinones (Q_A/Q_B) in the heterodimeric D1/D2 protein subunit pairs (Fig. 1a).¹ Electronic excitation of Chl_D1 leads to the formation of a charge-separated state [P_D1P_D2]˙⁺Pheo_D1˙⁻;² electron transfer subsequently occurs via Q_A to Q_B. Q_B reduction is linked to Q_B protonation. After the second electron transfer, doubly protonated Q_BH₂ is released from the binding site toward the bulk water region: the quinone pool. In [P_D1P_D2]˙⁺, the cationic state is more populated on P_D1 than P_D2 (P_D1˙⁺ : P_D2˙⁺ = 80 : 20^3–6). P_D1˙⁺ has a significantly high redox potential (E_m) for one-electron oxidation (>1100 mV^6–9), which makes P_D1˙⁺ abstract electrons from the substrate water molecules at the Mn₄CaO₅ cluster (700–800 mV¹⁰) via redox-active D1-Tyr161 (TyrZ; Fig. 1b).

Fig. 1 (a) Electron transfer pathways in PSII (PDB code: 3ARC). Low-barrier H-bond pairs are labeled in red. Black arrows indicate electron transfer between cofactors. The dotted line indicates pseudo-C₂ axes. (b) Low-barrier H-bond between TyrZ and D1-His190 near Mn₄CaO₅ cluster. (c) Typical potential energy profile of standard H-bonds. O_A and O_D stand for the acceptor and donor oxygen moieties, respectively. (d) Typical potential energy profile of low-barrier H-bonds (LBHBs). In low-barrier H-bonds, O_A and O_D cannot be distinguishable due to the same pK_a values.

TyrZ must form a hydrogen bond (H-bond) with D1-His190 to serve as a redox-active cofactor. Notably, TyrZ and D1-His190 form a significantly short (2.46 Ź¹) low-barrier H-bond.^12,13 Low-barrier H-bonds can form only when the pK_a values of the H-bond donor and acceptor moieties are nearly equal.^14,15 The shape of the potential energy curve of a low-barrier H-bond is symmetric, while the curve of a standard H-bond is asymmetric because pK_a(donor) > pK_a(acceptor)¹⁶ (Fig. 1c and d). The pK_a value of tyrosine (∼10) is higher than that of histidine (∼7) in water. In the PSII protein environment, a cluster of water molecules (whose positions are fixed by the components of the Mn₄CaO₅ cluster) donates a stable H-bond to the phenolic O site of TyrZ and decreases pK_a(TyrZ) to the level of pK_a(D1-His190), forming a low-barrier H-bond.^12,13

In PSII, low-barrier H-bond formation is also observed at the terminus of the electron transfer pathway (Q_B and the H-bond partner D1-His215) when the second electron transfer from Q_A˙⁻ to Q_BH˙ occurs.¹⁷ D1-His215 is the ligand of the non-heme Fe complex, and the Nδ site forms an H-bond with the carbonyl O site of Q_B. The presence of the cationic Fe facilitates the deprotonation of D1-His215 and proton uptake by Q_BH⁻. Specifically, decreasing the pK_a value of ∼14¹⁸ for deprotonation of singly protonated histidine to the level of the pK_a value of ∼11¹⁹ for deprotonation of doubly protonated plastoquinone. It seems likely that the low-barrier H-bond between D1-His215 and Q_B is not stable because it should not inhibit the release of the product Q_BH₂ from the PSII binding site toward the quinone pool. A lingering question is why the redox-active low-barrier H-bond between TyrZ and D1-His190 is stable enough to function repeatedly in each S-state transition, while the bond between D1-His215 and Q_B functions only once before Q_BH₂ is released.

Initially, the redox-active group was identified as D1-His190 based on observations of radical formation in the Ca²⁺-depleted PSII.^20,21 Later, it was proposed that TyrZ, not D1-His190, was the origin of the radical state.²² Since then, TyrZ has been regarded as a redox-active group in the electron transfer pathway. Indeed, E_m(TyrZ) is lower than E_m(D1-His190), however the E_m difference is not significantly large when the proton is at the TyrZ moiety.¹⁰E_m values of the H-bond donor and acceptor moieties along low-barrier H-bonds have not been reported. Here, we report E_m of the two redox-active low-barrier H-bonds, considering the entire electron and proton transfer pathways quantum-chemically in the presence of the PSII protein environment.

Results

Low-barrier H-bond switching from an electron donor to an acceptor

The potential-energy profile of the H-bond shows that TyrZ and D1-His190 form a low-barrier H-bond. The difference between pK_a(TyrZ) and pK_a(D1-His190) is nearly zero along the low-barrier H-bond. However, the difference between E_m(TyrZ) (i.e., Tyr–O⁻/Tyr–O˙) and E_m(D1-His190), (i.e., N–His–N⁻/N–His–N˙ or HN–His–N/HN–His–N˙⁺) depends on the H⁺ position (Fig. 2, right panel). E_m(TyrZ) decreases and E_m(D1-His190) increases as the proton moves from TyrZ to D1-His190. When H⁺ arrives at the lower E_m moiety (D1-His190), E_m(TyrZ) is the lowest and the driving force for electron transfer from TyrZ to oxidized chlorophyll P_D1˙⁺ is the largest. The highest occupied molecular orbital (HOMO) is predominantly localized at the TyrZ moiety. Thus, the TyrZ moiety can serve as an electron donor to P_D1˙⁺ most effectively when the proton is at the D1-His190 moiety. Note that the H⁺ atom position does not affect the E_m values of other cofactors in the electron transfer pathway.

Fig. 2 Potential energy profile of the H-bond between TyrZ and D1-His190 (right vertical axis in kcal mol⁻¹) and the cofactor E_m values in the S₁ to S₂ transition (left vertical axis in mV). The Mn oxidation states are (Mn1, Mn2, Mn3, Mn4) = (III, IV, IV, III) in S₁ and (III, IV, IV, IV) in S₂, i.e., dangling Mn4 (Fig. 1b) is the oxidation site in the S₁ to S₂ transition. Left panels show the HOMO when (a) H⁺ is at the D1-His190 moiety (O_TyrZ⋯H = 1.38 Å) with the HOMO population of TyrZ : D1-His190 = 96.1 : 0.4 and (b) at the TyrZ moiety (1.02 Å) with the HOMO population of TyrZ : D1-His190 = 51.4 : 44.7. Red arrows indicate electron transfer from TyrZ to P_D1. Blue arrows indicate electron transfer (ET) from Mn4(III) to the TyrZ/D1-His190 pair. Cyan arrows indicate movement of the proton between the D1-His190 and TyrZ moiety. The gray dotted vertical line indicates the H⁺ position where E_m(TyrZ) and E_m(His190) are nearly equal and the TyrZ/His190 pair serves as an electron acceptor, finalizing the electron transfer process.

In contrast, when H⁺ is at the TyrZ moiety, E_m(TyrZ) and E_m(D1-His190) are equal, and the highest occupied molecular orbital (HOMO) is delocalized over the two moieties (Fig. 2b, left panel). In this case, the driving force for election transfer from the oxygen-evolving complex (i.e., Mn4^III/IV) to the TyrZ⋯D1-His190 pair is largest (Fig. 2b, right panel). Thus, the entire TyrZ⋯D1-His190 moiety cooperatively serves as an electron acceptor when electron transfer occurs from the oxygen-evolving complex.

Low-barrier H-bond maximizing electronic coupling and driving force

As an electron transfers to Q_BH˙, a low-barrier H-bond forms between D1-His215 and Q_BH˙/Q_BH⁻,¹⁷ facilitating Q_BH₂ release from the binding site in PSII. We analyzed E_m(N–His–N˙/N–His–N⁻) and E_m(Q_BH˙/Q_BH⁻) along the low-barrier H-bond between D1-His215 and Q_B. Note that E_m(Q_BH˙/Q_BH⁻) can be analyzed based on the (fully occupied) HOMO of Q_BH⁻, which is the same as the (singly occupied) HOMO (SOMO) of Q_BH˙.

(i) Driving force. E _m(Q_BH˙/Q_BH⁻) increases and E_m(D1-His215) decreases as H⁺ moves from D1-His215 to Q_BH˙/Q_BH⁻ along the low-barrier H-bond (Fig. 3b, right panel). E_m(Q_BH˙/Q_BH⁻) increases and the difference between E_m(Q_A^0/˙⁻) and E_m(Q_BH˙/Q_BH⁻) (the driving force for electron transfer) also increases when H⁺ moves to the lower E_m moiety (Q_BH˙/Q_BH⁻).

Fig. 3 Potential energy profile of the H-bond between D1-His215 and Q_B (right vertical axis in kcal mol⁻¹) and the Q_B MO energy [≈E_m(Q_BH˙/Q_BH⁻)] in the S₁ to S₂ transition (left vertical axis in eV). Left panels show the HOMO when (a) H⁺ is at the D1-His215 moiety (N_His215⋯H = 1.15 Å) with the HOMO population of Q_B : D1-His215 = 97.4 : 1.9 and (b) at the Q_B moiety (1.54 Å) with the HOMO population of Q_B : D1-His215 = 37.6 : 41.8. Blue arrows indicate electron transfer (ET) from Q_A˙⁻ to Q_BH˙. Cyan arrows indicate movement of the proton from the D1-His215 to the Q_B moiety. The gray dotted vertical line indicates the H⁺ position where E_m(His215) and E_m(Q_B) are nearly equal and the His215/Q_B pair serves as a terminal electron acceptor, finalizing the electron transfer process.

(ii) Electronic coupling. The HOMO is predominantly localized at the Q_BH˙/Q_BH⁻ moiety when H⁺ is at the D1-His215 moiety (Fig. 3a, left panel). E_m(D1-His215) decreases to the same level as E_m(Q_BH˙/Q_BH⁻) and the HOMO is delocalized over the D1-His215 and Q_B pair when H⁺ arrives at Q_BH˙/Q_BH⁻ (Fig. 3b, left panel). Both electron transfer from Q_A˙⁻ to Q_BH˙ and proton transfer from D1-His215 to Q_BH⁻ occur concertedly via the D1-His215 HOMO (= the Q_B HOMO), i.e., proton-coupled electron transfer. Note that the SOMO is not delocalized over the Fe complex, including other ligand groups (D2-His214 and bicarbonate) irrespective of the entire Fe complex being quantum-chemically considered, which is consistent with the absence of the redox change of the Fe complex in electron transfer from Q_A to Q_B. This suggests that the entire [D1-His215⋯H⁺⋯Q_BH˙] moiety serves as an electron acceptor for Q_A˙⁻ and decreases the substantial donor–acceptor distance from Q_A˙⁻⋯Q_BH˙ (12.7 Å) to Q_A˙⁻⋯D1-His215 (8.3 Å), increasing electronic coupling (Fig. 3b, left panel).

Discussion

The redox-active low-barrier H-bond between TyrZ and D1-His190 plays a dual role in electron transfer by serving as both an electron donor and an electron acceptor. The movement of the proton by only 0.35 Å from the D1-His190 moiety along the low-barrier H-bond alters E_m by ∼120 mV and switches its role from an electron donor to acceptor (Fig. 2). As the proton moves back to the TyrZ moiety in the low-barrier H-bond, E_m(TyrZ) and E_m(D1-His190) become equal within the barrier-less potential (O_TyrZ⋯H = 1.02 Å; Fig. 2b, right panel). Thus, TyrZ⋯D1-His190 can serve as both an electron donor and acceptor in the middle of the electron transfer pathway during the S-cycle without being destabilized by proton movement. Until now, it has been thought that D1-His190 is not redox-active, but that TyrZ is²² (however see ref. 20 and 21). However, this result suggests that both TyrZ and D1-His190 are redox-active components in the electron transfer pathway, in particular when abstracting an electron from the oxygen-evolving complex (Fig. 2b).

On the other hand, movement of the proton by only 0.35 Å along the low-barrier H-bond between D1-His215 and Q_BH˙/Q_BH⁻ increases E_m(Q_BH˙/Q_BH⁻) to the level of E_m(D1-His215), leading to delocalization of the HOMO over the D1-His215 and Q_B moiety. The involvement of D1-His215 in the electron acceptor decreases the donor-to-acceptor distance from [Q_A]-to-[Q_B] to [Q_A]-to-[D1-His215] and increases the electronic coupling between the donor Q_A and the acceptor Q_BH˙/Q_BH⁻, because the electronic coupling increases as the electron donor-to-acceptor distance decreases.^23,24 The role of the redox-active low-barrier H-bond in increasing electronic coupling is observed specifically for D1-His215⋯Q_BH˙/Q_BH⁻, since the axis of the low-barrier H-bond is consistent with the axis of the electron transfer pathway and delocalization of the HOMO over D1-His215 and Q_B decreases the substantial donor–acceptor distance (Fig. 3, left panel).

In contrast to the TyrZ⋯D1-His190 H-bond, the difference between E_m(D1-His215) and E_m(Q_B) is small but not zero within the barrier-less potential as the proton moves to the Q_B moiety (H_D1-His215⋯H = 1.43 Å; Fig. 3b, right panel). Indeed, the HOMO is still predominantly localized at Q_B (Fig. S1, ESI†). Only after the proton exceeds the barrier-less potential by 0.1 Å is the HOMO evenly delocalized over D1-His215 and Q_B and the difference between E_m(D1-His215) and E_m(Q_B) reaches zero (H_D1-His215⋯H = 1.53 Å, Fig. 3b, right panel). In reality, the low-barrier H-bond is unlikely to exist when H_D1-His215⋯H = 1.53 Å because the potential energy profile was obtained assuming the presence of the H-bond (see Methods). Consistently, the donor–acceptor N⋯O distance (N_D1-His215⋯O_QB) is 2.47 Å when the proton is located at the D1-His215 moiety (N_D1-His215–H⋯O_QB), and it is lengthened to 2.55 Å when the proton arrives at the Q_B moiety (N_D1-His215⋯H–O_QB). Thus, Q_BH₂ is released from the binding site (D1-His215) when H⁺ leaves the barrier-less potential (after cleavage of the low-barrier H-bond). The irreversibility of the reaction is characteristic of, and required for, Q_B serving as the terminal electron acceptor in the electron transfer pathway. Intriguingly, the corresponding change in the donor–acceptor distance in response to the H⁺ movement is absent in the TyrZ⋯D1-His190 H-bond (O_TyrZ–H⋯N_D1-His190 = 2.52 Å and O_TyrZ⋯H–N_D1-His190 = 2.50 Å). This is consistent with the role of the TyrZ⋯D1-His190 pair in the middle of the electron transfer pathway, which involves reversibly donating and accepting an electron during the entire S-cycle. As reported, environmental fluctuations may affect the potential energy profile,²⁵ particularly for standard H-bonds. Proton transfer is energetically uphill in the standard H-bond (Fig. 1c). To overcome this energy barrier, fluctuation/rearrangement of the protein environment is a prerequisite for unstable protonation of the acceptor moiety. In contrast, proton transfer is barrier-less in the low-barrier H-bond, which does not require the corresponding fluctuation/rearrangement of the protein environment (Fig. 1d). Consistently, the B-factors of TyrZ and D1-His190 are specifically low¹¹ because the H-bond network is fixed by the PSII protein electrostatic environment, namely the Mn₄CaO₅ cluster and the ionized ligand residues.¹²

For the redox-active low-barrier H-bonds to serve as a redox-active cofactor in the electron transfer pathway, the presence of the protein electrostatic environment is a prerequisite. First, TyrZ and D1-His190 cannot form a low-barrier H-bond in the absence of the PSII protein environment¹² because of the difference in pK_a (∼3) between tyrosine and histidine. Second, the protein electrostatic environment is required for the formation of the downhill electron transfer pathway that proceeds from the Mn₄CaO₅ cluster via TyrZ, as it increases E_m(P_D1) by ∼400 mV with respect to E_m(Chla).^6,26 Thus, the environment, in which the pK_a values of the proton donor and acceptor moieties are equal and the E_m values of the electron donor and acceptor moieties are in the tunable range by proton transfer, is ultimately provided by the common protein electrostatic environment.

In summary, the pK_a values of the two moieties are equal, but the E_m values depend on the H⁺ position in redox-active low-barrier H-bonds. The results also suggest that redox-active low-barrier H-bonds can differ in their characteristics, whether they serve as a redox-active cofactor in the middle of the electron transfer pathway (TyrZ⋯D1-His190) or in the terminus of the electron acceptor pathway (D1-His215⋯Q_B). These findings provide a key to understanding how nature optimizes electron and proton transfer in biological systems using abundantly available protons.

Methods

The atomic coordinates of PSII were obtained from the PSII crystal structure (PDB code, 3ARC).¹¹ The atomic charges of the other cofactors in the MM region were taken from a previous study.²⁷

E _m calculations for the electron transfer pathway via TyrZ⋯D1-His190

The HOMO energy level is largely correlated with E_m for one-electron oxidation (e.g., ref. 28 and 29). For the Mn₄CaO₅ cluster, the HOMO in S_n corresponds to the molecular orbital, which predominantly contributes to the release of an electron in the S_n to S_n+1 transition. We included all redox-active cofactors (Mn₄CaO₅ cluster, TyrZ, P_D1, and P_D2) simultaneously in the QM region^30,31 (see below), identified HOMOs of Mn₄CaO₅, TyrZ, P_D1, and P_D2 in S₁ on the basis of the Mulliken population analysis³² (results provided in a previous study¹⁰), and obtained the E_m values (note, E_m(S₁/S₂) for the Mn₄CaO₅ cluster). The E_m values of the redox sites are calculated using eqn (1):

E_m = −0.15499E_HOMO + 1205.65, (1)

where E_HOMO is the HOMO energy level (meV).¹⁰

QM/MM calculations

The unrestricted density functional theory method was employed with the B3LYP functional using the QSite program.³³ B3LYP is a widely used functional to calculate redox potentials and HOMO and LUMO energy levels (e.g., ref. 34 and 35). Using B3LYP can also make a comparison with our previous studies for Mn₄CaO₅,^10,36 TyrZ,^12,13 P_D1,^6,37 and Q_B¹⁷ in PSII. In the QM region, all the atomic coordinates were fully relaxed. In the MM region, the positions of the H atoms were optimized using the OPLS2005 force field,³⁸ while the positions of the heavy atoms were fixed.

To investigate the H-bond between TyrZ and D1-His190, the LACVP* basis set was employed. The Mn₄CaO₅ cluster was considered to be in S₁ with antiferromagnetically coupled Mn ions; the resulting Mn oxidation states (Mn1, Mn2, Mn3, Mn4) and the total spin S were S = 8/2 (↑ ↓ ↑↑) in S₁. Note that the difference in S (for example, S = 0 in S₁³⁹) did not affect the calculated E_m values (≤3 mV).¹⁰ O1–O5 was considered to be unprotonated (O²⁻). The four water ligand molecules, W1–W4, were considered to be water molecules (H₂O). The Mn₄CaO₅ geometry was obtained from our previous studies.^40,41 The initial-guess wavefunctions were obtained using the ligand field theory⁴² implemented in the QSite program. The QM region was defined as the Mn₄CaO₅ cluster (including the ligand side-chains of D1-Asp170, D1-Glu189, D1-His332, D1-Glu333, D1-Asp342, and CP43-Glu354, and backbone of D1-Ala344), ligand water molecules of W1–W4, O4–water chain (W539, W538, and W393), Cl-1 binding site (Cl-1, W442, W446, and side-chains of D1-Asn181 and D2-Lys317), second-sphere ligands (side-chains of D1-Asp61 and CP43-Arg357), the H-bond network of TyrZ (side-chains of D1-Tyr161, D1-His190, and D1-Asn298), diamond-shaped cluster of water molecules¹² (W5, W6, and W7), and P_D1/P_D2. The MM region was defined as the entire PSII protein, as in a previous study.¹⁰ D1-His337 was considered to be protonated,⁴³ whereas all other titratable groups were in the standard protonation states. See ref. 10 for the QM/MM-optimized atomic coordinates.

To investigate the H-bond between D1-His215 and Q_B, the LACVP**+ basis set was employed. The QM region was defined as [Q_B, D1-His252, D1-Ser264, bicarbonate, Fe, D1-His215, D1-His272, D2-His214, and D2-His268]. We assumed a high-spin state (S = 2) of Fe²⁺ and set the spin multiplicity of the system to S = 2 in the calculations for Q_BH⁻ and Q_BH₂. The MM region was defined as the D1 and D2 protein subunits, as in a previous study.¹⁷ D1-His252 was considered to be protonated before Q_BH₂ formation,⁴⁴ whereas all other titratable groups were in the standard protonation states. See ref. 17 for the QM/MM-optimized atomic coordinates.

To obtain the potential energy profiles of the O⋯H⋯N H-bond, the QM/MM optimized geometry was used as the initial geometry. The H atom under investigation was moved between O and N by 0.05 Å before the geometry was optimized by constraining the O–H and H–N distances, and the energy was calculated. This procedure was repeated until the H atom reached the O and N atoms. This approach, which is based on the single QM/MM-optimized geometry, provides the unique minimum-energy pathway, in particular for the H-bond between TyrZ and D1-His190, as demonstrated by analyzing proton-transfer pathways in the different protein conformations.¹²

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research was supported by JST CREST (JPMJCR1656 to H. I.), JSPS KAKENHI (JP18H05155, JP18H01937, JP20H03217, and JP20H05090 to H. I., JP18H01186 to K. S., and JP16H06560 to K. S.), and the Interdisciplinary Computational Science Program in CCS, University of Tsukuba (K. S.).

References

1. J. R. Shen, Annu. Rev. Plant Biol., 2015, 66, 23–48.
2. H. Tamura, K. Saito and H. Ishikita, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 16373–16382.
3. S. E. J. Rigby, J. H. A. Nugent and P. J. O'Malley, Biochemistry, 1994, 33, 10043–10050.
4. B. A. Diner, E. Schlodder, P. J. Nixon, W. J. Coleman, F. Rappaport, J. Lavergne, W. F. J. Vermaas and D. A. Chisholm, Biochemistry, 2001, 40, 9265–9281.
5. T. Okubo, T. Tomo, M. Sugiura and T. Noguchi, Biochemistry, 2007, 46, 4390–4397.
6. K. Saito, T. Ishida, M. Sugiura, K. Kawakami, Y. Umena, N. Kamiya, J.-R. Shen and H. Ishikita, J. Am. Chem. Soc., 2011, 133, 14379–14388.
7. V. V. Klimov, S. I. Allakhverdiev, S. Demeter and A. A. Krasnovskii, Dokl. Akad. Nauk SSSR, 1979, 249, 227–230.
8. A. W. Rutherford, J. E. Mullet and A. R. Crofts, FEBS Lett., 1981, 123, 235–237.
9. F. Rappaport, M. Guergova-Kuras, P. J. Nixon, B. A. Diner and J. Lavergne, Biochemistry, 2002, 41, 8518–8527.
10. M. Mandal, K. Kawashima, K. Saito and H. Ishikita, J. Phys. Chem. Lett., 2020, 11, 249–255.
11. Y. Umena, K. Kawakami, J.-R. Shen and N. Kamiya, Nature, 2011, 473, 55–60.
12. K. Saito, J.-R. Shen, T. Ishida and H. Ishikita, Biochemistry, 2011, 50, 9836–9844.
13. K. Kawashima, K. Saito and H. Ishikita, Biochemistry, 2018, 57, 4997–5004.
14. A. Warshel, A. Papazyan and P. A. Kollman, Science, 1995, 269, 102–106.
15. C. N. Schutz and A. Warshel, Proteins, 2004, 55, 711–723.
16. C. L. Perrin and J. B. Nielson, Annu. Rev. Phys. Chem., 1997, 48, 511–544.
17. K. Saito, A. W. Rutherford and H. Ishikita, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 954–959.
18. T. C. Bruice and G. L. Schmir, J. Am. Chem. Soc., 1958, 80, 148–156.
19. R. Hasegawa, K. Saito, T. Takaoka and H. Ishikita, Photosynth. Res., 2017, 133, 297–304.
20. A. Boussac, J.-L. Zimmermann, A. W. Rutherford and J. Lavergne, Nature, 1990, 347, 303–306.
21. A. Boussac and A. W. Rutherford, Biochemistry, 1992, 31, 7441–7445.
22. B. J. Hallahan, J. H. A. Nugent, J. T. Warden and M. C. W. Evans, Biochemistry, 1992, 31, 4562–4573.
23. C. C. Moser, J. M. Keske, F. Warncke, R. S. Farid and P. L. Dutton, Nature, 1992, 355, 796–802.
24. C. C. Page, C. C. Moser, X. Chen and P. L. Dutton, Nature, 1999, 402, 47–52.
25. T. Vasilevskaya, M. G. Khrenova, A. V. Nemukhin and W. Thiel, J. Comput. Chem., 2016, 37, 1801–1809.
26. H. Ishikita, W. Saenger, J. Biesiadka, B. Loll and E.-W. Knapp, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 9855–9860.
27. K. Kawashima and H. Ishikita, Chem. Sci., 2018, 9, 4083–4092.
28. T. Watanabe and M. Kobayashi, in Chlorophylls, ed. H. Scheer, CRC Press, Boca Raton, FL, 1991, pp. 287–303.
29. D. D. Méndez-Hernández, P. Tarakeshwar, D. Gust, T. A. Moore, A. L. Moore and V. Mujica, J. Mol. Model., 2013, 19, 2845–2848.
30. S. G. Abuabara, L. G. C. Rego and V. S. Batista, J. Am. Chem. Soc., 2005, 127, 18234–18242.
31. P. Joshi, V. Shewale, R. Pandey, V. Shanker, S. Hussain and S. P. Karna, J. Phys. Chem. C, 2011, 115, 22818–22826.
32. R. S. Mulliken, J. Chem. Phys., 1955, 23, 1833–1840.
33. QSite, version 5.8, Schrödinger, LLC, New York, NY, 2012.
34. G. Zhang and C. B. Musgrave, J. Phys. Chem. A, 2007, 111, 1554–1561.
35. D. Coskun, S. V. Jerome and R. A. Friesner, J. Chem. Theory Comput., 2016, 12, 1121–1128.
36. K. Saito, M. Mandal and H. Ishikita, Biochemistry, 2020, 59, 3216–3224.
37. H. Tamura, K. Saito and H. Ishikita, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 16373–16382.
38. J. L. Banks, H. S. Beard, Y. X. Cao, A. E. Cho, W. Damm, R. Farid, A. K. Felts, T. A. Halgren, D. T. Mainz, J. R. Maple, R. Murphy, D. M. Philipp, M. P. Repasky, L. Y. Zhang, B. J. Berne, R. A. Friesner, E. Gallicchio and R. M. Levy, J. Comput. Chem., 2005, 26, 1752–1780.
39. D. Koulougliotis, D. J. Hirsh and G. W. Brudvig, J. Am. Chem. Soc., 1992, 114, 8322–8323.
40. K. Saito, A. W. Rutherford and H. Ishikita, Nat. Commun., 2015, 6, 8488.
41. K. Kawashima, T. Takaoka, H. Kimura, K. Saito and H. Ishikita, Nat. Commun., 2018, 9, 1247.
42. G. Vacek, J. K. Perry and J. M. Langlois, Chem. Phys. Lett., 1999, 310, 189–194.
43. S. Nakamura and T. Noguchi, J. Am. Chem. Soc., 2017, 139, 9364–9375.
44. H. Ishikita and E.-W. Knapp, J. Am. Chem. Soc., 2005, 127, 14714–14720.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0cp04265j

‡ These authors contributed equally to this work.

---