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 PD1 and PD2 “special” pair flanked by chlorophylls ChlD1 and ChlD2, and two pheophytin molecules, PheoD1 and PheoD2. The RC chromophores are arranged in a symmetric fashion along the D1 and D2 protein subunits of PSII (Figure ) that are highly conserved across oxygenic photosynthetic organisms. The RC of PSII receives excitation energy from integral chlorophyll-rich polypeptides (CP43 and CP47) and from external light-harvesting antenna systems that vary among different species.

Figure 1

(a) Arrangement of the photosystem II dimer viewed from the stromal side of the membrane; marked regions are the core-antenna chlorophyll proteins CP43 and CP47, and the reaction center. (b) Reaction center chromophores and selected additional components within the D1 and D2 polypeptides, with a scheme indicating the flow of electrons within PSII.

Following excitation and charge separation within the RC chromophores, the radical anion is localized within 0.3–3 ps[8−12] on PheoD1 and the hole is distributed over the PD1/PD2 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 YZ) 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 PheoD1 to plastoquinone QA and finally to the terminal mobile electron acceptor plastoquinone QB.[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 PD1–PD2 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: [PD1PD2]* → PD1–PD2+, and the negative charge is subsequently transferred to ChlD1, i.e., [PD1PD2]+ChlD1–,[28−30] before proceeding to PheoD1. However, the prevailing view is that ChlD1 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 [ChlD1PheoD1]* → ChlD1+PheoD1–. 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 PD1+PheoD1– charge-separated 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 b559, 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 b559.[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 far-red light.

Methodology

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-oleoyl-sn-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 ).

Figure 2

Side and top (stromal) view of the molecular mechanics model of the PSII monomer embedded in an equilibrated POPC lipid bilayer. The simulation box has dimensions of 176 × 176 × 160 Å3.

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 single-point 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 (Mn4CaO5- 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 S1 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 H2O molecules. Similarly, the NHI site is modeled as Fe(II) with the ligands HCO3–, 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, PD1 and PD2 of the reaction center are axially ligated to histidine residues and ChlD1 and ChlD2 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]

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]

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 ChlD1, ChlD2, PheoD1, and PheoD2, while the PD1/PD2 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 PD1/PD2, 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 PD1/PD2 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 PD1 and PD2 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.

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 TCutPNOsingles 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

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 PD1/PD2 pair are weakly stacked in the middle of the PSII monomer and both of them are axially ligated with histidine residues, His198 for PD1 and His197 for PD2 (Figure ). The accessory chlorophylls ChlD1 and ChlD2 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-132–COOCH3 substituent. However, important differences are found in the second coordination sphere of ChlD1 and ChlD2. The axially ligated water forms further hydrogen bond with Thr179 in case of ChlD1, but the axial water of ChlD2 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-131 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 PheoD1 and PheoD2. The C-131 keto group of PheoD1 is hydrogen bonded to the Gln130, whereas, the C-131 keto group of PheoD2 can establish two hydrogen bonds simultaneously with Gln129 and Asn142. The C-132–COOCH3 and C-173–COOR groups of PheoD1 hydrogen-bonds to two tyrosines, Tyr147 and Tyr126. However, C-132–COOCH3 and C-173–COOR groups of PheoD2 are surrounded by hydrophobic residues (Phe255, Phe125, and Phe146).

Figure 3

Protein environment around chromophores associated with the reaction center in Photosystem II. Hydrogen bonding interactions are indicated with dashed lines.

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 ). 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 and S2–S7), we observe that the distance between PD1 and PD2 spans only a very narrow range around the average of ca. 3.5 Å along the simulation. Contrasting behavior is seen between ChlD1 and ChlD2 with respect to the hydrogen-bonding interactions of their C-131 keto groups with the nearby water molecule. This hydrogen bond appears quite constrained in the case of ChlD1 but samples a wider range of distances in the case of ChlD2. Finally, we observe that the protein environment around PheoD2 is highly dynamic, with significant degree of rotations for the side-chains of Gln129 and Asn142. In contrast, the hydrogen bonding interaction of PheoD1 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.

Figure 4

Time evolution of the root-mean-square deviation (RMSD in Å) of the Cα atoms of the complete PSII complex and of the D1 and D2 core polypeptides during the production MD run. Disordered regions of the protein are not considered in the RMSD calculations.

Figure 5

Time evolution of selected distances involving RC chromophores during the production MD run. Further details are provided in the SI.

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 S0 → S1 (Q) vertical excitation energies for all chromophores, within the narrow range of 1.920–1.943 eV (see Figure and Tables S1–S4). Similar results were 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 protein-induced strain on the macrocyclic rings is not responsible for inducing functional asymmetry in the localization of the lowest excited state in the RC.

Figure 6

