Asymmetry in Charge Transfer Pathways Caused by Pigment–Protein Interactions in the Photosystem II Reaction Center Complex

: This article discusses the photoinduced charge transfer (CT) kinetics within the reaction center complex of photosystem II (PSII RC). The PSII RC exhibits a structural symmetry in its arrangement of pigments forming two prominent branches, D1 and D2. Despite this symmetry, the CT has been observed to occur exclusively in the D1 branch. The mechanism to realize such functional asymmetry is yet to be understood. To approach this matter, we applied the theoretical tight-binding model of pigment excitations and simulated CT dynamics based upon the framework of an open quantum system. This simulation used a recently developed method of computation based on the quasi-adiabatic propagator path integral. A quantum CT state is found to be dynamically active when its site energy is resonant with the exciton energies of the PSII RC, regardless of the excitonic landscape we utilized. Through our investigation, it was found that the relative displacement between the local molecular energy levels of pigments can play a crucial role in realizing this resonance and therefore greatly affects the CT asymmetry in the PSII RC. Using this mechanism phenomenologically, we demonstrate that a near 100-to-1 ratio of reduction between the pheophytins in the D1 and D2 branches can be realized at both 77 K and 300 K. Our results indicate that the chlorophyll Chl D1 is the most active precursor of the primary charge separation in the D1 branch and that the reduction of the pheophytins can occur within pico-seconds. Additionally, a broad resonance of the active CT state implies that a large static disorder observed in the CT state originates in the ﬂuctuations of the relative displacements between the local molecular energy levels of the pigments in the PSII RC.


Introduction
Photosystem II (PSII) is the only biological system that has the unique ability to oxidize H 2 O to O 2 [1]. The oxygen produced in the photosystem II reaction center (PSII RC) provides the source of oxygen in Earth's atmosphere providing the essential basis of life.
The PSII RC is comprised of D1 and D2 proteins, core antenna proteins CP43 and CP47, and several small subunits [1,2]. Similar to reaction centers found in purple bacteria, the D1 and D2 proteins comprise the core of the PSII RC along with a number of noncovalently associated cofactors, including pheophytins, two quinones, and an iron ion (Fe 2 + ), as shown in Figure 1 [3]. Located near the PSII RC, the oxygen-evolving complex (OEC), made up of a Mn 4 Ca cluster, is responsible for water oxidization [4,5] which makes the PSII RC unique.
The entire structure of the PSII RC includes two closely related proteins which form one large unit dimeric in structure, resulting in a symmetrical complex. [2,6]. This symmetry in the PSII RC can  [3]. The oxygen-evolving complex (OEC) and the bicarbonate ligand (BCT) are also shown. Molecular mechanisms containing functional asymmetry have been investigated, but the cause and function of this asymmetry are still under debate [1,11]. Spectral analysis of the charge transfer kinetics in the PSII RC has been examined thoroughly and has shown that the electron donor for primary charge separation in PSII RC involves Chl D1 [12] as well as the central pair, P D1 and P D2 [13][14][15][16]. Thus even within the D1 branch, multiple pathways for charge transfer are observed to be active. This also implies that the energetic landscape of the D1 is different from that of D2 which possibly serves as the cause of the asymmetry in the photochemical pathway [17]. Recent theoretical studies of molecular dynamics (MD) using the quantum mechanics/molecular mechanics (QM/MM) approach support the concept of an energy landscape that favors active charge transfer in the D1 branch [11,18]. However, to what extent this energy landscape could affect the CT dynamics in the complex requires further investigation.
In this work, a theoretical model of a mechanism by which the charge transfer kinetics in the PSII RC are highly asymmetric is presented. The pigments in the PSII RC non-covalently interact with local protein residues, causing variations in molecular excitation energies (site energies) among the pigments [11,18]. We postulate that such interactions cause variances in the relative displacements among the local molecular energy levels of the pigments as well as their excitation energies. To investigate effect of the molecular energy levels on the charge transfer phenomena, we adopt the tight-binding model of the PSII RC [19][20][21] as a starting point and then incorporate the energy levels as parameters to the model. We then simulate the dynamics of photo-excited charge transfer for various parameter settings using a non-Markovian and non-perturbative method of computation [22]. We find that charge transfer states are active only when their excitation energies are resonant with those of delocalized excitons of the PSII RC regardless of the particular excitonic landscape. Thus we show that the activity of charge transfer states is controlled by the shifting of the local molecular energy levels and the strong asymmetry of charge separation between the D1 and D2 branches which can be explained quantitatively in a phenomenological manner.

