The Electronic Origin of Far‐Red‐Light‐Driven Oxygenic Photosynthesis

Abstract Photosystem‐II uses sunlight to trigger charge separation and catalyze water oxidation. Intrinsic properties of chlorophyll a pigments define a natural “red limit” of photosynthesis at ≈680 nm. Nevertheless, charge separation can be triggered with far‐red photons up to 800 nm, without altering the nature of light‐harvesting pigments. Here we identify the electronic origin of this remarkable phenomenon using quantum chemical and multiscale simulations on a native Photosystem‐II model. We find that the reaction center is preorganized for charge separation in the far‐red region by specific chlorophyll–pheophytin pairs, potentially bypassing the light‐harvesting apparatus. Charge transfer can occur along two distinct pathways with one and the same pheophytin acceptor (PheoD1). The identity of the donor chlorophyll (ChlD1 or PD1) is wavelength‐dependent and conformational dynamics broaden the sampling of the far‐red region by the two charge‐transfer states. The two pathways rationalize spectroscopic observations and underpin designed extensions of the photosynthetically active radiation limit.

For the RESP fitting purposes, only Fe(II) and HCO3were considered in the cluster because we observed underestimated charges for the Fe(II) site when the complete cluster is considered. Therefore, the remaining amino acids were given the corresponding charges from the AMBER force-field library. The Fe(II) and HCO3based cluster models were first optimized at B3LYP/Def2-TZVP and then RESP fitting is performed at the B3LYP/6-31G* level. In both cases, we restrained the charge of the backbone atoms of the residues according to the original AMBER force field [6b] as such procedure is known to produce better back-bone dynamics during the simulation. [9] The bonded parameters of all organic cofactors were based on the General AMBER Force Field (GAFF2). The bonded parameters for the non-heme iron site, [10] chlorophyll a [11] and heme [12] site were taken from the literature. On the other hand, custom bonded parameters [13] were used for the oxygen evolving cluster, Mn4CaO5. The parameters for the protein residues, water and lipid bilayer were defined based on the Amber14SB, [14] TIP3P [15] and LIPID17 force-fields, [16] respectively.

Molecular dynamics simulations
A systematic minimization protocol is employed to remove energetically unfavorable geometric clashes. In the equilibration phase, we systematically heated the system from 10 to 100 K in a succession of 5 ps in the NVT ensemble. After this, the system was further heated from 100 to 303 K in the NPT ensemble for a total of 125 ps. The temperature during this step is maintained using the Langevin dynamics [17] with a collision frequency of 5 ps -1 . During the equilibration phase, the Cα atoms of amino acids were restrained with a force constant of 20 kcal mol -1 Å -2 . Subsequently, the restrains on the Cα atoms of amino acids were systematically released (2 kcal mol -1 Å -2 /400 ps). Thereafter, the Monte Carlo/Molecular Dynamics (MC/MD) module [18] of the AMBER package is employed for a controlled hydration and de-hydration from bulk to protein, and vice versa. Finally, the system is further equilibrated for 63 ns (without S6 any restraints) to properly equilibrate the bilayer. Thereafter, a production run was initiated for a total of 200 ns, the temperature and pressure were set at 303 K and 1 atm, respectively ( Figure   S3). The temperature was controlled using Langevin dynamics [17] with a collision frequency of 1 ps -1 . During the entire procedure, the pressure of the system is controlled using the Berendsen barostat [19] with anisotropic pressure scaling with a relaxation time of 2 ps. We employed the SHAKE algorithm [20] to constrain the hydrogen atom involving bonds, therefore, a time step of 2 fs could be used. The frames were saved every 2 ps. The electrostatic interactions were treated using the Particle Mesh Ewald (PME) [21] approach with a 10 Å cut-off. The AMBER18 package [22] was used to perform the energy minimizations and molecular dynamics simulations.
The production simulations were performed GPU version of the pmemd module (pmemd.cuda) of the AMBER18 package.

