Protein Matrix Control of Reaction Center Excitation in Photosystem II

Photosystem II (PSII) is a multisubunit pigment–protein complex that uses light-induced charge separation to power oxygenic photosynthesis. Its reaction center chromophores, where the charge transfer cascade is initiated, are arranged symmetrically along the D1 and D2 core polypeptides and comprise four chlorophyll (PD1, PD2, ChlD1, ChlD2) and two pheophytin molecules (PheoD1 and PheoD2). Evolution favored productive electron transfer only via the D1 branch, with the precise nature of primary excitation and the factors that control asymmetric charge transfer remaining under investigation. Here we present a detailed atomistic description for both. We combine large-scale simulations of membrane-embedded PSII with high-level quantum-mechanics/molecular-mechanics (QM/MM) calculations of individual and coupled reaction center chromophores to describe reaction center excited states. We employ both range-separated time-dependent density functional theory and the recently developed domain based local pair natural orbital (DLPNO) implementation of the similarity transformed equation of motion coupled cluster theory with single and double excitations (STEOM-CCSD), the first coupled cluster QM/MM calculations of the reaction center. We find that the protein matrix is exclusively responsible for both transverse (chlorophylls versus pheophytins) and lateral (D1 versus D2 branch) excitation asymmetry, making ChlD1 the chromophore with the lowest site energy. Multipigment calculations show that the protein matrix renders the ChlD1 → PheoD1 charge-transfer the lowest energy excitation globally within the reaction center, lower than any pigment-centered local excitation. Remarkably, no low-energy charge transfer states are located within the “special pair” PD1–PD2, which is therefore excluded as the site of initial charge separation in PSII. Finally, molecular dynamics simulations suggest that modulation of the electrostatic environment due to protein conformational flexibility enables direct excitation of low-lying charge transfer states by far-red light.


INTRODUCTION
Photosystem II (PSII) is a dimeric multisubunit protein− pigment complex responsible for the light-driven oxidation of water into molecular oxygen and for the supply of reducing equivalents in oxygenic photosynthesis. 1−7 Excitation-induced charge separation and the early steps of the electron transfer cascade take place within a cluster of six chlorin molecules known as the reaction center (RC). The RC comprises four chlorophylls (typically chlorophyll a), the central P D1 and P D2 "special" pair flanked by chlorophylls Chl D1 and Chl D2 , and two pheophytin molecules, Pheo D1 and Pheo D2 . The RC chromophores are arranged in a symmetric fashion along the D1 and D2 protein subunits of PSII ( Figure 1) that are highly conserved across oxygenic photosynthetic organisms. The RC of PSII receives excitation energy from integral chlorophyllrich polypeptides (CP43 and CP47) and from external lightharvesting antenna systems that vary among different species.
Following excitation and charge separation within the RC chromophores, the radical anion is localized within 0.3−3 ps 8−12 on Pheo D1 and the hole is distributed over the P D1 /P D2 pair. 13−15 The resulting radical cation, known as P680 + , is the strongest oxidant in biology. With an estimated reduction potential of 1.1−1.3 V, P680 + is able to drive the oxidation of water at the oxygen-evolving complex (OEC) via a redox active tyrosine residue (Tyr161 or Y Z ) that interfaces the two sites. A distinctive feature of PSII is the utilization of the D1 branch, which also harbors the OEC, for electron transfer following productive charge separation. On the acceptor side the negative charge proceeds from Pheo D1 to plastoquinone Q A and finally to the terminal mobile electron acceptor plastoquinone Q B . 16 Key questions include the nature and localization of initial excited states, the nature and energetics of charge-transfer (CT) excited states that may lead to productive charge separation, and the factors that determine the asymmetry of RC chromophores and directionality of electron transfer. 13,17,18 The multimer model 19 assumes that the local excitation energies of all pigments are similar, thus favoring delocalized excitation, but other studies supported models where the chromophore energies are distinct. 12,20−23 The P D1 −P D2 pair is a prime candidate for the initial excitation event in analogy to the "special pair" of bacterial reaction centers. 24,25 One of two major charge-transfer pathways 26,27 considered includes excitation and charge separation within this pair: [P D1 P D2 ]* → P D1 − P D2 + , and the negative charge is subsequently transferred to Chl D1 , i.e., [P D1 P D2 ] + Chl D1 − , 28−30 before proceeding to Pheo D1 . However, the prevailing view is that Chl D1 is the most red-shifted pigment and hence may function as the primary electron donor 11,21,23,31,32 (uncertainties remain regarding the low-energy excited states of pheophytins 23,33−35 ). The second pathway is thus described as [Chl D1 Pheo D1 ]* → Chl D1 + Pheo D1 − . It has been proposed that room temperature structural perturbation in the protein can induce switching between the different pathways. 26 Both would eventually lead to the same P D1 + Pheo D1 − chargeseparated state, but the electronic nature of the excitation and all underlying events are fundamentally distinct.
Stark spectroscopy studies suggested the presence of mixed local excitation−charge-transfer excitation in the active D1 branch, 36,37 while Styring and co-workers 38−40 proposed the presence of low-energy CT states responsible for far-red charge-separation and furthermore that the photochemistry of PSII is wavelength-dependent. 39 An important fact is that although the working threshold of PSII is typically considered to be 680 nm, charge separation in the reaction center can be initiated with far-red light (700−780 nm) either by direct excitation of a low-lying charge transfer state in Chl a RCs or by excitation of far-red chlorophylls (Chl d and f). 38,39,41−50 These observations highlight the significance of low-energy charge-transfer states for RC function. 17 Remarkably, no quantum chemical study has so far identified interpigment CT excited states low enough in energy to be consistent with these observations.
It is useful to keep in mind that experimental studies are typically restricted to nonphysiological and perhaps nonfunctional PSII preparations that may yield varying observations depending on the type of preparation and conditions used. Even disregarding light-harvesting antennae, a PSII monomer comprises more than 20 proteins and dozens of chlorophylls. Core complex preparations (PSII-CC) reduce this complexity by maintaining only the D1, D2, Cytochrome b 559 , CP43, and CP47 proteins, but the study of RC would still be challenging due to spectral congestion by the core protein chlorophylls, 27 therefore most experimental work involves PSII "RC complexes" (PSII-RCC) that maintain D1, D2, and Cytochrome b 559 . 51−55 These are considered as a minimal structural scaffold for studying the RC, 26,27,33,56−58 but their actual structure and the extent to which they represent the physiological system are debatable. 48,59−61 Computational studies can potentially assist in bridging the gap between observations made on such preparations and the properties of physiological PSII.
Theoretical studies 14,15,23,62−82 of photosynthetic reaction centers are challenging due to the size and complexity of the system. Beyond site energies, detailed electronic structure analysis of excited states that may be delocalized among different chromophores necessitates the use of quantum chemical approaches. A landmark study by Frankcombe used time-dependent density functional theory (TD-DFT) for the excited states of all PSII RC chromophores in the absence of the protein environment. 62 This and subsequent similar studies find neither asymmetry in local excitations along the D1 and D2 branches, nor low-lying CT states that could be related to the charge-separation function of the RC. 62,65 Taking the protein matrix into account with a combined quantum− mechanics/molecular−mechanics approach (QM/MM) 83−85 is obviously necessary. However, this approach must not be viewed as a mere technical extension that can automatically deliver good results. Four distinct methodological requirements must be met simultaneously and successfully to ensure a meaningful and reliable outcome. These are (i) explicit atomistic representation of the complete protein matrix, potentially with consideration of conformational dynamics, (ii) high-quality quantum chemical geometric definition of the chromophores 86−89 as opposed to the direct use of crystallographic coordinates, (iii) reliable excited state calculations of single and coupled chromophores, because site energies alone are insufficient to address the electronic nature of multichromophore excitations, and (iv) high-level quantum chemical methods that ensure the correct response of excited state energetics to protein electrostatics and, above all, provide a reliable description of the nature and energetics of interpigment charge-transfer excited states. The present study addresses for the first time all of the above requirements in a definitive way, aiming to provide reliable quantitative insights into the excitation profile of the reaction center of PSII.
A large-scale model of an entire membrane-embedded PSII complex is used as the basis for multiscale QM/MM modeling on geometry-optimized pigments to uncover the influence of the protein environment on the excitation profile of reaction center chromophores. Highly accurate quantum chemical descriptions of both local and distributed excitations among pairs of chromophores are obtained by long-range corrected time-dependent DFT as well as by coupled cluster theory at a level employed for the first time in such simulations, namely the similarity transformed equation of motion coupled cluster theory with single and double excitations (STEOM-CCSD). Our results provide a complete view of the nature and energetics of local and, most importantly, of charge-transfer excitations among different chromophore pairs. The results explain the origin of the dual type of asymmetry in the RC and identify the pigment pair responsible for the primary charge transfer excitation. In combination with insights from molecular dynamics simulations we determine the static and dynamic factors that control the excitation profile of the reaction center and enable charge separation to occur with farred light.

