Efficiency of energy funneling in the photosystem II supercomplex of higher plants

Abstract

The investigation of energy transfer properties in photosynthetic multi-protein networks gives insight into their underlying design principles. Here, we discuss the excitonic energy transfer mechanisms of the photosystem II (PS-II) C₂S₂M₂ supercomplex, which is the largest isolated functional unit of the photosynthetic apparatus of higher plants. Despite the lack of a definite energy gradient in C₂S₂M₂, we show that the energy transfer is directed by relaxation to low energy states. C₂S₂M₂ is not organized to form pathways with strict energetically downhill transfer, which has direct consequences for the transfer efficiency, transfer pathways and transfer limiting steps. The exciton dynamics is sensitive to small changes in the energetic layout which, for instance, are induced by the reorganization of vibrational coordinates. In order to incorporate the reorganization process in our numerical simulations, we go beyond rate equations and use the hierarchically coupled equation of motion approach (HEOM). While transfer from the peripheral antenna to the proteins in proximity to the reaction center occurs on a faster time scale, the final step of the energy transfer to the RC core is rather slow, and thus the limiting step in the transfer chain. Our findings suggest that the structure of the PS-II supercomplex guarantees photoprotection rather than optimized efficiency.

---

1 Introduction

Photosynthesis, in which light is absorbed and converted into chemical energy, is the most important process in nature. In higher plants the light-harvesting machinery is assembled of C₂S₂M₂ supercomplexes and networks of LHCII pigment-proteins^1–4 located in the grana membrane. The C₂S₂M₂ supercomplex is formed by a dimeric photosystem II (PS-II) with loosely bound LHCII trimers and several minor complexes.⁵ Energy transfer to the reaction center (RC) core pigments of PS-II, in which the primary step of charge separation initializes an avalanche of photochemical reactions,^6–8 reaches remarkable efficiencies of up to 90%.⁹ However, it remains unclear of how such high efficiencies can be achieved in large and disordered systems. In contrast to the photosynthetic apparatus of green sulfur bacteria, in which fast transfer is guaranteed by efficient energy funneling,¹⁰ microscopically derived Hamiltonians do not predict a definite energy gradient among the individual proteins of the C₂S₂M₂ supercomplex.^11–15

Previous works describe the transfer kinetics with phenomenological models, and extract certain decay components, such as the migration time (average time that it takes for an excitation to reach the RC) and trapping time, by fitting to fluorescence decay lines.^16–19 Several rate limiting models are discussed in the literature.^17,18,20,21 Recent studies favor the so-called transfer-to-trap limited kinetic model^12,16–18 in which the transfer rate from the antenna complexes of PS-II to the RC is proposed to be the transfer-limiting step. However, different kinetic models can be fitted equally well to measured fluorescence decay curves,²² and structure-based models of energy transfer become necessary to shed light on the underlying transfer mechanisms. The first microscopic simulations of the exciton dynamics in the C₂S₂M₂ supercomplex show that the overall transfer is driven by a complex interplay of multiple rates rather than through a single transfer-limiting step.²²

In pigment–protein complexes, the directionality of energy transfer is driven by energy relaxation. Variations in the energy bands of the individual proteins in C₂S₂M₂ are not as distinct as in other photosynthetic systems. Nevertheless, the energy gradient in C₂S₂M₂ is not completely flat, and the pigments form a certain structure in the energetic layout. For example, CP43 and CP47 are lower in energy than the LHCII antenna complexes.^11,12,14 However, energy transfer in C₂S₂M₂ is not a cascade of downhill steps toward the reaction center. Actually, the pigments in the proximity of the RC core are the energetically lowest ones.²³ Therefore, the last transfer step to the trap needs to overcome an energy barrier which supports the proposed transfer-to-trap limited exciton dynamics in C₂S₂M₂. The transfer-limiting step to the RC core pigments, which is not anticipated in previous structure-based simulations,²² becomes more evident once we include the recently derived Hamiltonian of CP29.¹⁵ The latter is substituted in ref. 22 by an LHCII monomer. We show that the minor complex CP29 modifies the pathway of energy flow and yields a relaxation channel which drives energy from the peripheral antenna towards pigments closer to PS-II.

The transfer properties are sensitive to small structural modulations, which is an immediate consequence induced by the flat energy gradient. There are several mechanisms that induce changes in the energetic structure. These include static disorder, in which site energies are subjected to random fluctuations on much slower time scales than the exciton dynamics, or the reorganization process, in which vibrational coordinates relax to a new equilibrium position after a vertical Franck–Condon transition to the excited state energy potential surface.²⁴ During this process, the reorganization energy is dissipated in the protein environment. While the transfer times of an ensemble of individual disorder realizations are randomly distributed around some average value,²² the reorganization process is a systematic effect pertaining to the dynamics in all realizations in the same way. Also the ratio of chlorophyll a versus chlorophyll b affects energy gradients. For example, Chlamydomonas-LHCII, found in species of green algae, does have a reduced chlorophyll a/b ratio when compared to the LHCII of higher plants, implying a blue shift in the absorption spectra.²⁵