QM/MM calculations
For QM/MM computations, we extracted a total of 22 snapshots of the complete system from the molecular dynamics simulations. The first snapshot (i.e. snapshot 1), which corresponds to the structural configuration close to the crystal structure, is derived after clustering [23] a series of frames from the early equilibration phase. The other 21 snapshots were extracted from the unbiasedly evolved production simulations of the PSII-bilayer system. The protein structural configurations in this case are more equilibrated and completely hydrated. To ensure that these protein configurations were structurally uncorrelated, we chose 21 snapshots (starting from the first frame of the production simulations) that are equally spaced from each other with a consecutive time interval of 10 ns. For the QM/MM computation, we considered the entire PSII monomer along with all the waters molecules (including all waters within protein cavities) which are 7 Å around the protein (includes around 8000 water molecules). The overall system was kept neutral using the appropriate amount of Na + ions placed at the equilibrated positions.
All QM/MM optimizations were performed using the ChemShell 3.7 package. [24] The MM part of the computation was handled by the in-built DL-POLY module, [25] whereas the ORCA 4.2 package [26] was employed as the QM engine. All QM/MM computations are based on the electrostatic embedding technique. The covalent bonds across the QM and MM region were cut using the hydrogen link atom approach. Additionally, we used the charge-shift method [27] to avoid over-polarization of the QM region by closely placed MM charges. Due to the large spatial extent and number of atoms in such assemblies, a direct "brute-force" QM/MM optimization of the whole RC is extremely expensive. Instead, we performed optimizations of the individual pigments or pairs of pigments, and subsequently the optimized configurations were grouped to form larger QM oligomeric assemblies. For any given snapshot, we carried out individual QM/MM geometry optimization of the ChlD1, ChlD2, PheoD1 and PheoD2. The ChlD1-PheoD1 pair on the active side was optimized as a single unit, and the PD1-PD2 pair was optimized as a single unit in order to properly reproduce the symmetry breaking feature of the vinyl group. In addition, with such supramolecular approach, we ensure that we retain the macrocyclic ring curvature of the PD1 and PD2, which is critical in determining the energetics of the Q-band. [28] We always considered the axially ligated residues in the QM region during geometry optimization, for example PD1 and PD2 are axially ligated by the His198 and His197, respectively, whereas ChlD1 and ChlD2 are axially ligated by a single water molecule.
In the case of the ChlD1 and ChlD2, the water molecule hydrogen bonded to the axially ligated S9 water and ester group attached to the 13 2 carbon position on ring E is also included in the QM region. In all cases, the phytyl chain of the pigments is treated at the MM level and the link hydrogen atom is placed between the C17 1 and C17 2 atoms.
The geometry optimization protocol adopted in this work includes division of the entire system into two parts, active and static. The atoms in the active region are free to move during the geometry optimization step, whereas the rest of the atoms form part of the static region and remain fixed during geometry optimization, acting as an electrostatic environment to the active region. In geometry optimization of the ChlD1, ChlD2, PheoD1 and PheoD2, protein components within 13 Å (ca. 1300-1600 atoms) around the pigments were considered in the active region.
In case of ChlD1-PheoD1 pair, an active region of 15 Å (ca. 2560 atoms) was chosen. In case of the optimization of the "special-pair" PD1-PD2, a large active region (~15 Å, ca. 2650 atoms) was selected. All geometry optimizations were performed with the PBE functional [29] and Def2-TZVP basis set, [30] along with the D3(BJ) dispersion corrections. [31] We employed a higher DFT integration grid (Grid6 in ORCA 4.2 convention) during our calculations. The computations of the Coulomb integrals used the resolution of identity (RI) approximation [32] with matching auxiliary basis sets. [33] The QM/MM geometry optimization convergence criterion for the maximum gradient component was set to 4.5×10 -4 Hartree/Bohr. The convergence criterion for energy change was set to 1.0×10 -6 Hartree. The criterion for the maximum step components was 1.8×10 -3 Bohr, in case of the RMS gradient 3.0×10 -4 Hartree/Bohr and for the RMS step the tolerance was 1.2×10 -3 Bohr.
In the case of the PD1-PD2-ChlD1-PheoD1 tetramer, QM/MM optimized geometries from ChlD1-PheoD1 and PD1-PD2 pair optimizations were combined, whereas in the case of PD1-PD2-ChlD1-ChlD2, QM/MM optimized geometries from individual ChlD1 and ChlD2, and PD1-PD2 pair optimizations were combined. To confirm that this approach yields accurate oligomeric geometries and spectroscopic properties, we performed for one MD snapshot (snapshot 1) direct S10 geometry optimizations of the entire PD1-PD2-ChlD1-PheoD1 and of the PD1-PD2-ChlD1-ChlD2 assemblies. We used PBE/def2-TZVP-MM level of theory, yielding a total of 6366 and 6472 basis-functions for PD1-PD2-ChlD1-PheoD1 and PD1-PD2-ChlD1-ChlD2, respectively. A large active region of 24 Å was chosen in both cases. Structural overlay of the geometries (Figure S5) obtained from our approach and brute optimization yielded a root mean square deviation (RMSD) of 0.090 Å and 0.100 Å (hydrogen atoms not included in the calculation) for PD1-PD2-ChlD1-PheoD1 and PD1-PD2-ChlD1-ChlD2, respectively. Therefore, the geometries obtained by directly optimizing a tetramer compared to combining pairwise optimized structures are essentially identical, as expected from the fact that the specific pigment pairs are structurally uncoupled. We note that tetramer optimizations are exceedingly costly in terms of time primarily due to dramatically slower convergence of the gradient, a result precisely of the shallow potential energy surface resulting from the lack of well-defined internal coordinates between the central pair and the accessory pigments. As expected from the above, excited state calculations using tetramer-optimized geometries provide the same picture of low-lying excited states in terms of the nature and energetics of the lowest CT and excitonic transitions. Figure S5. Structure overlay of tetrameric geometries obtained from brute force optimization (blue) and geometry assembly approach (red) adopted in the current work. Hydrogen atoms are not shown for clarity.