METHODOLOGY
2.1. Preparation of Models. The lipid-bilayer embedded model of Photosystem II is based on the high-resolution dimeric crystal structure (1.9 Å) of Thermosynechococcus vulcanus (PDB ID: 3WU2). 5 In the present work, we have used one of the monomers to build the entire system. Missing structural elements were completed to reproduce the physiological intactness of the complex. All crystallographic water associated with the monomer was retained and additional water molecules were added using the 3D-RISM technique 90−94 (three-dimensional reference interaction site model) to achieve a physiological hydration state of the protein complex. The PSII monomer was embedded inside a POPC (1-palmitoyl-2-oleoylsn-glycero-3-phosphocholine) lipid bilayer of dimension 176 × 176 Å 2 using Packmol-Memgen. 95,96 A total of 784 POPC molecules were added in the upper and lower leaflets of the trans-membrane region. The membrane-embedded protein complex was placed inside a water box. Appropriate amounts of Na + and Cl − ions were added to neutralize the system and maintain a physiological concentration of 0.15 M. The complete dimensions of the system were 176 × 176 × 160 Å 3 and it consists of 512 341 atoms. (Figure 2).
The electrostatic charges for all the cofactors were computed based on the MK-RESP (Merz−Kollman Restrained Electrostatic Potential) methodology. 97,98 For the organic cofactors, first the hydrogens atoms were optimized at the B3LYP/Def2-SVP level 99,100 and then singlepoint calculations were performed at the HF/6-31G* level of theory 97,101,102 using the ORCA program 103 and RESP fitting of the charges was performed using the Multiwfn code. 104 A bonded model is employed for the computation of the RESP charges on the OEC (Mn 4 CaO 5 -Oxygen Evolving Complex) and NHI (nonheme iron) sites, as a first step, a small cluster model is built around metal sites including the side chain of the residues which form a direct coordination with the metal site. The OEC is modeled in its S 1 state of the Kok−Joliot cycle, i.e., the oxidation states are Mn1(III)− Mn2(IV)−Mn3(IV)−Mn4(III) and involved ligands are Asp170, Glu354, Ala344, Asp342, Glu189, His332, Glu333, and four H 2 O molecules. Similarly, the NHI site is modeled as Fe(II) with the ligands HCO 3 − , His214, His268, His215, and His272. These models were first optimized at B3LYP/Def2-TZVP and then RESP fitting is performed at B3LYP/6-31G* level. More importantly, we restrained the charge of the backbone atoms of the residues according to the original AMBER force field 102 as such a procedure is known to produce better backbone dynamics during the simulation. 105 The RESP charges of the chlorophylls and the heme iron site were calculated in a similar fashion. The chlorophylls and the heme-iron are ligated axially to amino acids and water molecules, wherever applicable. For example, P D1 and P D2 of the reaction center are axially ligated to histidine residues and Chl D1 and Chl D2 are axially ligated to a single water molecule. Similarly, both heme sites are bound axially with two histidine residues.
Bonded parameters for the chlorophyll a, 106 heme, 107 and nonheme iron site 108 were obtained from the literature. Custom bonded parameters were defined for the OEC based on the study by Guerra et al. 109 Parameters for the standard protein residues were described using the Amber14SB force field. 110 The TIP3P model 111 was chosen for the water. The bonded parameters for the organic cofactors were described with GAFF2. 112 The LIPID17 force field 113,114 was chosen to describe our POPC bilayer. The nonbonded parameters for the metal ions were based on their respective oxidation states using data sets 105,115,116 available for the TIP3P water model. For Na + and Cl − , we used the Joung−Cheatham parameters compatible with the TIP3P water model. 117,118 2.2. Classical Molecular Dynamics Simulations. The complete system was minimized systematically to remove unfavorable geometric clashes. During the equilibration phase, the system is slowly heated from 10 to 100 K during 5 ps in the NVT ensemble. In the next step, the temperature is slowly increased from 100 to 303 K in the NPT ensemble, while maintaining the positional restraints (20 kcal mol −1 Å −2 ) on the C α atoms of amino acids. The temperature during this procedure is controlled using the Langevin dynamics 119 with a collision frequency of 5 ps −1 . We released the restraints on C α atoms in a controlled fashion (2 kcal mol −1 Å −2 /400 ps) and subsequently invoked the Monte Carlo/Molecular Dynamics (MC/ MD) module 120 for a controlled hydration and dehydration from bulk to protein, and vice versa. The system was then further equilibrated for another 63 ns in the NPT ensemble to properly equilibrate the lipid bilayer. Thereafter, we initiated the production simulation for 12 ns in the NPT ensemble using the collision frequency (Langevin dynamics) of 1 ps −1 . The pressure was regulated anisotropically using the Berendsen barostat 121 with a relaxation time of 2 ps and maintained at 1 bar. Particle Mesh Ewald (PME) 122 approach is used to treat all electrostatic interactions with a 10 Å cutoff. The SHAKE algorithm 123 was used to constrain bonds involving hydrogen atoms, which allowed us to use a time step of 2 fs. The frames were saved every 2 ps. Minimizations were performed using the CPU version while equilibration and production simulations were performed using the GPU version of the pmemd engine 124−126 of the AMBER18 package. 127,128 2.3. QM/MM Protocol. For the purposes of the present work we extracted five snapshots of the complete system from the classical molecular dynamics to be used in the QM/MM calculations. The first snapshot, which represents a structural configuration that is close to the crystal structure of PSII, is derived from the early equilibration procedure, where we clustered 129 a series of frames from the trajectory using the CPPTRAJ 130,131 module of AmberTools19. The other four snapshots were derived from the production run, that is, from an unbiasedly evolved, completely hydrated, and thoroughly equilibrated system. These four snapshots are equally spaced from each other, i.e., captured with an interval of 4.0 ns, which ensures that these snapshots are structurally uncorrelated. Inspection of the overall protein structure overlay ( Figure S1 of the Supporting Information, SI) of these snapshots confirms that the structural configuration of the protein in these snapshots is distinct. For the QM/MM setup we considered the entire PSII monomer and a total of 8000 water molecules, which includes all the waters presents in the protein cavity, various channels, and ∼7 Å bulk-region around the protein. In order to keep the system neutral, we maintained the required amount of Na + ions at their equilibrated positions. The final QM/MM system contains a total of 76 056 atoms.
The QM/MM calculations were performed with ChemShell 3.7, 132−134 where the in-built DL-POLY 135 was used for MM computations, whereas ORCA 103 was the QM engine. Our calculations are based on the electrostatic embedding technique. Covalent bonds were cut using the hydrogen link atom approach. The charge-shift method implemented in ChemShell is used to avoid overpolarization of the QM region by the MM region. For a given snapshot, we carried individual QM/MM geometry optimization calculations of Chl D1 , Chl D2 , Pheo D1 , and Pheo D2 , while the P D1 /P D2 pair is treated as a single QM unit in QM/MM optimizations due to the close proximity of the two chlorophylls. With such a procedure, we retain the symmetry breaking feature of the vinyl moiety in P D1 / P D2 , which is otherwise poorly represented by the MM force field. Similar concerns were raised by Mennucci and co-workers, 136 who found overstabilization of the acetyl group of Bacteriochlorophyll a by the MM force-field. In addition, with this approach both chromophores of the P D1 /P D2 pair retain their macrocyclic ring curvature, important for the position of the Q-band, 137,138 induced by their spatial proximity. Therefore, the geometry for all the single point calculations on the individual P D1 and P D2 is derived directly from the pair-optimized geometry.
During all QM/MM geometry optimizations the model is divided into two parts, active and static. The active region is defined by a QM subregion and an active MM region, whereas the static region remains fixed in the optimization procedure and can only influence the active region through the electrostatic effect of the point charges. In our cases where single chromophore is studied, we chose all residues and cofactors which are ca. 13 Å around the QM region (approximately 1300−1600 atoms). A larger active MM region (15 Å) was chosen for the optimization of the special pair (ca. 2650 atoms). For all the geometry optimizations, we employed the PBE functional 139 and Def2-TVZP basis set, 140 along with D3(BJ) dispersion corrections. 141,142 To speed up the calculations of Coulomb integrals we employed the resolution of identity (RI) approximation, 143,144 in combination with Weigend's universal Def2/J auxiliary basis set. 145 We used tight convergence criteria and a higher than default integration grid (grid6 in ORCA nomenclature) for all optimizations.
2.4. Computation of Excitation Energies. The computation of vertical excitation energies (VEE) in this work is performed on the QM/MM optimized geometries of the chromophores. The explicit effect of the entire PSII monomer environment on the electronic properties of chromophores is included through point-charges, handled by the orca_pc module which adds the contribution of the charges in the one-electron matrices and nuclear repulsion. We elucidated the effect of specific protein components on the excitation energies by switching off the corresponding point charges. Full Time-Dependent Density Functional Theory (TD-DFT), i.e., without the Tamm−Dancoff approximation, was used to compute the excited state properties using the range-separated ωB97X-D3(BJ) 146−148 functional along with the Def2-TZVP basis set. The choice of functional was guided by its highly successful application in related studies of biological chromophores 149 and its established superior performance for the treatment of charge-transfer states with TD-DFT. 150 The appropriateness of the functional is additionally confirmed in the present work via direct comparisons with coupled cluster theory. The RIJCOSX approximation 151,152 is used to speed up the calculations. In addition, very tight SCF convergence criteria were set along with higher integration grids (grid6 and gridX7 in ORCA nomenclature) in all cases considered in this work. We have computed 20 roots in all cases. We have also performed calculations with the long-range corrected LC-BLYP 153 functional, with the same settings as described above. In addition to individual chromophores we computed the excited states of groups of adjacent RC chromophores. The selection of pairs and trimers is based on proximity and on existing hypotheses about the primary excitation and charge separation. 26,29 Independently optimized geometries of individual chromophores were used to set up the oligomeric assemblies.
Despite using the best available functionals, known problems in the performance of TD-DFT in general and for charge-transfer states in particular, 154−162 indicate that independent validation beyond DFT is necessary to build confidence in the results, 88,163,164 especially when excited state properties of coupled chromophores are considered. Thus, both in order to obtain results that overcome potential limitations of TD-DFT and to further probe the nature of the excited states of the chromophores with an orthogonal methodology, we employed, for the first time in such a large-scale simulation, the domain based local pair natural orbital (DLPNO) implementation of the similarity transformed equation of motion coupled cluster theory with single and double excitations, STEOM-CCSD. 164−170 A recent study showed this method to provide a highly accurate description of all features of the gas-phase absorption spectrum of chlorophyll a, 171 while another study 172 showed that it performs exceptionally well for charge-transfer states, on par with coupled cluster methods that include triple excitations, in contrast to more approximate approaches such as CC2. For the coupled cluster QM/MM calculations, the first of their kind to be conducted on a photosynthetic reaction center, we computed in total 6 roots for each chromophore using the Def2-TZVP(-f) basis set for all atoms. When applied to pairs of chromophores, to maintain feasibility of the coupled cluster calculations the chlorophyll and pheophytin models were truncated by omitting macrocyclic ring substituents. With this procedure, we keep the number of basis functions within the bounds of computational feasibility while maintaining the chemical information regarding the nature and energetic order of excited states in RC chromophore pairs. For the DLPNO-STEOM-CCSD calculations of chromophore assemblies we have resorted to the Def2-SVP basis set and computed 20 roots in total. The RIJCOSX approximation is used to speed up the calculations throughout. "TightPNO" settings were applied for the DLPNO calculations. The T CutPNOsingles cutoff was set to 6.6 × 10 −10 and the active space selection keywords "Othresh" and "Vthresh" were set to 5.0 × 10 −3 .