Frenkel Excitation and Charge Transfer States of Pigment-Protein Complex
We represent the reaction center of PSII as a pigment-protein complex (PPC), a molecular aggregate of M pigments being held in a protein scaffolding. Molecular excitations of the pigments in the PPC are described in terms of electrons and holes residing in the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) within each pigment [19]. The ground state of the PPC, |0 , is defined as the state in which all of the HOMOs are occupied and all of the LUMOs unoccupied (see Figure 2a). The excited state of the PPC, |m + n − is defined to have one hole in the HOMO at pigment m and one electron in the LUMO at pigment n (see Figure 2b). In what follows, we refer to the excited state for m = n as the charge transfer (CT) state and that for m = n as the Frenkel exciton (FE) state which is denoted by |m * = |m + m − (see Figure 2c). The total charge is conserved under these excitations so the PPC as a whole remains neutral. The pigments in the PPC are affected by vibrational fluctuations originating in the protein scaffolding and surrounding solvent. To incorporate this phenomenon, the Hamiltonian of the PPC is divided into three parts asĤ =Ĥ S +Ĥ B +Ĥ SB , whereĤ S is the system Hamiltonian describing the excitation dynamics of the pigment sites,Ĥ B is the bath Hamiltonian describing the environmental vibrations, andĤ SB is the system-bath Hamiltonian that determines the way that the pigments are influenced by the vibrations. For the system Hamiltonian, we employ the tight-binding model of PPC [19], whose matrix elements are specified by where t nm ) are the tunneling matrix elements [23], U mn is the (renormalized) Coulomb interaction potential between the hole in the HOMO at pigment m and the electron in the LUMO at pigment n, and V mn (= V nm ) is the Förster resonance coupling between the FEs at pigments m and n, satisfying V mm = 0. For later convenience, we define the site energy ε mn by the diagonal elements (k = m and = n) ofĤ s , This represents the amount of energy it takes to excite the ground state |0 to |m + n − . The matrix element between two FE states is obtained from Equation (1) as m * |Ĥ s |n * = ε mm δ mn + V mn , indicating that the coupling between the FE states are determined only by V mn . In contrast, coupling between CT states and between FE and CT states are determined only by the off-diagonal elements of the tunneling matrix, t (±) mn of m = n. These are associated with the overlap integrals of the electronic wavefunctions related to the molecular orbitals which are well approximated by the exponential function of distance [19,23]. We therefore model the tunneling matrix by where A (±) mn is the coefficient representing the overall strength of coupling satisfying A (±) mm = 0, R m is the representative position of pigment m, and a mn is a length scale that generally depends on the pigment pair [19].
The bath of vibrations described byĤ B originates in the protein structure and surrounding solvent [23]. The influence of these vibrations on the system described byĤ SB gives rise to a renormalization of the interaction potentials U mn in Equation (1). Such renormalization takes the form of where U (C) mn is the (bare) Coulomb potential of the hole-electron interactions andhλ mn is the reorganization energy of the excited state |m + n − (including |m * for m = n). The Coulomb interaction potential we use here is where is the dielectric constant, and η is a cutoff parameter yielding the exciton binding energy for n = m that we denote by U B = U , is by definition the amount of energy it takes to excite an electron from the HOMO to the LUMO state [1]. This value is identical to the amount of energy required to create the FE at pigment m at the ground state. Thus, the binding energy of FE, U B , is included in the HOMO-LUMO gap. Based on this, we assign t (+) mm the amount of energy required to create a hole in the HOMO of pigment m by removing an electron from the ground state to the vacuum level, that is t Here, the biding energy is subtracted from the energy level of the LUMO level because an electron, instead of a hole, resides in the HOMO.
Incorporating our argument above, we introduce an alternative expression for the site energy.
Substituting t where ε 0 mn = ε nn −hλ nn + U mn − U B which we refer to as the zero-shift site energy of |m + n − . Since ∆E 1 = 0 holds, the site energies of Equation (6)