Excited state calculations
Vertical excitation energies (VEE) were computed for all snapshots of the QM/MM optimized individual pigments and larger assemblies. The explicit electrostatic effect of the protein environment on the computed properties is included through point charges. The effect of the point charges is handled by the orca_pc module, which adds the contribution of the charges in the one-electron matrices and nuclear repulsion. All excited state calculations were performed using the ORCA 5.0 package. We performed full Time-Dependent Density Functional Theory (TD-DFT), i.e. without the Tamm-Dancoff approximation, to compute the excited state properties using the range-separated ωB97X-D3(BJ) [34] functional along with the Def2-TZVP basis set. The choice of this functional is motivated by the earlier successful direct comparison with the results obtained using the similarity transformed equation of motion coupled cluster theory. [35] This functional performs best in its class in reproducing the correct excitation asymmetry and electrochromic shifts within the reaction center. [35c] This functional was also successfully applied to study the excited state properties of pigments in the CP47 core antenna complex of Photosystem II. [36] In addition, various other studies have demonstrated the good performance of this functional for charge-transfer states. [37] We used the RIJCOSX approximation [38] is used to speed up the calculations. In all the cases considered in this work, very tight SCF convergence criteria were set along with the new DefGrid2 DFT integration grids of ORCA 5.0. The study of excited state properties of the individual pigments is crucial to study the site-energy (Qy excitation), therefore we computed the first 10 roots. However, in the case of the tetrameric assemblies PD1-PD2-ChlD1-PheoD1 and PD1-PD2-ChlD1-ChlD2 we computed the first 12 roots due to the high cost of these calculations (total of ~6500 basis functions). We observed that the first 12 roots of the tetrameric assembly essentially cover the Q-band energy range and provide all necessary information about the critical excited states with local excitation (LE), charge-transfer (CT), and mixed character. S12

Excited state properties of individual pigments and assemblies
In this section, the TD-DFT results obtained from the computation on individual pigments and the tetrameric assemblies, i.e. PD1-PD2-ChlD1-PheoD1 and PD1-PD2-ChlD1-ChlD2, are tabulated. In case of the individual pigments, we report the lowest excited state of the pigments (also known as 'site-energy'). In case of the tetrameric assemblies, all excited states along with their respective character are tabulated for snapshot 1. In addition, we also performed excited state calculations (on both tetramers) without the point charge field of the protein in order to quantify the electrostatic effect of the protein matrix (performed on snapshot 1).
The first four excited states corresponding to the low-energy spectrum (red and far-red region) of PD1-PD2-ChlD1-PheoD1 are listed for all the snapshots under investigation, since this assembly representing the active branch of PSII forms the central focus of the current work. Individual states were assigned a definitive character based on the analysis of the canonical and natural transition orbitals (NTOs). The oscillator strength (f) corresponding to the excited states are also shown. It becomes evident from these results that the conformational motion of the electrostatic matrix of the protein greatly influences primarily the energetics of the CT states in contrast to the local excitations. Table S1. The lowest energy excited state (in eV) of all individual pigments of the reaction center for all snapshots considered in this work. In each case, the geometry is optimized in the QM/MM framework at the PBE/Def2-TZVP-MM level of theory. The excited state calculations are performed within the TD-DFT framework using the ωB97X-D3(BJ)/Def2-TZVP-MM level of theory. The effect of the protein environment is included using point charges. Corresponding oscillator strengths are indicated in parentheses. The redmost excitation energy within a snapshot is marked in bold. *Although the site-energy (Qy) of ChlD1 pigment is the lowest among all pigments for a majority of snapshots, however, it is noted that our past [35a] and current results suggest that Qy band of ChlD1 is generally mixed with a CT state of ChlD1 ＋ PheoD1character, which leads to a red-shift (see excited state results from the tetramer assemblies in the following). Therefore, drawing direct conclusions from site energies of the individual pigments alone is partly unjustified.     Table S6. Low-energy spectrum of PD1-PD2-ChlD1-PheoD1 tetramer (snapshot 2) with the protein environment. This snapshot corresponds to the first frame of the production simulations.