RESULTS AND DISCUSSION
3.1. Structural Aspects of Pigment−Protein Interactions. The intrinsic photophysical properties of the chromophores are engineered in the protein matrix for efficient light harvesting. In Photosystem II various domains of multichromophoric systems exist, such as internal antenna systems (CP43 and CP47) and the reaction center. The chromophores in the RC, i.e., four chlorophyll a and two pheophytin a molecules, are symmetrically placed along the D1 (344 residues) and D2 (342 residues) polypeptide chains. These chromophores are situated deep inside the transmembrane region of PSII, and far away from the solvent exposed stromal and lumenal regions. Due to differences in the amino acid sequence of D1 and D2, 14 and overall PSII structural organization, the chromophores may experience distinct strain and electrostatic effects.
The chlorophylls of the P D1 /P D2 pair are weakly stacked in the middle of the PSII monomer and both of them are axially ligated with histidine residues, His198 for P D1 and His197 for Journal of the American Chemical Society pubs.acs.org/JACS Article P D2 (Figure 3). The accessory chlorophylls Chl D1 and Chl D2 are axially ligated with a single water molecule, while another water molecule is present in the PSII crystal structure that hydrogen-bonds simultaneously with the axially ligated water and the C-13 2 −COOCH 3 substituent. However, important differences are found in the second coordination sphere of Chl D1 and Chl D2 . The axially ligated water forms further hydrogen bond with Thr179 in case of Chl D1 , but the axial water of Chl D2 does not participate in any further hydrogen bonding and a hydrophobic Ile178 is found instead in its vicinity. Other important structural differences exist in the peripheral region, where the C-13 1 keto group of the macrocyclic ring is hydrogen bonded with a single water molecule. This water molecule is also found to be highly conserved across many high-resolution crystal structures. Interestingly, a recent experimental study using absorption spectroscopy and Resonance Raman techniques by Robert and co-workers 173 showed that hydrogen bonds to the keto group of Chl a can fine-tune the absorption properties of the LHCII (Light Harvesting Complex) antenna system. Another investigation by Collini and co-workers 174 showed how such water-mediated interaction with the peripheral substituents can influence the overall conjugation in Chl a/b found in the water-soluble chlorophyll protein (WSCP). A small number of differences also exist in the immediate environment of the Pheo D1 and Pheo D2 . The C-13 1 keto group of Pheo D1 is hydrogen bonded to the Gln130, whereas, the C-13 1 keto group of Pheo D2 can establish two hydrogen bonds simultaneously with Gln129 and Asn142. The C-13 2 − COOCH 3 and C-17 3 −COOR groups of Pheo D1 hydrogenbonds to two tyrosines, Tyr147 and Tyr126. However, C-13 2 − COOCH 3 and C-17 3 −COOR groups of Pheo D2 are surrounded by hydrophobic residues (Phe255, Phe125, and Phe146). Analysis of our molecular dynamics trajectories indicates that the D1 and D2 chains are extremely stable throughout the simulations compared to the dynamic evolution of the complete system (Figure 4). This strongly suggests that the D1 and D2 proteins do not undergo any large-scale conformational change in the PSII core complexes within the time scale of the simulation. The stability of the D1 and D2 poly peptides ensures that the relative orientations of the chromophores stay essentially intact. In terms of individual chromophores ( Figures 5 and S2−S7), we observe that the distance between P D1 and P D2 spans only a very narrow range around the average of ca. 3.5 Å along the simulation. Contrasting behavior is seen between Chl D1 and Chl D2 with respect to the hydrogen-bonding interactions of their C-13 1 keto groups with the nearby water molecule. This hydrogen bond appears quite constrained in the case of Chl D1 but samples a wider range of distances in the case of Chl D2 . Finally, we observe that the protein environment around Pheo D2 is highly dynamic, with significant degree of rotations for the side-chains of Gln129 and Asn142. In contrast, the hydrogen bonding interaction of Pheo D1 with Gln130 remains stable throughout the MD simulation and no large fluctuations in side-chain conformations are observed. Overall, production MD simulations over 12 ns strongly indicate that the environment of the D1 branch chromophores is more rigid and under tighter steric control of the protein matrix than the D2 branch.
3.2. Excitation Profiles of Individual Chromophores. The protein environment is known to influence the excited state properties of the chromophores in various ways, such as geometric strain, hydrogen bonding, and electrostatic effects. To understand the influence of these factors in a bottom-up fashion, first we computed the excited state properties of individual chromophores using their QM/MM optimized geometries but in the absence of any representation of the protein environment. These gas-phase TD-DFT calculations with the ωB97X-D3(BJ) functional lead to approximately the same S 0 → S 1 (Q y ) vertical excitation energies for all chromophores, within the narrow range of 1.920−1.943 eV (see Figure 6 and Tables S1−S4). Similar results were  Journal of the American Chemical Society pubs.acs.org/JACS Article obtained in past studies 62,65 where a uniform dielectric medium was used to mimic the protein environment. As the geometries of the pigments already incorporate the structural effects of the protein matrix, we conclude that the proteininduced strain on the macrocyclic rings is not responsible for inducing functional asymmetry in the localization of the lowest excited state in the RC. Subsequently, excited state calculations using the same QM/ MM optimized geometries were performed in a TD-DFT− QM/MM fashion, i.e., in the presence of the protein electrostatic environment (values in green in Figure 6). Two striking observations can be made: (a) all four chlorophyll a chromophores of the RC are red-shifted upon embedding in the protein environment, whereas both pheophytins are blueshifted (∼0.1 eV), and (b) Chl D1 is the most strongly redshifted pigment. Interestingly, the "special pair" P D1 −P D2 displays negligible internal asymmetry in terms of site energies either with or without the protein matrix. The identification of Chl D1 as the pigment with the lowest energy excited state supports previous interpretations that account for the effects of the protein matrix. [21][22][23]67,76,175 Different quantum chemical methods agree on the absence of noticeable asymmetry in calculations that omit the electrostatic effect of the protein, even though they provide different absolute values for the excitations, which is anticipated (see SI). What is important for the fundamental question of the emergence of excitation asymmetry is the reliable reproduction of the nature and extent of red or blue shift of the lowest excitation when the electrostatic effect of the protein matrix is included in the calculations. Comparison of the shifts produced by different methods (Table 1) shows that ωB97X-D3(BJ) and LC-BLYP produce very similar shifts, in line with those of DLPNO-STEOM-CCSD calculations. The coupled cluster results agree best with ωB97X-D3(BJ). The only significant difference between the two methods is that STEOM-CCSD predicts an even stronger blue shift for Pheo D2 .
It is evident that two types of protein matrix electrostatic asymmetry develop within the transmembrane region: (a) transverse asymmetry, which red-shifts the chlorophylls and blue-shifts the pheophytins, and (b) lateral asymmetry, which differentiates the D1 and D2 branches. This is also reflected in the electrostatic potential as experienced by the chromophores in the protein matrix. The potential is indeed distinct along D1 and D2, and the pheophytins reside in a relatively more positive electrostatic potential pocket (Figure 7). Therefore, the intrinsic protein electrostatic environment is the principal factor in modulating the distribution of excitation energies of RC chromophores and giving rise to both chlorophyll− pheophytin asymmetry (transversely) and D1−D2 branch asymmetry (laterally).   The asymmetry induced by protein electrostatics can be associated with differences between the D1 and D2 sequences, spatial proximity of charged residues and redox active cofactors, and the overall organization of extrinsic proteins, which differ in cyanobacteria, algae, and higher plants. Each chromophore experiences the intrinsic protein electrostatics differently because of their location and orientation with respect to the transmembrane region. Therefore, the asymmetry in the reaction center is not an intrinsic property related to the spatial arrangement of the chromophores. This explains why past computational studies that used at most a uniform dielectric to mimic the protein environment 62,65 reported essentially the same Q y energies for the individual RC chromophores but did not reproduce the excitation asymmetry that is a feature of the real system. It is noted that previous works 20,23,33 employing spectral modeling have assigned nearly equal Q y energies of chlorophylls and pheophytins. This is not in line with the present results and further work is required for a confident analysis of the optical spectrum. It cannot be excluded that a more elaborate theoretical treatment of the protein environment, for example with inclusion of polarization effects, may provide a better balance between the spectral fitting and the quantum chemical results.
Identification of the key protein components that give rise to asymmetry in the reaction center is important for our overall understanding of unidirectional electron transfer. Protein electrostatics are responsible for red-shifting Chl D1 by 0.064 eV (516 cm −1 ), which is by far the most drastic protein matrix effect on a single chromophore. The major contributors to this red-shift include Met172 (76 cm −1 ) and Phe158 (48 cm −1 ). Interestingly, the pseudosymmetric partners of these Chl D1 "red-shifters" are dif ferent on the Chl D2 side, which suggests that localized factors control the asymmetry around these chlorophylls. Additional contributions arise from P D1 (79 cm −1 ), and from the chloride ions in the vicinity of the OEC (61 cm −1 ). However, major contributors toward the 0.141 eV (1137 cm −1 ) blue shift of Pheo D1 were found to be mostly the closely lying amino acids and cofactors, including Tyr147 (115 cm −1 ), Pro150 (89 cm −1 ), Chl D1 (66 cm −1 ), Leu151 (52 cm −1 ), and Ile213 (28 cm −1 ).
It has been suggested that the OEC or residues that coordinate the Mn 4 CaO 5 cluster contribute to the red shift of Chl D1 . 67 Our calculations show that this is not the case and we attribute this to incorrect MM setup. Specifically, the use of integer charges for the Mn ions of the OEC in accordance to formal oxidation states results in exaggerated concentration of charge on the inorganic core and individually on its ligands. The physically motivated approach is to use distributed charges as in the present work, deducing them from RESP calculations that treat the Mn 4 CaO 5 cluster and all its covalently bonded ligands as a single chemical entity. This eliminates longrange Coulombic artifacts. Beyond the fundamental technical aspect, attribution of a major site energy determining role to the OEC is conceptually problematic because normal function of the reaction center is required for photoassembly of the OEC 177−179 and hence must be independent of its presence.
Although we pinpointed certain major contributors to the observed shifts, the total shifts in each case are not completely produced by a limited list of contributions, but instead an ever increasing number of residues and cofactors with ever smaller contributions can be found. Therefore, not only localized but also global electrostatics 23 play a key role to produce the total shift compared to the gas-phase result, which shows that the evolutionary optimization of the enzyme operates on all scales. The present description has a parallel in the photosynthetic bacterial reaction center, where the chromophores were also found to be embedded in a dielectrically asymmetric environment. 180,181 Some comments on methodological issues must be made at this point. First of all, the above observations underline the necessity of the QM/MM approach for obtaining meaningful results. We confirm that enlarging a gas-phase QM model by including several selected amino acid residues around each chromophore does not even begin to approximate the full QM/MM results. Another critical methodological issue relevant to studies that employ QM or QM/MM methods is the use of experimental (crystallographic) structures for the pigments. This choice is detrimental for quantum chemical approaches (alternative approaches 23,175,182 for studying site energies are largely unaffected) because experimental structures typically do not exhibit correct bond length alternation of conjugated systems, 86,87 often not even qualitatively. This is, however, the key geometric parameter that governs the electronic structure of the ground and excited states because it directly determines the nature and energetics of frontier orbitals. Therefore, the use of experimental geometries introduces a fundamental inconsistency between the quantum chemical approach and the structural model on which it operates, leading to randomization of results through uncontrolled errors. When this is coupled with the neglect of protein matrix electrostatics, the combined methodological deficiencies practically guarantee that the quantum chemical Journal of the American Chemical Society pubs.acs.org/JACS Article results are of little relevance to the real system (for example, a study 183 that satisfies none of these conditions finds the lowestenergy excitation of the RC to be localized on the Pheo D1 ). The same holds for the direct use of force-field (MM) geometries. 67 These considerations apply equally to the quantum chemical treatment of multiple RC chromophores, an even more delicate case because of the sensitivity of interpigment charge-transfer states on both the geometries of interacting chromophores and on protein matrix electrostatics.