Due to the lack of the computational capability to carry out accurate calculations of the exciton dynamics, previous simulations of transfer time-scales in light-harvesting complexes (LHCs) employ a combined modified Redfield/generalized Förster rate equation approach.^22,26–29 However, the combined modified Redfield/generalized Förster lacks the ability to simulate the reorganization process. In addition, those models provide an ad hoc description of dynamic localization, and depend on an empirical cut-off parameter. Recently, a non-Markovian (ZOFE) quantum master equation description has been employed to investigate the robustness of transfer efficiency and the importance of vibrationally enhanced transfer in PS-II.³⁰ Here, we perform accurate simulations based on the hierarchically coupled equations of motion approach (HEOM),^31–35 which accurately incorporates the reorganization process and works for a wide parameter range for the coupling strength to the environment.

Since the computational complexity of HEOM scales exponentially with increasing system size, novel algorithms based on optimized parallelization schemes have been developed.^36–38 The most efficient implementation³⁶ employs the high computational throughput provided by modern graphics processing units (GPUs), for which a cloud computing version is hosted on http://nanohub.org.³⁹ GPU-HEOM is bound to the available GPU memory, and simulations are limited to intermediate-sized systems. Here, we overcome the memory limitation by using QMaster,³⁸ which runs on various hardware architectures including GPUs and high memory multi-core CPU architectures. We make use of the large CPU memory to benchmark the convergence of the hierarchy depth and use the high computational throughput of the GPUs for production runs.

In Section 2 we outline the structure of the Frenkel exciton model for energy transfer in C₂S₂M₂. The technical aspects of the HEOM approach are stated in Section 3. After that, we continue with the discussion of the time-scales of inter-protein transfer in the PS-II supercomplex (see Section 4). Finally, we investigate the impact of variations in the energetic structure of C₂S₂M₂ on the transfer pathways and the transfer efficiency.

2 Exciton model

The orientation of the individual proteins of the C₂S₂M₂ supercomplex is determined by a projection map at 12 Å resolution.⁵ The C₂S₂M₂ supercomplex, the structure of which is depicted in Fig. 1a, comprises four LHCII trimers, six minor light-harvesting complexes and a dimeric PS-II core complex. Absorbed light in the outer LHCII antenna complexes is transferred via the minor complexes CP24, CP26, and CP29 to CP47 and CP43 of PS-II. The transfer process is completed by irreversible charge separation triggered in the RC core. Electron transfer in the RC core is described phenomenologically by radical pair states RP1, RP2 and RP3.^12,40 We assume that the primary step of charge separation is initiated through the electronically excited core pigment Chl_D1 (the location of Chl_D1 in PS-II is illustrated in Fig. 1b) and described by the rate equation

Fig. 1 (a) Sketch of the protein structure of the C₂S₂M₂ supercomplex. The multiprotein complex contains 4 LHCII trimers and the minor complexes CP24, CP26 and CP29 which are connected to the PS-II.⁵ (b) Assembly of the pigments of PS-II, composed of CP43, CP47 and the RC core. The primary step of charge separation is initiated through excitation of pigment Chl_D1 (see eqn (1)).

We neglect backward rates since fluorescence decay lines suggest that charge recombination occurs on a much slower time scale than primary electron transfer.¹² Within this limit we model primary charge separation as irreversible exciton trapping. In the literature, more sophisticated models are also discussed, which include multiple pathways of charge separation.^41,42

We describe energy transfer in the C₂S₂M₂ supercomplex within a Frenkel exciton Hamiltonian, for which we assume that only one of the pigments is excited at once. The Hamiltonian of the single exciton manifold reads

Here |m〉 denotes the state in which pigment m is excited while the other pigments remain in the electronic ground state. For the inter-site couplings J_mn we distinguish between intra-complex and inter-complex coupling terms, depending on whether or not pigments m and n are located within the same protein. We use the same parameters for the exciton Hamiltonian as in a previous study by Bennett et al. in ref. 22, in which the Hamiltonian for C₂S₂M₂ is constructed as follows: the exciton system for the individual proteins of the LHCII-trimer, CP43, CP47, and the RC-core are taken from the literature,^11–14,43 while the inter-protein couplings are calculated using a dipole–dipole approximation. Since the structures of the minor complexes CP24 and CP26 have not yet been resolved, they are substituted by LHCII monomers. Although the exciton dynamics within the minor complexes are shown to be roughly the same as the dynamics within LHCII,^11,44 we remark that this replacement may affect the transfer pathways. In ref. 22, CP29 is also replaced by an LHCII monomer (without pigment Chl 605). However, recently the Hamiltonian for CP29 has been resolved.¹⁵ In order to isolate how much the new CP29 Hamiltonian influences transfer and to compare the HEOM results with previous approximate modified Redfield/generalized Förster simulations,²² we carry out calculations for both models: (i) with the CP29 Hamiltonian and (ii) with the LHCII monomer substitution.

The pigments are coupled to the protein environment modeled by a set of independent harmonic oscillators

and we assume a linear coupling of the exciton system to the vibrations