State VEE
f Nature and NTO Contribution S1      Table S3.  Table S3. It can be observed that this CT state has admixture of local excitation character resulting in slightly higher oscillator strength compared to the PD1 + PheoD1 − CT State.  Table S3. This excited state has a majority of a local excitation character mixed with small CT character.  Table S5. Detailed list of all computed excited states can be found in Table S5.

Parallel charge separation pathways
Various ideas about the existence of more than one charge separation pathway in the RC of PSII have been discussed in the past, [39] typically based on investigations within the "red-limit" (~680 nm). Conformational dynamics [39c] or the wavelength of the incident photon [39e] were usually considered as independent factors, whereas here we give weight on both. In the following we mention some selected previous proposals and clarify differences and similarities to the model resulting from the present work.
One such model has been proposed by Romero [42] and results of the current study, can be used as basis for future experiments to understand wavelength dependence of charge separation pathways, especially for the PD1 case.

YD (D2-Tyr160) as probe for the charge separation pathway
Direct identification of the charge separation pathways triggered by far-red light is difficult due to the dark nature of the corresponding charge-transfer states. One indirect way is to exploit the spectroscopic signals of the secondary electron donors (e.g. YD, ChlZ, Cytb559 etc.) of PSII, as shown by Styring and coworkers. [42] The secondary electron donors play critical roles such as photo-protection and regulation. [43] Upon electron donation towards the RC, these species become EPR-active showcasing unique signals important for their identification. The work by Styring and coworkers [42] used spectroscopic signals from YZ, Cytb559, ChlZ and CarD2 S32 secondary electron donors. Similarly, Ahmadova and Mamedov [44] exploited the redox properties and EPR signals of the YZ and YD radicals using various wavelengths to probe charge separation. Based in our results, it can be proposed that EPR signals from the YD radical can be used in the identification of far-red light driven charge separation pathways. The mechanism of the oxidation of the YD is understood since high-level spectroscopy-oriented quantum chemical analysis explained why YD exhibits biphasic slow and fast oxidation under low and high pH conditions, respectively. [44][45] It has been observed that (a) under high pH (~8.5) conditions YD outcompetes YZ in terms of electron donation rates to the PD1-PD2 pair, [46] and (b) YD can be oxidized even at 1.7 K under high pH, which gives rise to a unique EPR signal with gx=2.0064. [45c, 47] Based on this, it will be interesting to conduct experiments with PSII in high pH and cryogenic conditions investigating the EPR signal of YD in combination with far-red illumination. It can be speculated that light in the 700-750 nm range would most likely activate the ChlD1 pathway, which would probably give rise to few PSII centers with a YD radical due to the large distance between ChlD1 and YD (see Figure S12), while if the PD1 pathway can be activated within the 750-800 nm range one might expect more PSII centers with a YD radical.