Lowest-energy excitations (site energies), in eV, of PSII reaction center chromophores without (red) and with (green) the protein electrostatic environment. The calculations were performed with the ωB97X-D3(BJ) functional using QM/MM geometries.

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 ). 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 blue-shifted (∼0.1 eV), and (b) ChlD1 is the most strongly red-shifted pigment. Interestingly, the “special pair” PD1–PD2 displays negligible internal asymmetry in terms of site energies either with or without the protein matrix. The identification of ChlD1 as the pigment with the lowest energy excited state supports previous interpretations that account for the effects of the protein matrix.[21−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 ) 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 PheoD2.

Table 1

Comparison of Electrochromic Shifts (in eV) of Site Energies (Q) from the Gas-Phase to the Protein Matrix for Individual Reaction Center Chromophores, Computed with TD-DFT Using Two Different Range-Separated Density Functionals and with DLPNO-STEOM-CCSD

	ωB97X-D3(BJ)	LC-BLYP	DLPNO-STEOM-CCSD
PD1	–0.018	–0.016	–0.020
PD2	–0.023	–0.021	–0.015
ChlD1	–0.064	–0.062	–0.067
ChlD2	–0.027	–0.025	–0.025
PheoD1	+0.141	+0.136	+0.142
PheoD2	+0.145	+0.140	+0.177

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 ). 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).

Figure 7