The reorganization energy with is added to the exciton Hamiltonian eqn (2). We define the site energies as ε_m = ε⁰_m + λ_m. The phonon mode-dependent coupling strength is captured by the spectral density

Frequently, the reorganization energy and the spectral density are assumed to be site independent. However, for the C₂S₂M₂ supercomplex each individual protein has its own reorganization energies and own form for the spectral density.²² The spectral density for LHCII is extracted from fluorescence line narrowing spectra. Since the experimental spectra cannot distinguish between the Chla and Chlb pigments, we assume the same spectral density, composed of 48 vibrational peaks, for both.^45,46 Transfer times are not much affected by the structures in the spectral density, and a coarse-grained Drude–Lorentz spectral density is appropriate to describe energy transfer in LHCII.³⁸ More details about how the parameters for the coarse-grained spectral density are obtained can be found in ref. 38. Microscopic details for the spectral densities of the minor complexes and PS-II are not known. The structure of CP29 is similar to that of an LHCII monomer. Thus, we assume that the spectral density of CP29 can be substituted with the LHCII spectral density.^15,22 For CP47 and the RC core pigments, λ = 38.64 cm⁻¹ and λ = 50.23 cm⁻¹ are respectively suggested as reasonable values for the reorganization energy.¹² For the RC core pigments a higher reorganization energy is also discussed.^41,47 The explicit form and parameters for the spectral densities used in this manuscript are given in the ESI.†

We include the primary step of charge separation phenomenologically as irreversible population trapping, which we incorporate by anti-Hermitian parts in the Hamiltonian

[script letter H]_trap = −iħΓ_RP1/2|Chl_D1〉〈Chl_D1|, (6)

where Γ_RP1 defines the rate of the primary charge separation. In a similar way we incorporate exciton losseswhere we assume exciton lifetimes of (Γ_loss)⁻¹ = 2 ns. We characterize transfer properties by the transfer efficiencyand average transfer time

For numerical evaluations, we replace the upper integration limit by t_max, which is chosen such that the total population within the pigments of the C₂S₂M₂ supercomplex has dropped below 0.001.

3 Method

We evaluate the exciton dynamics within the hierarchically coupled equation of motion (HEOM) method.^31,34,35 HEOM is an open quantum system approach which treats the coupling to the vibrational modes as a bath. The time evolution of the total density operator R(t), which characterizes the degrees of freedom of the exciton system as well as the ones of the phonon bath, is governed by the Liouville equation

We assume that the total density operator R(t) factorizes at initial time t₀ = 0 in excitonic and vibrational degrees of freedom R(t₀) = ρ(t₀) ⊗ ρ_phon(t₀). We trace out the vibrational degrees of freedom, and get the time evolution for the reduced density matrix ρ(t) describing the exciton degrees of freedom only:

By employing a second order cumulant expansion, using a Drude–Lorentz spectral density in combination with a high temperature approximation ħγ_m/k_BT < 1, we cast the time non-local eqn (11) into a hierarchy of coupled time local equations of motion:

where we define ρ(t) = σ^(0,…,0)(t), V^×_m• = [V_m,•], V^o_m• = {V_m,•}, and V_m = |m〉〈m|. The hierarchy can be truncated for a sufficiently large hierarchy level . Convergence of the hierarchy can be tested by comparing deviations in the dynamics with increasing truncation level.

To increase the accuracy of the high temperature approximation (HTA) of HEOM,⁴⁸ we include additional correction terms⁴⁹ for which we replace

Here ν₁ is defined as the first Matsubara frequency ν₁ = 2π/βħ.

For structured spectral densities, a similar hierarchical expansion has been derived which relies on a decomposition of the spectral density in terms of shifted Drude–Lorentz peaks^50,51 or underdamped Brownian oscillators.³⁵ The accuracy of the improved high temperature correction can be validated by various means. For example, convergence could be tested by including more Matsubara frequencies. However, the slow convergence of the Matsubara expansion renders this approach numerically challenging. Less numerically involved approaches include the validation of the stationary state, which is expected to follow a thermal Boltzmann distribution for the parameter regime of C₂S₂M₂. Also, comparison of monomer line-shapes computed with HEOM to analytically known results have been proposed as an appropriate test for the HTA.³⁸

4 Discussion

Together with LHCII complexes, the C₂S₂M₂ supercomplex aggregates as a large photosynthetic network in the grana membrane. For each C₂S₂M₂ supercomplex, there are about six additional loosely bound LHCII trimers,⁴ which form a large antenna system with densely packed chlorophylls. Energy is either absorbed in the pool of loosely bound LHCII trimers and then transferred to one of the peripherical LHCIIs of the C₂S₂M₂ supercomplex or absorbed directly in the LHCII trimers of C₂S₂M₂. Further, to some extent, energy is absorbed in the minor complexes and PS-II. We expect that the contribution of light absorption in the minor complexes and PS-II to the photosynthetic yield is of less importance, since most of the photoactive area in the grana membrane is covered by the LHCIIs. Thus, to reach a high photosynthetic yield, fast and efficient transfer from the LHCIIs toward the RC core pigments of PS-II becomes indispensable.