Far-red light acclimated photosynthetic organisms
Photosynthetic organisms can enhance their far-red light (FRL) absorbing capabilities (beyond 700 nm) by (a) replacing some pigments with the long-wavelength absorbing pigments, and (b) by tweaking protein electrostatics, enabling direct excitation of the low-lying CT states as shown in the current work. The former case is known only in certain photosynthetic species, whereas the latter may apply to nearly all photosynthetic species. For example, a few cyanobacterial species when subjected to far-red light synthesize the far-red light absorbing Chl d and Chl f. [48] Acaryochloris marina [49] is one of the cyanobacteria which uses Chl d as a major pigment. With its pigment composition, it can easily utilize the long-wavelength photons up to ca. 720 nm. On the other hand, the recently discovered cyanobacterium Chroococcidiopsis thermalis [50] produces Chl f [48a, 51] as a minor pigment, in addition to Chl a and Chl d, when Chroococcidiopsis thermalis is known to employ a mixture of Chl a, Chl d and Chl f in PSII and PSI to harvest far-red light. However, the location and identity of these pigments in the RC remains under investigation. Nürnberg et al. [50] proposed that the Chl f occupies the ChlD1 position in the RC for facilitating FRL photochemistry. Judd et al. [52] proposed ChlD1 or PD2 position in the RC of PSII to be occupied by a Chl f pigment. Zamzam et al. [53] proposed that the initial electron donor in the reaction center is a FRL absorbing pigment (Chl f) with a likelihood of occupying the ChlD1 position. The reasoning [50] behind placing the far-red pigments in the reaction centers is that it would be otherwise an uphill excitation energy (~100 S34 meV) transfer from the Chl f antenna pigments to the all Chl a based reaction center.
Interestingly, Croce and coworkers [54] showed that FRL adapted PSII samples show a much slower charge separation rate. Intuitively, if Chl f is indeed present in the RC of PSII, then overall charge separation rate must be at least comparable to the simple PSII (white light grown, all Chl a based photosystems) due to global optimization of the excitation energy transfer and charge separation.
Similarly, variable suggestions and observations exist in case of the FRL light driven photochemistry of PSI. Nürnberg et al. [50] proposed that one or both accessory pigments positions, i.e. AA and AB, in the RC are occupied by Chl f (see Figure S13 for the PSI-RC cofactor nomenclature). Similar proposals were made by Hastings et al. [55] based on Fourier transform visible and infrared difference spectroscopy performed on Fischerella thermalis PCC 7521 and by Zamzam et al. [56] using ultrafast transient infrared absorption spectroscopy. In contrast, Chl f was only assigned an antenna role in the study by Cherepanov et al., [57] who found that it takes nearly 230 ps to develop the charge separated state, proposing that this is due to uphill energy migration from the Chl f present in the antenna to the Chl a based RC. Similar observations regarding slow charge separation rates were made by Mascoli et al. [54] and Tros et al. [58] The work by Tros et al. [58] also weighs in favor of an all Chl a based PSI-RC, suggesting that Chl f in the antenna transfers excitation energy to low energy charge-transfer states of the RC (i.e. below the P700 energy band). The presence of such low-energy CT states in the PSI RC of simple PSI (white light grown, all Chl a based photosystems) has been proposed earlier by Mokvist et al. [59] and Schlodder et al. [60] The crystal structures [61] of the FRL acclimated PSI also confirms the absence of Chl f from the RC. However, it must be pointed out that cryoelectron microscopy (cryo-EM) [62] also suffers from intrinsic limitations in differentiating between Chl a and Chl f.

S35
Based on the brief discussion presented on the FRL driven PSI and PSII photochemistry, there exist many unanswered questions, such as: (1)

Photosystem II compared to the bacterial reaction center
The photosynthetic reaction center from purple bacteria (Rhodobacter sphaeroides and The growing consensus in the study of charge separation in PSII considers ChlD1 to be the primary electron donor. [ Table S27). Despite fundamental similarities between the PSII and BRC, the overall charge separation picture deviates significantly from a conceptual point of view. Therefore, rather than using the BRC as a model to understand PSII charge separation, it is better to consider the two systems independently.

Differences in the protein matrix of Photosystem II from different species
The results in the present work are based on the structure of cyanobacterial PSII derived from Thermosynechococcus vulcanus. The protein sequences of the D1 and D2 proteins from the higher plants, green algae and cyanobacteria have slight differences (Figures S15 and S16).
Some of these differences are near the RC pigments, for example D1-Gln130 in cyanobacteria is replaced by D1-Glu130 in green algae and higher plants. Cyanobacteria also carry different forms of the D1 protein (encoded by the psbA gene) leading to various changes in the protein matrix. [70] Based on earlier [35a, 35c, 36] and current work, protein matrix electrostatics are critical for creating and maintaining asymmetry within the RC, modulating the site energies and CT energetics. Therefore, from a fundamental of point of view it would be interesting to correlate the effect of key differences in the D1 and D2 sequence among organisms to excited state energetics and charge separation kinetics. This analysis is currently pursued in our lab in the hope that it will not only enhance our understanding of the response of organisms to variable environment effects, but also provide opportunities for enzyme engineering.

MK-RESP charges for the Oxygen Evolving Complex and its protein based components
Note: The internal connectivity within each component is also defined along with the assigned atom types.

MK-RESP charges for the non-heme iron site and its protein based components
Note: The internal connectivity within each component is also defined along with the assigned atom types.