The Model of the PSII Reaction Center
The reaction center of the PSII contains four chlorophylls (P D1 , P D2 , Chl D1 and Chl D2 ), two pheophytins (Pheo D1 and Pheo D2 ), and two plastoquinones (Q A and Q B ). The timescale of the kinetics of the charge transfer between the pigments is on the order of a few pico-seconds (ps) whereas the timescale of Q A reduction is on the order of 100 ps [24]. For this reason we exclude the plastoquinones from our model, focusing on the first few pico-seconds (ps) after the photo-excitation of the complex. In what follows, for the sake of notational simplicity, we also refer to each of the six pigments by numbers (1 through 6) assigned as shown in Figure 4.  For the FE-FE sector of Equation (1), we adopt the site energies, ε mm for m = 1, 2, · · · , 6, from Gelzinis et al. [21] and the FE couplings V mn from Shibata et al. [25]. The relevant values are illustrated in Figure 4b and are also summarized in Table A1. For the CT-CT and FE-CT sectors of Equation (1), we must obtain the zero-shift site energies ε 0 mn and the off-diagonal elements of tunneling matrix, t (±) mn . The parameters specifying ε 0 mn are ε mm , , η, and λ mn , among which ε mm are the previously determinded site energies of the FEs. For the protein environment of the PSII RC, we use = 1.5 as suggested by Müh and Renger [11]. For η, we examined the zero-phonon lines given by Novoderezhkin et al. [15] and found that η = 7.4 Å is an optimum value. This η yields the binding energy of FE, U B = U (C) mm = −10463 cm −1 −1.30 eV, for all of the pigment sites. The positions R m have been adopted from Gelzinis et al. [21]. These values are summarized in Table A2. The reorganization energy for the FE states are obtained once the spectral density of the pigment-protein interaction is determined. As discussed in more depth later on in the Methods section, we sethλ mm = 50 cm −1 for all of the FE states and (m = n)hλ mn = 3 λ mm = 150 cm −1 for all of the CT states. The zero-shift site energies are then evaluated for each of the 36 states (6 FE states and 30 CT states), the results from which are shown in Figure 5a. Note that the site energies of the FEs are identical to their zero-shift site energies, ε mm = ε 0 mm , due to Equation (6). To specify the CT-CT and FE-CT couplings, we further require the inclusion of A (±) mn and a mn in the tunneling matrix elements of Equation (3). For these parameters, we aim to adjust the coupling between P * D2 and P + D1 P − D2 to be −75 cm −1 as given by Novoderezhkin et al. [26]. This can be achieved by taking A (±) mn = −415 cm −1 assuming a mn = 4.5 Å [19]. These values are used for all m and n of m = n for simplicity. The results of the coupling strength, | m + n − |Ĥ S |k + − | for |m + n − = |k + − , are shown in Figure 5b.   We now investigate how the population dynamics of the excited states depend on the level shifts of the pigments in the PSII RC. The populations of excited states are represented as the diagonal elements of the reduced density matrix, m + n − |ρ(t)|m + n − . In order to see correlations among the FE and CT populations, we quantify the transfer efficiency by the time average of population for each state, where τ is the time over which the average is computed. It has been observed in past experimental studies that the PSII RC strongly prefers the reduction to occur at Pheo D1 , the pheophytin in the D1 branch. To characterize this asymmetry, we define the efficiency of electron transfer (ET) to Pheo D1 by ET D1 = ∑ 6 m =5 p m5 , and likewise the efficiency to Pheo D2 by ET D2 = ∑ 6 m =6 p m6 . Their difference, is a measure of the asymmetry. ∆ET > 0 stipulates that ET is more likely to occur at the D1 branch, as observed experimentally, while ∆ET < 0 specifies that ET is expected to occur at the D2 branch.