In the following, we investigate average transfer times and the efficiency of energy transfer from the peripherical LHCII-m monomeric unit, labeled as unit-1 in Fig. 1, to the reaction center in which the primary step of charge separation takes place. Since the simulation of the complete C₂S₂M₂ supercomplex with HEOM is beyond the current capabilities of QMaster, we employ a certain amount of symmetry along the x-axis and y-axis and reduce the system to a multi-protein network composed of LHCII-m, CP24, CP29, CP47 and the RC-core, comprising 93 pigments in total. This reduction does not take into account transfer pathways between neighboring quadrants of the C₂S₂M₂ supercomplex, which could open additional transfer pathways.

4.1 Energy gradient drives directionality

First, we monitor the population dynamics in the absence of trapping and energy losses. We highlight how energy spreads among the different protein complexes, which, as we will analyze in detail, is driven by energy gradients in the pigments of the C₂S₂M₂ supercomplex. Further, we explore the influence of the reorganization process on the exciton dynamics.

In the following, we investigate the deficiency of the approximate combined modified Redfield/generalized Förster method by comparing the population dynamics obtained from the combined method with the HEOM results. The combined modified Redfield/generalized Förster approach divides the exciton Hamiltonian eqn (2) into a strongly coupled part [script letter H]_strong (associated with strongly coupled domains) and a weakly coupled part. Hereby, [script letter H]_strong is given by strongly coupled clusters with inter-site couplings J_nm larger than or equal to a certain threshold value. We follow ref. 11 and use a threshold of 15 cm⁻¹. The intra-domain dynamics is then modeled by modified Redfield, while the inter-domain transfer is described by generalized Förster theory. We use the same implementation for the combined modified Redfield/generalized Förster approach as in ref. 38. Therein, transfer time-scales of the Chlb/Chla inter-band relaxation in an LHCII monomer are discussed, and it is shown that the results of the combined modified Redfield/generalized Förster approach are in reasonable agreement with the HEOM results. Since the choice of initial conditions of the combined method is restricted to eigenstates of certain domains in [script letter H]_strong, we set the highest energy state of the domain which predominantly populates pigment Chlb 606 of the LHCII-m unit-1 as the initial condition. To allow for comparison, we use the same initial condition for the HEOM calculations.

Fig. 2a depicts the aggregated population at the individual protein complexes obtained within HEOM. Convergence of the hierarchy depth is verified by comparison with a higher truncation level. Further, the stationary state is in good agreement with a thermal Boltzmann distribution, which supports that the high temperature approximation eqn (13) is valid. Details about the convergence of HEOM are given in the ESI.† Overall energy transfer and directionality are driven by energy relaxation along energy gradients within the C₂S₂M₂ supercomplex. LHCII and the minor complexes (modeled by LHCII monomers without Chl 605) exhibit the highest energies, while CP47 and the RC core pigments are lower in energy. The exciton, initially located at the unit-1 LHCII-m monomer, spreads over the complete LHCII-m trimer and populates the minor complexes. The fast initial spread shows as maxima in the aggregated populations at LHCII-m units-2 and 3. The highest population at unit-3 is obtained after about 18 ps, while the maximum population at unit-2 is reached a bit later at about 43 ps. The high population of LHCII-m unit-2 indicates that energy transfer does not exclusively proceed along pathways corresponding to the shortest distance to PS-II and the RC. The minor complexes are populated on a similar time-scale as the monomeric LHCII-m units. On a slower timescale, CP47 and the RC core of PS-II get populated, and finally after about 250–350 ps the system reaches a steady state in which energy relaxation drives the population to the energetically low exciton states at CP47. A schematic sketch summarizing the rough estimates for the energy transfer time-scales is given in Fig. 2d.

Fig. 2 Aggregated populations at T = 277 K in absence of trapping in the 93 site network comprising LHCII-m, CP24, CP29, CP47 and the RC-core. The initial state is given by the highest exciton state within the domain of [script letter H]_strong, which dominantly excites pigment Chlb 606 of the LHCII-m unit-1. Depicted are the population dynamics within the (a) HEOM and (b) combined modified Redfield/generalized Förster approach. (c) Sketch of the layout of the exciton energy bands before (dashed boxes) and after (solid boxes) dissipation of the reorganization energy. (d) Illustrates rough estimates for the time scales of how energy spreads across the individual proteins.

The dynamics of the combined modified Redfield/generalized Förster approach (Fig. 2b) diverges from the HEOM results in several aspects. Overall relaxation is overestimated by the combined modified Redfield/generalized Förster approach. In particular, transfer to LHCII-m unit-2 is about seven times faster, and already at 6 ps there is 0.29 population at unit-2. Further, unit-2 gets populated ahead of unit-3. Therefore, the pathway of how energy spreads over the monomeric units of the LHCII-m trimer is reversed when compared to the HEOM calculation, and thus the combined method does not predict reliable pathways of energy flow during the first picoseconds. However, the main difference is in the stationary population, which is not only approached faster (at about 150–250 ps) but predicts a much higher aggregated population at CP47 and the RC. The combined modified Redfield/generalized Förster approach overestimates the efficiency of energy funneling towards the PS-II.