Electrostatic potential experienced by the reaction center chromophores (in kT/e) inside the protein matrix. The calculations were performed using the APBS (Adaptive Poisson–Boltzmann Solver) program[176] on the crystal-structure-like configuration of the protein (snapshot 1).

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 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 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 ChlD1 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 ChlD1 “red-shifters” are different on the ChlD2 side, which suggests that localized factors control the asymmetry around these chlorophylls. Additional contributions arise from PD1 (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 PheoD1 were found to be mostly the closely lying amino acids and cofactors, including Tyr147 (115 cm–1), Pro150 (89 cm–1), ChlD1 (66 cm–1), Leu151 (52 cm–1), and Ile213 (28 cm–1).

It has been suggested that the OEC or residues that coordinate the Mn4CaO5 cluster contribute to the red shift of ChlD1.[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 Mn4CaO5 cluster and all its covalently bonded ligands as a single chemical entity. This eliminates long-range 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 results are of little relevance to the real system (for example, a study[183] that satisfies none of these conditions finds the lowest-energy excitation of the RC to be localized on the PheoD1). 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., PD1–PD2, PD1–ChlD1, PD2–ChlD2, ChlD1–PheoD1, and ChlD2-PheoD2.

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 PD1 and PD2 individually. Analysis of the excitations in terms of natural transition orbitals (NTOs) is provided in Figure . The first excitation with CT character (root 5), PD1–PD2+, 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 (S2, 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 (PD1+PD2–) but do not stabilize it significantly compared to the gas-phase result (2.999 eV).

Figure 8

A detailed description of the identity and nature of the lowest excited state (S1) and of the first root with significant charge transfer character for the PD1–PD2 and ChlD1–PheoD1 pairs in terms of Natural Transition Orbitals (NTOs) and relative contributions to a given excitation. In the case of ChlD1–PheoD1 both states depicted have significant CT character. Vertical excitation energies (in eV) and oscillator strengths (f) are provided for each state depicted (from ωB97X-D3(BJ) TD-DFT calculations), comparing the results in the absence (left) and in the presence (right) of the electrostatic effect of the complete PSII monomer.

Next, we have investigated the symmetry-related PD1-ChlD1 and PD2-ChlD2 pairs along the active and inactive chains. The lowest excited state (1.905 eV) computed from gas-phase calculations on PD1-ChlD1 is a linear combination of local excitations on ChlD1 and PD1, whereas the first CT state (root 9, PD1– ChlD1+) 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 ChlD1. In addition, the lowest CT state (root 5) changes in directionality (PD1+ ChlD1–) and is found much lower in energy, at 3.130 eV. Similar observations are made for the PD2–ChlD2 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 PD2), lowers the energy and alters the directionality of the first CT state (PD2+ ChlD2–, root 5, 3.232 eV, see also Table S7).

In the case of the ChlD2-PheoD2 pair (see Figure S8), the gas-phase calculations indicate that the lowest energy state (1.926 eV) is a local excitation on ChlD2, whereas the lowest CT state (ChlD2+PheoD2–, 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 ChlD2) and a significant stabilization of the CT state (ChlD2+PheoD2–, root 3) at 2.092 eV.

The most profound demonstration of protein matrix control is observed in the case of the ChlD1–PheoD1 pair. The nature of the lowest excited state at 1.905 eV obtained from the gas-phase calculations on the ChlD1–PheoD1 is primarily a local excitation on the ChlD1. The first CT state (root 9, ChlD1+PheoD1–) 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 ChlD1+PheoD1– charge transfer character (Figure ). 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 ChlD1 to PheoD1. 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 ChlD1+PheoD1– charge-transfer character (vide infra).

The above results are in line with observations from Stark spectroscopy,[36,37] where the ChlD1+PheoD1– CT state was found to be mixed with LE on ChlD1. This mixing was proposed to be crucial for initiating charge separation in the RC. A similar observation was made by Valkunas and co-workers,[66] where a good fit of the Stark spectrum (using complex time-dependent Redfield theory) was obtained with inclusion of the ChlD1+PheoD1– CT state, which was found to be more important in reproducing the Stark spectrum than the CT state originating from PD1–PD2+.

Calculation of the excited states of the ChlD1–PheoD1 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 (S9, 3.728 eV) is only marginally stabilized by CPCM (S8, 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 ChlD1 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 ChlD1–PheoD1 pair.

Excited States of the PD1–PD2–ChlD1 Trimer

The results presented above firmly support a CT excited state of the type ChlD1+PheoD1– 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 PD1 or the PD1–PD2 pair with charge transfer to ChlD1. 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 (PD1, PD2, and ChlD1). The TD-DFT QM/MM results presented and analyzed in terms of NTO compositions in Table S17 show that the lowest excited state of the trimer (S1 at 1.836 eV) is fully localized on ChlD1, which is consistent with the attribution of the lowest site energy of the RC to this pigment. The second and third (S2 at 1.880 eV and S3 at 1.910 eV) are localized excitations on PD1 and PD2. Charge transfer states begin to appear above 3 eV. The lowest is a CT state within the PD1–PD2 pair (PD1+PD2–, S7 at 3.050 eV), while the first CT state with PD1+ChlD1– CT character is S8 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 ChlD1, and strongly disfavor the hypothetical participation of an anionic ChlD1 species in native PSII charge separation.

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 excitation energies remains the same, i.e., ChlD1 has the lowest site energy (Figure ). Focusing on chromophore pairs, we find that the lowest-lying excited state of the RC remains on the ChlD1–PheoD1 pair and retains its CT character (ChlD1+PheoD1–) 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 depicts difference densities for the lowest excitation of the ChlD1–PheoD1 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.

Figure 9

Effect of protein conformational flexibility on the first excited state (in eV) of the individual chromophores (a) and on the two lowest excited states of the PD1–PD2 and ChlD1–PheoD1 pairs of chromophores (b) from ωB97X-D3(BJ) TD-DFT QM/MM calculations on QM/MM geometries optimized individually for each snapshot. Snapshot 1 represents the “crystal-like-conformation” of the protein, whereas snapshots 2–5 were derived from the production molecular dynamics simulations with regular intervals of 4 ns.

Figure 10

Difference densities describing the lowest excitation of the ChlD1–PheoD1 pair with exclusion (top row) and inclusion (bottom row) of protein electrostatics. With explicit consideration of the protein matrix the lowest excited state of the RC has dominant ChlD1+PheoD1– charge transfer character in the “crystal-like-conformation” of snapshot 1 and exclusive charge transfer character in the production MD snapshots 2–5.

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 PD1–PD2 pair (Figures and S9), whereas there is a significant impact on the ChlD1–PheoD1 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 S1 states with dominant (snapshot 1) to pure (all snapshots along the MD) ChlD1+PheoD1– 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 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.

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 far-red 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 ChlD1 cation radical is trapped at cryogenic conditions (5 K) and that Tyr161 (YZ) is the preferred electron donor in this case over the Cyt-b559/ChlZ/CarD2 pathway (see Figure ) in far-red light at 5 K. This is fully consistent with our results, which show the “dark” low-lying CT states of the ChlD1–PheoD1 pair to be in the red and far-red region. This implies that after charge separation within this pair takes place, the cation radical would initially be formed on ChlD1. The hole would subsequently migrate to PD1 under physiological conditions. Our results are also in very good agreement with the suggestion that YZ becomes the preferred electron donor under far-red light at 5 K. The migration of the ChlD1 hole to the special pair would require a reorganization of the protein environment that is probably inhibited at cryogenic temperature, resulting in YZ becoming an electron donor to ChlD1. Further EPR experiments[39] showed that charge separation upon light excitation is wavelength-dependent, leading to the hypothesis that PD1 is excited with visible light (532 nm), whereas ChlD1 is excited with far-red light.[39] Our results on the PD1–PD2–ChlD1 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 Nürnberg et al.[49] assigned Chl f to be the ChlD1 position, while recent work by Judd et al.[191] demonstrated the likelihood of PD2 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 ChlD1–PheoD1 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 O2 molecule, there is reduced likelihood that the OEC can advance regularly under these conditions, i.e., that subsequent far-red charge separations can occur on-time to advance the catalytic cycle outcompeting recombination. The fact that low O2 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 long-range 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 ChlD1 and the lowest excited state within the reaction center is a state of charge transfer character localized at the ChlD1–PheoD1 pair. Therefore, the present results support assigning this pair as the site of initial charge separation in PSII. The central PD1–PD2 chlorophylls do not present low-energy charge transfer excitations internally, while the possibility of charge transfer excitation from this pair to ChlD1 is excluded. Finally, protein dynamics have only weak influence on the localization of low-energy excitations, but enable charge transfer excitations within the ChlD1–PheoD1 pair to occur with far-red light.