Excitation
Profiles of Chromophore Pairs. Intermolecular charge-transfer (CT) excitations, i.e., those where the electron donor and acceptor are two different chromophores, are central in the function of reaction centers across photosynthetic organisms. 36,37,39 Although site energies already reveal a lot about the RC of PSII, understanding the initiating events of productive primary photoexcitation requires direct insight into the excitation profiles of multiple chromophores. To obtain this information we have systematically studied the excited states of pairs of chromophores in the D1 and D2 branches, i.e., P D1 −P D2 , P D1 −Chl D1 , P D2 −Chl D2 , Chl D1 −Pheo D1 , and Chl D2 -Pheo D2 .
The "special pair" holds a special status both in bacterial reaction centers 13,184 and in Photosystem II, 13,26,29 therefore we focus on this one first. TD-DFT calculations performed in the gas phase (using the QM/MM optimized geometry of the pair) reveal that the nature of the lowest excitation, at 1.917 eV, is a linear combination of local excitations (LE) on P D1 and P D2 individually. Analysis of the excitations in terms of natural transition orbitals (NTOs) is provided in Figure 8. The first excitation with CT character (root 5), P D1 − P D2 + , is situated considerably higher in energy, at 3.091 eV. In comparison to the gas-phase results, TD-DFT QM/MM computations performed with full account of protein electrostatics result in an overall red-shift of the lowest excited state, which is now predicted at 1.884 eV, however the nature of the excited state remains the same, i.e., superposition of local excitations. We note that this is what one would expect in a Frenkel excitation picture and that the second excited state (S 2 , at 1.911 eV, see also Table S5) is also a linear combination of the two local excited states. Protein electrostatics affect the "directionality" of the lowest charge transfer state (P D1 + P D2 − ) but do not stabilize it significantly compared to the gas-phase result (2.999 eV).
Next, we have investigated the symmetry-related P D1 -Chl D1 and P D2 -Chl D2 pairs along the active and inactive chains. The lowest excited state (1.905 eV) computed from gas-phase calculations on P D1 -Chl D1 is a linear combination of local excitations on Chl D1 and P D1 , whereas the first CT state (root Journal of the American Chemical Society pubs.acs.org/JACS Article 9, P D1 − Chl D1 + ) was found much higher in energy (3.639 eV). Protein electrostatics induce an overall red shift in the lowest excited state (1.843 eV) and character alteration of the lowest excited state, which becomes a local excitation on Chl D1 . In addition, the lowest CT state (root 5) changes in directionality (P D1 + Chl D1 − ) and is found much lower in energy, at 3.130 eV. Similar observations are made for the P D2 −Chl D2 pair, where in relation to the gas-phase results the protein electrostatics induces a red-shift of the lowest excited state (1.895 eV, LE on P D2 ), lowers the energy and alters the directionality of the first CT state (P D2 + Chl D2 − , root 5, 3.232 eV, see also Table S7). In the case of the Chl D2 -Pheo D2 pair (see Figure S8), the gas-phase calculations indicate that the lowest energy state (1.926 eV) is a local excitation on Chl D2 , whereas the lowest CT state (Chl D2 + Pheo D2 − , root 9) is located at 3.726 eV. Within the protein matrix we computed a slight red-shift of the lowest excited state (1.903 eV, LE on Chl D2 ) and a significant stabilization of the CT state (Chl D2 + Pheo D2 − , root 3) at 2.092 eV.
The most profound demonstration of protein matrix control is observed in the case of the Chl D1 −Pheo D1 pair. The nature of the lowest excited state at 1.905 eV obtained from the gasphase calculations on the Chl D1 −Pheo D1 is primarily a local excitation on the Chl D1 . The first CT state (root 9, Chl D1 + Pheo D1 − ) is much higher in energy (3.728 eV). In this case, however, protein electrostatics drastically reorganize the excitation profile of the pair. The lowest excited state is computed at 1.828 eV and has significant Chl D1 + Pheo D1 − charge transfer character ( Figure 8). Crucially, this particular CT state marks the lowest-energy excited state among all chromophore pairs of the PSII RC. In fact, not only the first but also the second excited state computed for this pair display significant or dominant CT character favoring excitation from Chl D1 to Pheo D1 . The picture obtained from the NTOs is consistent with the difference densities computed for these roots, as will also be shown for the additional MD snapshots discussed in the present work and which have pure Chl D1 + Pheo D1 − charge-transfer character (vide infra). The above results are in line with observations from Stark spectroscopy, 36,37 where the Chl D1 + Pheo D1 − CT state was found to be mixed with LE on Chl D1 . This mixing was proposed to be crucial for initiating charge separation in the RC. A similar observation was made by Valkunas and coworkers, 66 where a good fit of the Stark spectrum (using complex time-dependent Redfield theory) was obtained with inclusion of the Chl D1 + Pheo D1 − CT state, which was found to be more important in reproducing the Stark spectrum than the CT state originating from P D1 − P D2 + . Calculation of the excited states of the Chl D1 −Pheo D1 pair using the conductor-like polarizable continuum model (CPCM) 185,186 with a constant dielectric (ε = 4, appropriate for the transmembrane region) does not produce significant differences in the excitation spectrum of the chromophore pair compared to the gas-phase results. The lowest gas-phase CT state (S 9 , 3.728 eV) is only marginally stabilized by CPCM (S 8 , 3.592 eV). Overall the nature of all excited states remains essentially indistinguishable from the gas-phase TD-DFT spectrum, and hence the continuum dielectric approach is no substitute for the electrostatically asymmetric protein matrix.
Even by using one of the best available DFT methods, the fact that we predict the lowest excited state to have dominant CT character creates the need for further confirmation, because the correct prediction of CT character has been historically a challenge for TD-DFT. The only definitive way to achieve this is to go beyond DFT. Excited states of truncated pigment pairs computed by DLPNO-STEOM-CCSD (see discussion in the SI and Tables S15 and S16) confirm the nature of excited states obtained with ωB97X-D3(BJ), and therefore fully support the above conclusions beyond any conceivable uncertainty arising from the level of theory.
In conclusion, our QM/MM results on monomer and pair excitation energies confirm Chl D1 as the pigment with the lowest site energy and furthermore identify the lowest excitation of the RC as associated with a CT state of the Chl D1 −Pheo D1 pair.
3.4. Excited States of the P D1 −P D2 −Chl D1 Trimer. The results presented above firmly support a CT excited state of the type Chl D1 + Pheo D1 − as the lowest energy excited state among pairs of RC chromophores, and hence suggest that actual charge separation would occur accordingly within this pair. An alternative pathway mentioned in the introduction as one of the possibilities under discussion involves excitation of P D1 or the P D1 −P D2 pair with charge transfer to Chl D1 . To investigate this possibility we conducted the same type of QM/MM calculations with simultaneous inclusion of all three relevant chromophores in the QM region (P D1 , P D2 , and Chl D1 ). The TD-DFT QM/MM results presented and analyzed in terms of NTO compositions in Table S17 show that the lowest excited The lowest is a CT state within the P D1 −P D2 pair (P D1 + P D2 − , S 7 at 3.050 eV), while the first CT state with P D1 + Chl D1 − CT character is S 8 at 3.093 eV. These results are in line with those obtained for monomers and dimers. Therefore, the present data on the trimer exclude the possibility of an energetically accessible CT excited state within the RC that involves delocalization of negative charge onto Chl D1 , and strongly disfavor the hypothetical participation of an anionic Chl D1 species in native PSII charge separation.
3.5. Dynamic Control of Low-Energy Charge-Transfer States. In view of the key role of the electrostatic environment described above, it is interesting to investigate if the dynamic evolution of the protein conformation influences the excited state properties of the RC established so far within a single structural configuration of PSII. For this purpose, we performed the same set of excited state calculations, each with individually optimized QM/MM geometries, on structurally independent snapshots obtained from the molecular dynamics simulations. These snapshots were obtained from unbiased production simulations of the PSII−membrane complex, i.e., with no restraints or constraints. We chose four distinct structural configurations of the PSII−membrane complex with a consecutive interval of 4 ns. This has the advantage that the excited state properties are computed on uncorrelated protein configurations that are removed from the crystal structure minima and are properly hydrated and equilibrated with the surrounding environment.
The overall trend in the respective blue and red-shift of the individual RC chromophores and the relative ranking of the Q y excitation energies remains the same, i.e., Chl D1 has the lowest site energy (Figure 9). Focusing on chromophore pairs, we find that the lowest-lying excited state of the RC remains on the Chl D1 −Pheo D1 pair and retains its CT character (Chl D1 + Pheo D1 − ) irrespective of the dynamics of the protein. The above observations suggest that the nature of the intrinsic electric field of the protein matrix, and the resulting excitation asymmetry are essentially unperturbed by the conformational dynamics of the protein. Figure 10 depicts difference densities for the lowest excitation of the Chl D1 −Pheo D1 pair, which also demonstrate that the effect of the protein matrix is the same both in the crystallographic conformation of the protein and in the selected MD snapshots. Overall, we conclude that asymmetry does not arise as a result of the random conformational fluctuations. However, our findings indicate that the conformational flexibility of the PSII complex plays another important role.
Protein dynamics are seen to affect chromophore pairs in different ways. The conformational flexibility of the protein has little impact on the excited state properties of the inactive branch or of the P D1 −P D2 pair (Figures 9 and S9), whereas there is a significant impact on the Chl D1 −Pheo D1 pair of the active branch, where we observe high sensitivity in the energy of the first excited state. For the conformations of the protein studied here, we find the S 1 states with dominant (snapshot 1) to pure (all snapshots along the MD) Chl D1 + Pheo D1 − CT character in the range 1.828 eV (678 nm) to 1.595 (777 nm). Although the absolute computed values are not suggested to be "exact", since there is a dependence of the absolute values on the choice of QM method, it is interesting to note that this is beyond the nominal threshold for oxygenic photosynthesis of 680 nm (1.82 eV).
We stress that the number of snapshots used here is very small and a considerably more extensive sampling of the MD trajectory would be needed for quantitative analysis. Nevertheless, the present observations serve adequately as demonstration of principle, namely that the intrinsic electrostatic environment and flexibility of the protein are responsible for enabling access to low-energy charge-transfer states. This implies that protein matrix dynamics can push the red limit of Journal of the American Chemical Society pubs.acs.org/JACS Article oxygenic photosynthesis even in species that do not benefit from alternative types of chlorophyll (d or f). In this sense, specific variants of core PSII proteins available to different organisms might be utilized to adjust the red limit in response to environmental conditions not only by presenting alternate localized electrostatic contributors to critical pigments but also by favoring different distribution of global protein conformations. Therefore, electrostatic control by the dynamically evolving protein matrix must be considered equally important to the intrinsic absorption properties of participating chromophores in determining the red limit of photosynthesis. 3.6. Implications for Charge Separation Pathways. The computational results presented above form a solid basis for exploring the physiological function of the reaction center in PSII and connecting with various experimental observations. Studies on sunflower and bean leaves, 43 spinach, 38 green algae, 46 and cyanobacteria 47 showed that the known threshold for charge-separation in oxygenic photosynthesis can be pushed to the far-red region. 38,40, 45 Pettai et al. 43,44 reported that higher plants can evolve oxygen using wavelengths as long as 780 nm. Hughes et al. 45 showed that charge separation in PSII can be induced with light of 690−730 nm wavelength (1.7−1.8 eV) at 1.7 K. Similarly, Styring and co-workers 38,40 documented generation of the cation radical in RC using farred light (up to 750 nm), however a decrease in the overall charge-separated centers was observed with increasing wavelength. Furthermore, it was suggested that a significant population of the Chl D1 cation radical is trapped at cryogenic conditions (5 K) and that Tyr161 (Y Z ) is the preferred electron donor in this case over the Cyt-b 559 /Chl Z /Car D2 pathway (see Figure 1) in far-red light at 5 K. This is fully consistent with our results, which show the "dark" low-lying CT states of the Chl D1 −Pheo D1 pair to be in the red and farred region. This implies that after charge separation within this pair takes place, the cation radical would initially be formed on Chl D1 . The hole would subsequently migrate to P D1 under physiological conditions. Our results are also in very good agreement with the suggestion that Y Z becomes the preferred electron donor under far-red light at 5 K. The migration of the Chl D1 hole to the special pair would require a reorganization of the protein environment that is probably inhibited at cryogenic temperature, resulting in Y Z becoming an electron donor to Chl D1 . Further EPR experiments 39 showed that charge separation upon light excitation is wavelength-dependent, leading to the hypothesis that P D1 is excited with visible light (532 nm), whereas Chl D1 is excited with far-red light. 39 Our results on the P D1 −P D2 −Chl D1 trimer do not support the presence of any energetically accessible CT state that could lead directly to productive charge separation along the D1 branch, but the conditions under which the special pair may function as the primary donor will need further studies to be clarified.
The present results are also of relevance to understanding RC function in organisms that employ alternative types of chlorophyll. Specifically, some cyanobacteria acclimatized to far-red light synthesize Chl d (Acaryochloris marina) 187 or Chl f (Chroococcidiopsis thermalis) 49,188 (absorption maxima at 710 and 750 nm, respectively) with alterations of chlorophyll pigments in the light harvesting antennae and possibly in the reaction centers themselves. The existence and location of Chl d pigments in the RC is a subject of debate. 189,190 In the case of C. thermalis Nurnberg et al. 49 assigned Chl f to be the Chl D1 position, while recent work by Judd et al. 191 demonstrated the likelihood of P D2 being occupied by Chl f. Future work will evaluate these possibilities using the QM/MM approach presented in this study.
Finally, our results indicate that protein conformational flexibility would play a critical role in charge separation and subsequent oxygen evolution with far-red light. Due to the high dependence of the far-red absorption capability of Chl D1 − Pheo D1 on conformational dynamics, only a fraction of PSII centers 38 could lead to formation of the charge-separated state in far-red. Since four "productive" charge separation events are required to make one O 2 molecule, there is reduced likelihood that the OEC can advance regularly under these conditions, i.e., that subsequent far-red charge separations can occur ontime to advance the catalytic cycle outcompeting recombination. The fact that low O 2 evolution is observed in cyanobacterial and higher-plant PSII excited with far-red light 43,49 is consistent with this scenario.