To understand the discrepancy in the stationary population, we need to investigate the energetic layout of the C₂S₂M₂ supercomplex. The stationary state typically follows a thermal Boltzmann distribution. However, the situation becomes more complicated in the presence of the reorganization process, in which the reorganization energy dissipates during the dynamics, which modifies the energetic layout and affects the thermal population. The boxes in Fig. 2c indicate the extension of the exciton bands of the isolated proteins. The dashed lines correspond to the situation where the site energies comprise the bare excitation energy plus the reorganization energy. Due to the reorganization process, the energetic structure changes during the dynamics and the energy of the proteins is lowered by the reorganization energy. In particular, the band of the monomeric LHCII-m units and the band of the minor complexes shift to lower energies, while the small reorganization energies at CP47 and RC induce only minor modifications. In total, the already flat energy gradient gets even more flattened. This has a significant impact on the thermal population. Without the reorganization process (dashed lines), we expect a thermal population of about 0.75 at the pigments of PS-II. Taking into account the reorganization process (solid line) reduces the efficiency of energy funneling and only a population of 0.56 accumulates at PS-II. Our analysis is consistent with the findings for the population dynamics and explains the strong deviations in the stationary state between HEOM and the combined modified Redfield/generalized Förster method. We note that for the combined modified Redfield/generalized Förster method, the effect of the reorganization energy on the thermal population can be corrected by subtracting the reorganization energy from the exciton Hamiltonian prior to the dynamics. This is based on the assumption that the reorganization energy dissipates on an infinitely fast time-scale. With this inclusion of the reorganization process, the combined modified Redfield/generalized Förster approach gives the correct thermal population. Nevertheless, the time-scale of energy relaxation is still overestimated. Further, despite the inclusion of the bath reorganization, the combined modified Redfield/generalized Förster approach does not give an accurate prediction of energy transfer time-scales in the early stages of the transfer process (first ∼ 50 ps). For more details we refer the reader to the ESI.†

4.2 The reorganization process affects transfer efficiency

In the following, we investigate how minor changes in the energetic layout of the C₂S₂M₂ supercomplex influence transfer properties such as transfer efficiency and average transfer time. As we have discussed in detail in the previous section, one mechanism is the reorganization process. Here, we continue the discussion and examine how much the reorganization process affects transfer efficiency. Another aspect is the influence of the replacement of the minor complexes with LHCII monomeric units on the transfer properties. For instance, the recently derived exciton Hamiltonian of CP29 shows various differences from the exciton system of an LHCII monomer.¹⁵

We incorporate the primary step of charge separation by irreversible energy trapping as described in Section 2. Different values for the rate constant of primary charge separation Γ_RP1 have been predicted from fits to fluorescence decay lines, ranging from Γ_RP1⁻¹ = 0.1 ps (ref. 12) to Γ_RP1⁻¹ = 0.64 ps.²² Pump-probe spectra predict even larger time constants for the pheophytin (Pheo) reduction of about 3 ps.⁴⁰ We do not explicitly take into account mechanisms of photoprotection and quenching and phenomenologically describe exciton losses by assuming an exciton lifetime of (Γ_loss)⁻¹ = 2 ns.

In the following, we carry out HEOM simulations, in which we include trapping and energy losses. To investigate the effects of the reorganization process on the energy transfer times, we slightly modify the Hamiltonian of the C₂S₂M₂ supercomplex in a benchmark calculation, for which we artificially restore the original energy gradient across the pigment proteins by adding the reorganization energy of 220 cm⁻¹ to the site energies of LHCII and the minor complexes. We neglect the minor energetic changes induced by the reorganization process at the pigments of CP47 and the RC and denote the modified Hamiltonian as [script letter H]_{ex,add reorg.}. Relaxation time scales in the population dynamics are hardly affected by the shifts in the site-energies, but the thermal state adjusts now according to the modified energy gradient. For [script letter H]_{ex,add reorg.} we obtain a similar thermal state in the population dynamics with high population at the PS-II pigments (0.81), as predicted by the calculations with the combined modified Redfield/generalized Förster method. The small deviations largely result from the fact that we did not add additional reorganization energies to the site energies of CP47 and RC.

The transfer time as a function of trapping rate follows a linear trend for the considered parameter regime, as illustrated in Fig. 3. We assume that initially the exciton is located at LHCII-m unit-1, and we populate the initial density matrix according to eigenstates of the isolated LHCII monomeric unit. The shown results correspond to transfer times averaged over all 14 exciton states used as the initial conditions. The chosen initial conditions show similar transfer times, which is also reflected by the small standard deviations (5–6 ps). Transfer times (efficiency) for the C₂S₂M₂ supercomplex (marked by the red circles) are in the range between 242 ps (88.0%) and 302 ps (85.2%), depending on the trapping rate Γ_RP1. [script letter H]_{ex,add reorg.} exhibits more efficient energy funneling towards the pigments at PS-II, and therefore transfer is faster by about 36 ps to 53 ps. Previous calculations based on the combined modified Redfield/generalized Förster method predict transfer times of about 200 ps for transfer from the peripherical domains to the Chl_D1 in the RC.²² This is in good agreement with our results for [script letter H]_{ex,add reorg.}, which yields a transfer time of 211 ps for a trapping rate of Γ_RP1 = 0.5 ps, which is similar to the one used in ref. 22.