Level Shifts in the D1 Branch
As the starting point of the analysis of charge transfer phenomena in the PSII RC, we first set all of the level shifts to zero so that mn = 0 mn . Then, a simulation of the FE and CT states after the photo-excitation of the PSII RC by a δ pulse was run to visualize the population dynamics within the system. The result for the first ps of the simulation is shown in Figure 6. Although each of the 36 states (see Figure 5) are included in the computation, the populations are dominated by FE states while the contribution from CT states is limited. Within the FE states, a major transition is observed to occur during the first 0.4 ps. After this event, |3 * and |6 * steadily increase while the others decrease. This is mainly because |3 * and |6 * have the lowest site energies (see Figure 4b) within the D1 and D2 branches respectively and therefore act as sink sites of the exciton in the PSII RC. Similarly, the two CT populations of |1 + 2 − and |2 + 1 − seen in Figure 6 have the first and the second lowest zero-shift site energies respectively amongst all of the CT states (see Figure 5a). Considering this similarity, we therefore expect the CT states to become more excited as their site energies are lowered by manipulating the level shifts. Population 1 * 2 * 3 * 4 * 5 * 6 * 1 + 2 2 + 1 Figure 6. Population dynamics of FE and CT states at T = 77 K. All of the levels shifts are set to zero, ∆E m = 0 for m = 2, · · · , 6. All of the 36 states in Figure 5 are used for computation, but those of populations less than 0.01 are not displayed.
First, we evaluate which CT states can be excited if the level shifts are changed in the D1 branch. We begin by observing how the level shift at Pheo D1 (m = 5), that is ∆E 5 , affects the transfer efficiency by maintaining all other level shifts at zero, ∆E m =5 = 0. Simulations were run to obtain a range of ∆E 5 values to secure transfer efficiencies for τ = 1 (in Equation (7)). The results of the simulations are shown in Figure 7a. A remarkable feature can be seen in the data collected of the transfer efficiencies in the form of a large peak of p 35 (dashed green) at ∆E 5 = −0.422 eV (= −3400 cm −1 ) accompanied by large fractional decreases in p 33 (solid green) and in p 55 (solid blue). This peak indicates that the CT state |3 + 5 − is being activated by the two FE states, |3 * and |5 * . This is interpreted as a formation of the radical pair Chl + D1 Pheo − D1 by the photo-induced charge transfer process occurring between the excited molecules Chl * D1 and Pheo * D1 . In the upper panel of Figure 7a, there is a prominent positive peak in ∆ET associated with the formation of Chl + D1 Pheo − D1 , representing an efficient electron transfer to Pheo D1 . In addition to the large peak in p 35 , there is a very small peak of p 15 (dashed red) at ∆E 5 = −0.707 eV (= −5700 cm −1 ) corresponding to a slight ET increase in the D1 branch. The peaks in the transfer efficiencies can be interpreted in terms of resonance and coupling strength among the excited states. The site energies of the states for the range of ∆E 5 values are shown in Figure 7b. Resonance between the states can be induced when their site energies are close together, illustrated by the crossing points of the site energies in Figure 7b. The five CT states |m + 5 − for m = 5 whose site energies depend on ∆E 5 as ε m5 = ε 0 m5 + ∆E 5 , have site energies which are seen to cross with all of the FE states within the parameter window. Inspecting Figure 7a, however, only two of them, (|3 + 5 − and |1 + 5 − ), are found to be active in the actual dynamics. Such selection of CT states can be understood in terms of coupling strength between the excited states. The CT states strongly coupled with |3 + 5 − are shown in the upper panel of Figure 7b. The FE states, |3 * and |5 * , are both strongly coupled with |3 + 5 − , thus active transfers between these states can occur around the crossing point, ∆E 5 −ε 0 35 = −0.395 eV. This is exactly what is observed in Figure 7a where ∆E 5 is lower than this value by the amount of ∼0.02 eV(∼150 cm −1 ). This shift of the resonance point is on the order of 100 cm −1 , indicting that it can be induced by pigment-pigment and pigment-protein interactions. The CT state |1 + 5 − is also strongly coupled with |3 + 5 − . However, because their site energies do not cross each other, the charge transfer will not occur between them. Likewise, a transfer between |2 + 5 − and |3 + 5 − does not occur because their site energies also do not cross each other. The small peak of p 15 in Figure 7a is another result of resonance between |1 + 5 − and the FE states, |1 * and |5 * , around ∆E 5 −ε 0 15 = −0.703 eV. However, the charge transfer here is not as active as that of |3 + 5 − because the populations of |1 * and |5 * are small (see Figure 7a). Additionally, the coupling strength between these FE states and |1 + 5 − is weak at 14 cm −1 (see Figure 5b), further restricting the CT activity. The rest of the linearly changing states in Figure 7b, which are |2 + 5 − , |4 + 5 − , and |6 + 5 − , remain inactive because none are coupled strongly enough with the FE states, as can be noted in Figure 5b. There is an additional small peak of p 35 at ∆E 5 = −0.285 eV(= −2300 cm −1 ) in Figure 7a whose origin is yet to be identified.
Next, the effects of the level shift at Chl D1 (m = 3) on the transfer efficiencies are examined. Simulations were run for a range of ∆E 3 keeping all other level shifts unchanged at zero, ∆E m =3 = 0. The results from these simulations are shown in Figure S1a of the Supplementary Materials. Remarkably, three large peaks of the CT states are seen in the transfer efficiency, p 13 (dash-dotted red) at ∆E 3 = −0.360 eV(= −2900 cm −1 ), of p 53 (dash-dotted blue) at ∆E 3 = −0.409 eV (= −3300 cm −1 ), and of p 23 (dash-dotted yellow) at ∆E 3 = −0.533 eV (= −4300 cm −1 ). As shown in Figure S1b of the Supplementary Materials, these peaks correspond to the resonance of the three CT states, |1 + 3 − , |5 + 3 − , and |2 + 3 − , with FE states at their crossing points. Thus three radical pairs P + D1 Chl − D1 , Pheo + D1 Chl − D1 , and P + D2 Chl − D1 are formed by this resonance. The site energies of |4 + 3 − and |6 + 3 − are also crossing with those of FE states (shown in Figure S1b), but they remain inactive because they are coupled only weakly with FE states, as can be seen in Figure 5b. Now it must be considered how the population of the radical pair P + D1 Pheo − D1 (|1 + 5 − ) can dominate the terminal sate of the PSII RC. We do this by seeking charge transfer pathways that maximize the transfer efficiency p 15 . As can be seen in Figure 5b, the FE states do not couple strongly with |1 + 5 − , so |1 + 5 − needs to be excited by other CT states. We refer to such states as the primary CT states. Inspecting Figure 5b, |1 + 5 − is strongly coupled with |1 + 3 − , |2 + 5 − , |3 + 5 − , and |4 + 5 − . Meanwhile, |2 + 5 − and |4 + 5 − are not as active because they are only slightly coupled with the FE states. Since no coupling exists between |1 + 3 − and |3 + 5 − , they independently couple with |1 + 5 − . Thus, we have narrowed down the possibilities to the following two pathways: FEs → |1 + 3 − → |1 + 5 − and FEs → |3 + 5 − → |1 + 5 − . The first pathway depends on the activity of |1 + 3 − , which can be excited efficiently by the FE states around the resonance peak at ∆E 3 = −0.360 eV −ε 0 13 as seen in Figure S1a  The second possible pathway to excite P + D1 Pheo − D1 is FEs → |3 + 5 − → |1 + 5 − . This pathway is expected to be active when these CT states are mutually in resonance by ε 35 ε 15 . From Figure 7a, we have gathered that the transfer efficiency of |3 + 5 − is maximized when ε 35 = ε 0 35 − 0.422 eV. Using Equation (6), this implies that the condition ∆E 3 = ∆E 5 + 0.422 eV will induce the highest activity of |3 + 5 − . To examine the activity of |1 + 5 − in the pathway, simulations were run for a range of ∆E 5 values under this condition. The results from these simulations are shown in Figure 8b. The peak of p 15 (dashed red) is observed at ∆E 5 = −0.744 eV (= −6000 cm −1 ) accompanied by the largest decrease in p 35 , indicating that the pathway is most active at this point. Interestingly, we have obtained the same ∆E 5 value as that of the first pathway. This leads to ∆E 3 = −0.322 eV (=−2600 cm −1 ) for the activation of the second pathway, which only deviates 0.038 eV (300 cm −1 ) from that of the first pathway.
The radical pair, P + D1 Pheo − D1 (|1 + 5 − ), can thus be excited efficiently from the FE states mediated by the two primary CT states, P + D1 Chl − D1 (|1 + 3 − ) and Chl + D1 Pheo − D1 (|3 + 5 − ). To visualize how these efficient pathways impact the population dynamics, simulations were run for an extended time up to 10 ps at T = 77 K. The result for the first pathway, FEs → |1 + 3 − → |1 + 5 − , is shown in Figure 9a. A very quick excitation of the radical pair is observed as its population (dashed red) reaches the greatest value within the first 2 ps. Paying close attention on the first 0.4 ps, the population dynamics of the FE states shown in the inset of Figure 9a are observed to be similar to those in Figure 6, except that the population of |3 * in Figure 9a (solid green) does not increase monotonically but begins to decrease around 0.1 ps. This difference indicates that |3 * is the primary donor of the excitation energy for the three CT states, |1 + 5 − , |1 + 3 − , and |3 + 5 − . From 0.4 ps to 2 ps, the population of |1 + 5 − continues to grow until it reaches its maximum value of 0.38 at 1.9 ps, whereas the populations of |3 * and |1 + 3 − continue to decrease exponentially. Hence the radical pair P + D1 Pheo − D1 causes the dominant population. After the first 2 ps, the excited CT states remain active by keeping their populations mostly steady but with a slow decrease in |1 + 5 − over time. Note that the population of |6 * (solid cyan) appears to behave differently from other states because it belongs to the D2 branch. During the course of population change, ∆ET is always positive and increasing until t = 10 ps, indicating that the election continues to be transferred to the D1 branch through the interplay between |1 + 5 − and other states. As a result of this, ET D1 and ET D2 reaches 0.664 and 0.003, respectively, at 10 ps. This produces a ∆ET value of 0.661 in Figure 9a. Viewing the asymmetry by the ratio of these ET populations, that is (ET D1 − ET D2 )/(ET D1 + ET D2 ) = 99.1%, it is suggested that ET in the PSII RC is extremely biased toward the D1 branch by the charge transfer pathway.
At the higher temperature of 300 K, the populations shown in Figure 9b appear to behave similarly to those at 77 K. The major transition at both temperatures happens within the first 3 ps and the population of |1 + 5 − is the largest. However, the population of |1 + 5 − is much lower with a maximum value of only 0.16 at 2.7 ps while other CT states appear to be more active at the higher temperature. In particular, the populations of |1 + 2 − (dotted red) and |2 + 1 − (dotted yellow) are more pronounced compared to those at 77 K, possibly due to thermal excitations induced in the FE and CT states via the system-bath coupling. In regard to the ET dynamics, the overall value of ∆ET is also lower but is still positive and increasing over time. At 10 ps, ET D1 and ET D1 are 0.411 and 0.012 respectively, making the asymmetry ratio of (ET D1 − ET D2 )/(ET D1 + ET D2 ) = 94%. This large value indicates ET in the PSII RC is still substantially biased toward the D1 branch even at 300 K by the active pathway.
The result of our simulation for the second pathway, FEs → |3 + 5 − → |1 + 5 − , is shown in Figure S2 of Supplementary Material. The dynamics of |1 + 5 − are similar to those of the first pathway, but its initial ascent progresses slightly slower and overall population is about 3/4 as much of that of the first pathway. The population reaches its maximum value of 0.29 at 2.2 ps, and then decays slowly over time. Despite the lower population in the radical pair, the ET and its asymmetry are comparable to that of the first pathway. At 10 ps, ET D1 and ET D2 reaches 0.611 and 0.003 at 77 K, respectively, which amounts to the asymmetry ratio of (ET D1 − ET D2 )/(ET D1 + ET D2 ) = 98.9%. At 300 K, ET D1 and ET D2 reach 0.476 and 0.011 respectively which produces a asymmetry ratio of 95.3%. Thus, the ET is highly biased toward the D1 branch at both 77 K and 300 K.