CONCLUSIONS
We presented large scale MM−MD and QM/MM results on a complete membrane embedded PSII monomer, focusing on excitation energies of single and paired reaction center chromophores. The quantum mechanical level of theory in the QM/MM calculation of excitation energies include longrange corrected TD-DFT and the DLPNO implementation of the similarity transformed equation of motion coupled cluster theory with single and double excitations, STEOM-CCSD. The approach presented here serves as a reference for future studies of the reaction center. Any deviations of past reports from the results of the present work for either single or multiple chromophores can be directly traced to neglect of one or more of the methodological pillars defined in the present study.
Our results demonstrate explicitly that the excitation asymmetry in the reaction center of PSII is not an intrinsic property of RC chromophores and does not originate from their distinct geometric distortion or coordination. Asymmetry is not observed, and cannot be understood, in the absence of the protein environment. It arises exclusively through the electrostatic effect of the protein matrix. We demonstrated that the electrostatic field of the protein acts by shifting the intrinsic site energies of the chlorophylls and pheophytins in opposite directions. Red-shifting of chlorophylls versus blue-shifting of pheophytins creates transverse asymmetry with respect to the membrane normal. Preferential lowering of site energies and most importantly of charge-transfer excited states along the D1 side of the Photosystem II creates lateral asymmetry. In the presence of the protein matrix the pigment with the lowest site energy is Chl D1 and the lowest excited state within the reaction center is a state of charge transfer character localized at the Chl D1 −Pheo D1 pair. Therefore, the present results support assigning this pair as the site of initial charge separation in PSII. The central P D1 −P D2 chlorophylls do not present lowenergy charge transfer excitations internally, while the possibility of charge transfer excitation from this pair to Chl D1 is excluded. Finally, protein dynamics have only weak influence on the localization of low-energy excitations, but enable charge transfer excitations within the Chl D1 −Pheo D1 pair to occur with far-red light. The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/jacs.0c08526.
Figures S1−S11 and Tables S1−S20, overview of site energies of the PSII RC chromophores; excited state properties of interacting pigment pairs with ωB97X-D3(BJ); excited state properties of interacting pigment pairs with LC-BLYP; comparison of TDDFT and DLPNO-STEOM-CCSD; excited states of the P D1 − P D2 −Chl D1 trimer; low-lying excited states of interacting chromophores from additional protein configurations; D1 and D2 sequence alignment from different photosynthetic organisms; Cartesian coordinates of optimized chromophores; and additional references (PDF)