Fig. 3 Trapping time evaluated for various rate constants of primary charge separation Γ_RP1 at T = 277 K. The transfer time is given as an average over different initial conditions corresponding to eigenstates of the isolated LHCII-m unit-1 monomer. The error bars mark the standard deviations. We investigate changes in the transfer time induced by structural changes in the Hamiltonian. We compare three different scenarios, (i) [script letter H]_{ex,CP29-Renger}, for which we use the CP29 Hamiltonian of Renger et al., ref. 15, (ii) [script letter H]_ex, for which the CP29 is substituted by an LHCII monomer (without Chl 605) and (iii) [script letter H]_{ex,add reorg.}, for which we add the reorganization energy to the site energies of LHCII-m and the minor complexes of [script letter H]_ex.

4.3 Minor complex CP29 guides transfer

Since the reorganization process already alters the energy transfer efficiency, we expect that the substitution of the minor complexes by LHCII monomers may also significantly affect the energy transfer properties. In this section, we use the Hamiltonian of CP29 derived by Renger et al., ref. 15, instead of the LHCII monomer replacement. We denote the new Hamiltonian as [script letter H]_{ex,CP29-Renger}, while the previous situation with the LHCII monomer substitution is referred to as [script letter H]_ex.

For [script letter H]_{ex,CP29-Renger}, the pigments of CP29 form a lower energy band than the LHCII monomers. This has a two-fold implication for the transfer process. Firstly, the energy gradient between the outer LHCII antenna and the minor complex CP29 gives rise to an additional grade of directionality and supports fast transfer from the peripherical LHCII-m trimer to CP29. The minor complex CP29 presumably acts as an exit marker, which guides energy from the outer antenna towards pigments closer to the reaction center. Secondly, the pigments of CP29 and CP47 form a spatially extended region of low-energy states, and hence energy accumulates at the pigments in proximity to the RC, while the final transfer step to the RC core pigments is energetically uphill and therefore slow. Overall, the two effects result in a slightly slower energy transfer within the C₂S₂M₂ supercomplex while including the CP29 Hamiltonian, see Fig. 3. For large trapping rates, Γ_RP1 > 1 ps, the slow down of the energy transfer gets more pronounced.

Fig. 4 charts snapshots of the exciton dynamics. The upper (lower) panels correspond to [script letter H]_{ex,CP29-Renger} ([script letter H]_ex). The radius of the colored spheres represents the population at each pigment. For better visualization we use an arctan scale. The spheres are uncolored if the population remains below 0.0079. Initially, the highest eigenstate of LHCII-m unit-1 is excited. Both Hamiltonians show a fast spread of the energy, and at 12 ps the energy distributes across the whole LHCII-m trimer. While [script letter H]_ex distributes population equally among the minor complexes, [script letter H]_{ex,CP29-Renger} yields a more directed energy transfer towards CP29 and CP47. For longer times of 150 ps, energy accumulates at the low energy states at CP29 and CP47 for [script letter H]_{ex,CP29-Renger} and thus forms a bottleneck for transfer to the RC. The bottleneck is less pronounced for [script letter H]_ex. The RC pigments do not show significant population at any time, since as soon as energy enters the RC there is fast transfer to Chl_D1 and the fast time-scale of primary charge separation leads to the trapping of the population.

Fig. 4 Snapshots of the exciton dynamics in the presence of population trapping in the reaction center (Γ_RP1⁻¹ = 0.5 ps) at T = 277 K. The spheres represent the positions of the individual pigments while the radii reflect the population at each pigment (we employ an arctan scale). The color signifies the population at the individual protein complexes. Different colors are used for the individual protein complexes (the color encoding is analogues to Fig. 2). Spheres in gray indicate pigments with less than 0.0079 population. As the initial condition, we use the highest eigenstate of the isolated LHCII-m unit-1. The upper and lower panels show the results for the two different Hamiltonians, [script letter H]_{ex,CP29-Renger} and [script letter H]_ex, respectively. Both differ in the structure of the minor complex CP29.

The rate limiting step in the transfer chain is the energetically up-hill transfer to the RC core. This is illustrated best in the aggregated population dynamics in the presence of trapping in the reaction center, Fig. 5. We obtain a fast decay of population in LHCII-m and after 100 ps, more than 0.75 of the population has left the LHCII-m trimer. At the same time, about 0.44 of the population accumulates in CP29 and CP47. After 300 ps, 0.2 of the population still remains at CP29 and CP47.