Level Shifts in the D2 Branch
We now turn our attention to the other branch of the PSII RC, the D2 branch, which has been experimentally observed to be mostly inactive in ET dynamics. Here, we examine to what extent the D2 branch could be active in the PSII RC. In particular, we seek a pathway to activate the radical pair P + D2 Pheo − D2 , the counter part of P + D1 Pheo + D1 in the D2 branch. As we have done for the D1 branch, we change the level shifts of the pigments in the D2 branch to maximize the transfer efficiency p 26 . For simplicity, we keep ∆E m = 0 for m = 2, 3, 5, and adjust ∆E 4 and ∆E 6 to form a resonance that activates |2 + 6 − . The results of the simulations for a range of ∆E 4 values are shown in Figure S3a of the Supplementary Materials. By analogy with the FEs → |1 + 3 − → |1 + 5 − pathway of the D1 branch, we consider the FEs → |2 + 4 − → |2 + 6 − pathway here. The resonance peak of p 24 (dash-dotted yellow) is found at ∆E 4 = −0.384 eV (=−3100 cm −1 ), and thus we use this to search for the ∆E 6 value that maximizes p 26 . The results of the simulations are shown in Figure S3b of Supplementary Material. We locate the peak of p 26 (dashed yellow) at ∆E 6 = −0.694 eV (= −5600 cm −1 ).
With these values of level shifts, we see how the populations of the radical pair P + D2 Pheo − D2 are formed through dynamics. The first 10 ps of the simulation run at 77 K is shown in Figure 10.
Comparing this with Figure 9a, the behavior of active CT states are observed to be similar to their counter parts in the D1 branch, which are |2 + 6 − ⇔ |1 + 5 − , |4 + 6 − ⇔ |3 + 5 − , and |2 + 4 − ⇔ |1 + 3 − . In particular, the major transitions that occur during the first 3 ps of the simulation and the population is dominated by the radical pair, |2 + 6 − . However, the maximum population of |2 + 6 − is recorded at 0.29 in Figure 10, so it is lower than its D1 counterpart, |1 + 5 − in Figure 9a, by a factor of 0.76. This ratio roughly holds for their subsequent populations after 3 ps. This difference in the populations of the radical pairs can be caused by the difference in molecular species of their primary donor sites. The donor of |2 + 6 − is the FE of a pheophytin |6 * , whereas the donor of |1 + 5 − is the FE of a chlorophyll |3 * . As can be seen in Figure 6, the initial populations of |6 * and |3 * in our simulations are 0.115 and 0.193 respectively which implies that the population of |6 * is lower by a factor of 0.6. This is partly the consequence of the transition dipole moment of pheophytin and that of chlorophyll, which are 3.4 debye and 4.4 debye, respectively, making the ratio of oscillator strengths roughly (3.4/4.4) 2 = 0.6. Thus, the activity of the radical pairs is largely affected by the intrinsic photo-activity of their donor FE states. Although the population of P + D2 Pheo − D2 is low in Figure 10, the asymmetry in ET is strongly biased toward the D2 branch. This is quantified by the ET populations in the D1 and D2 branches, ET D1 and ET D2 which are 0.03 and 0.49 at 10 ps, respectively, yielding the asymmetry ratio of (ET D1 − ET D2 )/(ET D1 + ET D2 ) = −99%, indicating the electron is almost certainly found in the D2 branch. This is the complete opposite to what is observed in experiments.