Fig. 5 Aggregated populations for [script letter H]_{ex,CP29-Renger} in the presence of trapping (Γ_RP1 = 0.5 ps) at 277 K. Energy accumulates at low-energy bottleneck states at CP29 and CP47, limiting transfer to the RC.

5 Conclusion

With QMaster, a high-performance implementation of the HEOM method, accurate calculations of excitonic energy transfer in multi-protein photosynthetic functional units, such as the C₂S₂M₂ supercomplex, become feasible. We investigate the transfer times and transfer efficiency of energy conversion within the primary step of charge separation.

The general concept behind energy transfer in C₂S₂M₂ is given by energy relaxation. Due to the flat energy gradient across the proteins, even small variations induced by the reorganization of the molecular coordinates within the excited potential energy surface affect the energy transfer process. The impact of the reorganization process is rather significant, and energy relaxation drives much less population to CP47 and the RC than expected from the site energies of the Hamiltonian. The reorganization process induces a noticeable drop in the transfer efficiency of about 1.8% to 2.6% in absolute numbers for a 2 ps exciton lifetime, and thus cannot be neglected in simulations of energy transfer in large multi-protein complexes.

Our simulations suggest that the minor complex CP29 acts as an exit marker and adds directionality to the energy transfer from the peripherical LHCII to the proteins in proximity to the RC core. The C₂S₂M₂ supercomplex is not optimized for efficient transfer. Energy accumulates in low energy states at CP29 and CP47, while the final transfer step needs to overcome an energy barrier and therefore is slow. Thus, the energy transfer exhibits the character of a transfer-to-trap limited model. In conclusion, within our model, we show that the structural layout of C₂S₂M₂ is not optimized for efficient transfer and suggests that photoprotection considerations are very relevant. The extension of accurate HEOM models to this case is possible and is a promising direction for future research. Our accurate simulations provide a first step toward a consistent description of energy transfer in multi-protein networks such as C₂S₂M₂. The currently used parameter sets for the individual proteins are largely obtained empirically by fitting simulations based on approximate methods to experimental spectra, and a refinement of the individual Hamiltonians may be needed. To this end, future works need to focus on the comparison of the accurate HEOM calculations with experimental data, such as absorption spectra, circular dichroism spectra or hole-burning data.⁵²

Acknowledgements

The authors would like to thank Dr Kapil Amernath and Dr Doran Bennett for helpful discussion. We thank Dr Doran Bennett for providing the atomistic structure and the Hamiltonian of C₂S₂M₂. This work was supported as part of the Center for Excitonics, an Energy Frontier Research Center funded by the U. S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0001088. We thank Nvidia for support via the Harvard CUDA Center of Excellence. The computations in this paper were run on the Odyssey cluster, supported by the FAS Division of Science, Research Computing Group at Harvard University.

References