Discussion
Our results show that the activities of the CT states in the PSII RC depend highly on the level shifts representing the relative displacements in the molecular energy levels of the pigments. The level shifts control whether a CT state is resonant or not with the FE and other CT states. When a CT state is resonant and coupled with the photoactive FE states, it becomes the primary CT state. If the CT state is resonant and coupled instead with the primary CT state, it becomes the secondary CT state.
Our results also show that the energy landscape of the FE states is not the determining factor of the CT asymmetry in the PSII RC. It was demonstrated that the CT in the D2 branch could be activated by a specific arrangement of level shifts, even with the FE landscape favoring the D1 branch to be more excitonic than the D2 branch [11,18]. Thus our results suggest that in the naturally occurring PSII RC, the D1 branch is exclusively active as a consequence of the level shifts arranged in a way so that the CT states in the D1 branch are resonant with the FE states while those in the D2 branch are not.
The CT states thus activated in the D1 branch are seen to be robust around their corresponding resonances. Our results show that the widths of resonant peaks in the transfer efficiencies are around 0.05 eV ( 400 cm −1 ) (see Figure 8). This is in agreement with experimental results which show that the static disorders of the CT states are much larger than those of the FE states which are typically around 50 cm −1 [20,21]. Based on this, we speculate that the site energy disorders of CT states are induced by fluctuations of the level shifts (difference in HOMO levels between pigments, see Figure 3c) in the PSII RC. Meanwhile, the disorders of FE states are due to fluctuations of HOMO-LUMO gaps within individual pigments [11,18]. Therefore the hypothesis can be tested by examining the correlations between the disorders of the CT states based on Equation (6).
Concerning the timescale, our simulations suggest that the formation of the primary CT states, Chl + D1 Pheo − D1 and P + D1 Chl − D1 , occurs during the first 0.4 ps after the initial photo-excitation, followed by the formation of the secondary CT state, P + D1 Pheo − D1 , within the first 2 ps at 77 K and 3 ps at 300 K (see Figure 9). On the formation of the primary CT states, our result is largely consistent with a sub-pico-second timescale for the reduction rate of Pheo D1 obtained experimentally by Groot et al. [12] and then further quantified to 0.3 ps by Raszewski and Renger [9,27]. This is also similar to the timescale of the FE dynamics (see Figure 6) as the strength of coupling between the FEs and the primary CT states (∼40 cm −1 ) is similar to the FE-FE couplings in the D1 branch. However, several experiments using two-dimensional electronic spectroscopy have reported pico-second timescales; 1-3 ps at 77 K by Myers et al. [28] and 1.5 ps at room temperature by Duan et al. [29]. Furthermore, on the formation of the secondary CT state, Groot et al. [12] reported a time scale of ∼ 6 ps for the oxidation of P D1 at room temperature, which is about two times larger than our result. These experimental results indicate that the FE-CT and CT-CT couplings are likely to be weaker than those we used here, so the off-diagonal elements of the tunneling matrix (Equation (3) for m = n) needs to be refined in future work.
Another prospect to examine is an additional CT pathway involving the central pair, P + D1 P − D2 and P + D2 P − D1 , which has been proposed by van Grondelle and coworkers [13][14][15][16]. Although we did not explore this possibility here, such a pathway can be realized in our model by carefully adjusting the level shift ∆E 2 while keeping all other CT states in the D2 branch off resonant with the FE and CT states.
In this work, we focused on how the level shifts can contribute to producing the CT asymmetry seen in the PSII RC. For this naturally occurring complex, we estimate that the HOMO levels in the D1 branch are progressively lowered from P D1 to Chl D1 and Chl D1 to Pheo D1 roughly by the interval of ∼0.35 eV, which amounts to ∼19% of the FE site energy of each pigment. We theroize that such shifts may be realized by the pigment-protein interactions. However, it is still a question as to which particular interactions at each of the pigment sites can cause such ideal displacements between the energy levels. Additionally, we did not consider the site dependence of the off-diagonal elements of the tunneling matrix which can also contribute to the CT asymmetry. We anticipate that further theoretical investigations using ab initio approach, such as MD and QM/MM, can reveal the precise nature of the local molecular energy levels in the PSII RC.

Tight-Biding Model of Molecular Excitations and Pigment-Protein Interactions
In the tight-binding model, the excited states |m + n − of PPC are constructed by using two sets of fermionic operators,ĉ m andd m for m = 1, 2 · · · , M, satisfying the anti-commutation relations for arbitrary m and n, ĉ m ,ĉ n = d m ,d n = 0, and ĉ m ,ĉ † n = d m ,d † n = δ mn . Asĉ m andd m belong to different degrees of freedom, the commutation relations hold as ĉ m , wherehΛ 1 = 15 cm −1 ,hΛ 2 = 35 cm −1 ,hΩ 1 = 20 cm −1 , andhΩ 2 = 80 cm −1 . This yields the reorganization energy of FE stateshλ m = 50 cm −1 . For the scaling factors in Equation (12), we simplify them by s , yielding the reorganization energy of CT sateshλ mn = 3hλ m = 150 cm −1 .

Photo Excitation of the PPC
We consider initial state of PPC prepared by photo-excitation. The transition dipole operator where µ m is the transition dipole moment of pigment m. This operator acts on the ground state to yield the transition dipole, where |E a is the ath eigenstate ofĤ S , satisfyingĤ S |E a = E a |E a .

Computing Quantum Dynamics with the Scalable QuAPI Method
The quantum dynamics of our model is described with the theory of open quantum systems which uses the reduced density operatorρ(t) [23,32]. This contains all of the information about the system under the influence of a thermal bath of temperature T. For numerical computations ofρ(t), we employ the scalable QuAPI method (S-QuAPI) [22] which is one of the recent updates of the quasi-adiabatic propagator path integral (QuAPI) scheme [33][34][35][36] for enhanced scalability and memory efficiency to deal with large quantum systems such as our model of the PSII RC. The method has been designed to be most effective for the modern architecture of massively-parallel platforms. For our simulations, we use a high-performance computer cluster utilizing 12 units of NVIDIA Tesla K80 GPU.
S-QuAPI has three parameters to control its accuracy of dynamics, which are the time slice of path integral ∆t, the number of time steps for preserving the memory effect ∆k max , and the threshold ϑ for dropping insignificant propagators from computations. The method is both non-Markovian and non-perturbative converging to the exact result at the limit of ∆t → 0, ∆k max → ∞, and ϑ = 0. Thus an approximation is made by setting finite values for these parameters. For our simulations, we set ∆t = 25 fs and ∆k max = 3 for all computations and optimized the ϑ values in between 5.0 × 10 −6 and 1.0 × 10 −5 depending on time t for whichρ(t) is evaluated. In general, smaller ϑ is required for larger t to achieve better accuracy [22].

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:  Table A2. Data of pigment sites adapted from Gelzinis et al. [21]. The coordinates of geometric center of pigment, R = (R x , R y , R z ), are in units of Å, and the transition dipole moment of pigment, µ = (µ x , µ y , µ z ), are in units of debye.