1. J. P. Dekker and E. J. Boekema, Biochim. Biophys. Acta, Bioenerg., 2005, 1706, 12–39.
2. M. P. Johnson, T. K. Goral, C. D. Duffy, A. P. Brain, C. W. Mullineaux and A. V. Ruban, Plant Cell, 2011, 23, 1468.
3. R. Kouřil, J. P. Dekker and E. J. Boekema, Biochim. Biophys. Acta, Bioenerg., 2012, 1817, 2.
4. C. Duffy, L. Valkunas and A. Ruban, Phys. Chem. Chem. Phys., 2013, 15, 18752.
5. S. Caffarri, R. Kouřil, S. Kereïche, E. J. Boekema and R. Croce, EMBO J., 2009, 28, 3052.
6. J. P. Dekker and R. Van Grondelle, Photosynth. Res., 2000, 63, 195.
7. B. A. Diner and F. Rappaport, Annu. Rev. Plant Biol., 2002, 53, 551.
8. J. Kern and G. Renger, Photosynth. Res., 2007, 94, 183.
9. S. Caffarri, K. Broess, R. Croce and H. van Amerongen, Biophys. J., 2011, 100, 2094.
10. J. Huh, S. K. Saikin, J. C. Brookes, S. Valleau, T. Fujita and A. Aspuru-Guzik, J. Am. Chem. Soc., 2014, 136, 2048.
11. V. Novoderezhkin, A. Marin and R. van Grondelle, Phys. Chem. Chem. Phys., 2011, 13, 17093.
12. G. Raszewski and T. Renger, J. Am. Chem. Soc., 2008, 130, 4431.
13. G. Raszewski, B. A. Diner, E. Schlodder and T. Renger, Biophys. J., 2008, 95, 105–119.
14. F. Müh, M. E.-A. Madjet and T. Renger, Photosynth. Res., 2012, 111, 87.
15. F. Müh, D. Lindorfer, M. S. am Busch and T. Renger, Phys. Chem. Chem. Phys., 2014, 16, 11848.
16. K. Broess, G. Trinkunas, C. D. van der Weij-de Wit, J. P. Dekker, A. van Hoek and H. van Amerongen, Biophys. J., 2006, 91, 3776.
17. K. Broess, G. Trinkunas, A. van Hoek, R. Croce and H. van Amerongen, Biochim. Biophys. Acta, Bioenerg., 2008, 1777, 404.
18. R. Croce and H. van Amerongen, J. Photochem. Photobiol., B, 2011, 104, 142.
19. J. Chmeliov, G. Trinkunas, H. van Amerongen and L. Valkunas, J. Am. Chem. Soc., 2014, 8963.
20. Y. Miloslavina, M. Szczepaniak, M. Müller, J. Sander, M. Nowaczyk, M. Rögner and A. Holzwarth, Biochemistry, 2006, 45, 2436.
21. G. Tumino, A. P. Casazza, E. Engelmann, F. M. Garlaschi, G. Zucchelli and R. C. Jennings, Biochemistry, 2008, 47, 10449.
22. D. I. Bennett, K. Amarnath and G. R. Fleming, J. Am. Chem. Soc., 2013, 135, 9164.
23. Y. Shibata, S. Nishi, K. Kawakami, J.-R. Shen and T. Renger, J. Am. Chem. Soc., 2013, 135, 6903.
24. V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, Wiley-VCH, Weinheim, 2004.
25. B. Drop, M. Webber-Birungi, K.-N. Yadav Sathish, A. Filipowicz-Szymanska, F. Fusetti, J. Boekema Egbert and R. Croce, Biochim. Biophys. Acta, 2014, 1837, 63.
26. T. Renger, M. Madjet, A. Knorr and F. Müh, J. Plant Physiol., 2011, 168, 1497.
27. M. Yang and G. R. Fleming, Chem. Phys., 2002, 282, 163.
28. V. Novoderezhkin, Biochemistry (Moscow) Supplement Series A: Membrane and Cell Biology, 2012, vol. 6, p. 314.
29. D. Zigmantas, E. L. Read, T. Mančal, T. Brixner, A. T. Gardiner, R. J. Cogdell and G. R. Fleming, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 12672.
30. J. J. J. Roden, D. I. G. Bennett and K. B. Whaley, axiv preprint, arXiv:1501.06674.
31. Y. Tanimura and R. Kubo, J. Phys. Soc. Jpn., 1989, 58, 101.
32. Q. Shi, L. Chen, G. Nan, R.-X. Xu and Y. Yan, J. Chem. Phys., 2009, 130, 084105.
33. J. Hu, M. Luo, F. Jiang, R.-X. Xu and Y. Yan, J. Chem. Phys., 2011, 134, 244106.
34. A. Ishizaki and G. R. Fleming, J. Chem. Phys., 2009, 130, 234111.
35. Y. Tanimura, J. Chem. Phys., 2012, 137, 22A550.
36. C. Kreisbeck, T. Kramer, M. Rodríguez and B. Hein, J. Chem. Theory Comput., 2011, 7, 2166.
37. J. Strümpfer and K. Schulten, J. Chem. Theory Comput., 2012, 8, 2808.
38. C. Kreisbeck, T. Kramer and A. Aspuru-Guzik, J. Chem. Theory Comput., 2014, 10, 4045.
39. C. Kreisbeck and T. Kramer, Exciton Dynamics Lab for Light-Harvesting Complexes (GPU-HEOM), 2013.
40. A. Holzwarth, M. Müller, M. Reus, M. Nowaczyk, J. Sander and M. Rögner, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 6895.
41. V. I. Novoderezhkin, E. Romero, J. P. Dekker and R. van Grondelle, ChemPhysChem, 2011, 12, 681.
42. E. Romero, R. Augulis, V. I. Novoderezhkin, M. Ferretti, J. Thieme, D. Zigmantas and R. van Grondelle, Nat. Phys., 2014, 10, 676–682.
43. G. Raszewski, W. Saenger and T. Renger, Biophys. J., 2005, 88, 986.
44. A. Marin, F. Passarini, R. Croce and R. van Grondelle, Biophys. J., 2010, 99, 4056.
45. V. I. Novoderezhkin, M. A. Palacios, H. van Amerongen and R. van Grondelle, J. Phys. Chem. B, 2004, 108, 10363.
46. V. I. Novoderezhkin and R. van Grondelle, Phys. Chem. Chem. Phys., 2010, 12, 7352.
47. V. I. Novoderezhkin, E. G. Andrizhiyevskaya, J. P. Dekker and R. van Grondelle, Biophys. J., 2005, 89, 1464.
48. J. Olsina, T. Kramer, C. Kreisbeck and T. Mancal, J. Chem. Phys., 2014, 141, 164109.
49. A. Ishizaki and G. R. Fleming, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 17255.
50. C. Kreisbeck and T. Kramer, J. Phys. Chem. Lett., 2012, 3, 2828.
51. C. Kreisbeck, T. Kramer and A. Aspuru-Guzik, J. Phys. Chem. B, 2013, 117, 9380.
52. R. Jankowiak, M. Reppert, V. Zazubovich, J. Pieper and T. Reinot, Chem. Rev., 2011, 111, 4546.

---

Footnote

† Electronic supplementary information (ESI) available: The ESI contains the parameters for the used spectral densities and details the convergence of HEOM. See DOI: 10.1039/c5sc04296h

---