Modeling of excitation dynamics in photosynthetic light-harvesting complexes: exact versus perturbative approaches

The primary processes of photosynthesis occur in light-harvesting (antenna) complexes that absorb sunlight and transfer the excitation energy to the reaction centers (RC), where a charge separation is initiated [1–3]. The light-harvesting complexes are essentially heterogeneous, i.e. they are characterized by a big spread of pigment–pigment distances. They may contain monomeric pigments and tightly packed aggregates with strong excitonic interactions (dimers, trimers, etc). The excitation dynamics within strongly coupled clusters has the form of Redfield-type relaxation between delocalized exciton states, whereas intercluster migration resembles Förster-type hopping between weakly coupled compartments. In recent years, significant progress in our understanding of light harvesting in antenna complexes has been made using a theoretical approach where the modified Redfield model is combined with generalized Förster theory [4, 5]. The Redfield–Förster approach allows us to build a unified picture of the exciton structure and excitation dynamics in an antenna, where the linear and nonlinear (ps/fs) spectral responses can be modeled at a quantitative level after the structure-based exciton Hamiltonian and electron-phonon spectral density (extracted from experiment or calculated ab initio) are specified. Some of the parameters (like site energies that are difficult to calculate precisely) can be adjusted from an evolutionary-based search with a simultaneous fit of all the available spectra and kinetics. Such an approach has led to quantitative models of many antenna and core complexes, i.e. antenna complexes B800-850 [6], B800-820 [7], and B850 [8, 9] from purple bacteria, FMO complex from green bacteria [10, 11], PE545 complex from cryptophyte algae [12, 13], LHCSR3 complex from green algae [14], plant light-harvesting complexes LHCI (Lhca4) [15], LHCII [16, 17], CP29 [18, 19], core complexes of photosystem I (PSI) [4] and photosystem II (PSII) [20, 21], and the PSII reaction center (PSII-RC) [22–28].

Recently, a coherent modified Redfield (cmR) theory has been developed [29–31], containing transfers between the one-exciton populations together with the decay of the coherences. Unfortunately, nonsecular terms, i.e. transfers between coherences as well as transfers between populations and coherences, are not included. These terms are important to describe transfers between the weakly coupled and isoenergetic sites, where nonsecular transfers maintain long-lived coherences between the exciton (delocalized) eigenstates. These coherences keep the excitation localized, whereas in the secular approximation (where coherences quickly decay without being 'repumped' from the populations) initially localized excitations exhibit rather fast delocalization between the donor and acceptor states. In the site representation this looks like unrealistically fast energy transfer [32]. The cmR theory suffers from such a 'resonant artifact' similar to modified Redfield (where the coherences are not included at all). To circumvent it, the cmR can be combined with the generalized Förster theory, or with the small polaron quantum master equation [31]. In such combined approaches the most questionable point is how to determine the crossover between the two energy transfer regimes. The simplest way is to split the whole system into compartments, with the exciton couplings M_jj' exceeding some critical value, M_c [4, 5]. Then the dynamics within the clusters with M_jj' > M_c can be described with the Redfield theory, whereas the transfers between these clusters (with intercluster couplings M_jj' < M_c) are treated with the generalized Förster theory. The critical value M_c can be chosen intuitively [4, 5, 10, 16, 17, 20, 21, 33] or by comparing the results with those of exact methods, like the hierarchical equation of motion (HEOM) [32, 34–37].

Another limitation of the Redfield theory is connected with the positivity of the density matrix that generally cannot be guaranteed (depending on initial conditions, temperature, and details of the exciton–phonon coupling, like the correlation time etc) [38, 39]. This issue has been studied recently by comparing the Redfield (Bloch–Redfield) approaches to exact solutions [40, 41]. The problem has also been discussed in the context of the modeling of energy transfer in photosynthetic antenna complexes [42] and electron transfer reactions [43].

In this paper we calculate the excitation dynamics in the major light-harvesting complex LHCII from higher plants using the combined cmR-generalized Förster theory (cmRgF) and HEOM. In order to reduce the HEOM computing time we propose a simplified spectral density (in the form of a single overdamped Brownian oscillator) approximating a real spectral density and capable of reproducing the spectral profiles. In the cmRgF approach we split the complex into clusters, where the pigments j, j' within clusters are connected by strong couplings, i.e. M_jj' > M_c, implying that the biggest intercluster coupling is small, i.e. M_ij < M_c. We show that such compartmentalization depends on the energy gap E_ij between the most strongly coupled ith and jth sites from different clusters. If E_ij is small enough we have to invoke the generalized Förster theory for intercluster transfers. In the case of a big gap E_ij we can neglect nonsecular transfers and treat the intercluster dynamics with the modified Redfield theory, considering the two clusters as one (meaning that the new M_c' value is smaller than the former one M_c, because now M_ij > M_c'). Comparing the cmRgF kinetics obtained for different compartmentalization schemes (with different M_c values) and the exact solution (HEOM) we show that the optimal M_c really depends on the energy gap between the sites, and therefore can be different for different parts of the complex. We demonstrate how to define a reliable cmRgF model using the LHCII complex as an example, and we propose that a similar approach can be taken for other photosynthetic antennae and core complexes. The thus defined cmRgF schemes allow a fine adjustment of the parameters of the exciton model (site energies, disorder values, etc) using an evolutionary-based search (with a quantitative fit of linear and nonlinear spectral responses), which is difficult to do with numerically expensive exact methods.

The structure of the major light harvesting complex of higher plants (LHCII) is shown in figure 1. Each of the three subunits of the trimeric LHCII complex contains 14 chlorophylls (Chls) arranged in two layers. Analysis of the exciton couplings between the 14 Chls allows us to define several strongly coupled clusters within both the stromal- and lumenal-side layers (as shown in table 1). The arrangement of these clusters is shown in figure 2. One can expect the Redfield-type relaxation within strongly coupled clusters and the Förster-type migration between (weakly connected) clusters. The Förster theory also works in the case of large energy gaps, breaking the mixing of the donor and acceptor (even in the case of strong coupling between them). For example, in spite of the strong coupling between the Chls b and Chls a on the stromal side, they can be considered as separate clusters connected by Förster-type Chls b → a transfers.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Notice that partitioning into clusters (shown in table 1 and figure 2) is not unique. Below we analyze this scheme together with alternative ones by comparing the emerging Redfield–Förster pictures with the exact solution of the hierarchical equation of motion (HEOM) for the whole monomeric subunit.

In this paper we model the excitation dynamics in LHCII using different approaches, i.e. standard, modified, and coherent modified Redfield theories, the combined Redfield–Förster theory, and various versions of HEOM. The dynamics is determined by the spectral density of the system–bath (exciton–phonon) coupling that can be arbitrary in the Redfield (and Redfield–Förster) approach. In contrast, the HEOM method most often applies a simple spectral density in the form of an overdamped Brownian oscillator or a single underdamped vibration. Including just one more vibrational mode in the spectral density increases the number of equations immensely. Thus, using a complicated spectral density with 7–11 peaks (approximating the realistic shapes containing up to several dozen high frequency vibrations) requires powerful computing [36, 45]. In this section we define parameters of a simple spectral density (consisting of a single overdamped Brownian oscillator), that is still capable of describing quantitatively the spectral shapes and relaxation rates in LHCII at room temperature. The thus-defined spectral density can be used to compare the dynamics calculated with HEOM and various perturbative approaches.

A realistic spectral density for LHCII (based on measured low-temperature fluorescence line-narrowing profiles [46]) can be constructed as a sum of one low-frequency overdamped Brownian oscillator and 48 underdamped vibrations [47]. In figure 3 we use this spectral density to model the room-temperature absorption (OD) and fluorescence (FL) spectra of LHCII with the modified Redfield theory. We calculate the spectral profiles for one monomeric subunit, containing 14 Chls with the parameters listed in table 1. The site inhomogeneity value is σ = 90 cm⁻¹. The inhomogeneity and phonon couplings are the same for Chls a and b.

Zoom In Zoom Out Reset image size

Now we switch from the realistic to a simpler spectral density consisting of a single overdamped Brownian oscillator. In order to reproduce the spectral profiles we need the characteristic frequency (damping constant) γ_j = 500 cm⁻¹ and coupling (reorganization energy) of λ_j = 300 cm⁻¹. The spectra calculated with this spectral density and using the modified Redfield theory are shown in figure 3. Obviously, the fit is not as good for such a simplified model. In particular, a single Brownian oscillator is not sufficient to reproduce the 740 nm FL sub-band determined by intense vibrations around 1200 cm⁻¹. On the other hand, the FL profile in the region from 680 to 720 nm (corresponding to frequencies up to 800 cm⁻¹) is reproduced well. Bearing in mind that the biggest exciton splittings in LHCII (corresponding to Chls b → a transfers) do not exceed 600 cm⁻¹, we conclude that our simplified spectral density should give realistic rates for both intra- and inter-band transfers. Notice that the reorganization energy with such a spectral density (i.e. 300 cm⁻¹ in the site representation) is less than the reorganization with the original spectral density (which is about 480 cm⁻¹ due to the presence of intense high-frequency vibrations). To this end, all the site energies have been red-shifted by about 120 cm⁻¹ (to compensate for the reduced reorganization shift of the exciton components) in order to reproduce the observed OD/FL peak positions (grey curve in figure 3). Below, we use this single-component spectral density and the shifted site energies in all the calculations performed both with the cmRgF and HEOM approaches. Using a simplified spectral density makes it possible to calculate the kinetics with HEOM with routinely available computational facilities.

In all numerical examples (presented in the following sections) we use the scaled hierarchical equation of motion developed by Zhu et al [49]. The main equations are given in appendix A. We have found that the scaled HEOM gives converging results already at small cutoff values N_c, i.e. one can restrict to N_c = 3–4 (for N_c > 4 there are no visible changes in the kinetics). There are some small but visible difference between the K = 0 and K = 1 cases (where K is the number of temperature correction terms), whereas switching to K = 2 does not change the results any more. Therefore we calculate the scaled HEOM kinetics with N_c = 3, K = 1 (for small fragments of the complex consisting of N = 5–8 molecules) and N_c = 3, K = 0 (for the whole monomeric subunit with N = 14).

Notice that quantitatively the same kinetics can be obtained with the original (nonscaled) HEOM, but in this case bigger cutoff values are required, i.e. N_c = 7–10. We have found that the original version of HEOM with the lowest-order temperature correction [50] gives the same results as the scaled HEOM with K = 0, whereas the original HEOM in the high temperature limit [51] gives sizable deviations (which can be expected, because our damping constant significantly exceeds the room temperature thermal energy).

Below we compare the HEOM solution with various versions of the Redfield theory (including combined Redfield–Förster approaches). It should be noted that we did not find any violation of the positivity of the density matrix in our numerical examples, where we have checked various initial conditions and various system–bath correlation times (by changing the damping constant γ_j).

4.1. Equilibration within the stromal-side Chls a (N = 5 model)

As a first example we consider the dynamics within the two stromal-side Chls a clusters, i.e. a602-603 and a610-611-612 upon excitation of the blue-most a603 site (figure 4). We use the spectral density in the form of a single overdamped Brownian oscillator with γ_j = 500 cm⁻¹ and various reorganization energies, i.e. λ_j = 10, 30, 100 and 300 cm⁻¹. The population dynamics is calculated for one realization of the disorder (with the unperturbed site energies from table 1). For each reorganization energy value we compare site populations calculated with four theories, i.e. standard Redfield (sR) with the full relaxation tensor, coherent modified Redfield theory (cmR) in the Markovian (time-independent) limit (see appendix B), time-independent coherent modified Redfield theory combined with generalized Förster theory (cmRgF), and scaled HEOM. Some of the coherences between the sites are shown as well, i.e. intra-cluster a602-603 and inter-cluster a602-610 coherences (figure 4).

Zoom In Zoom Out Reset image size

At weak coupling to phonons (λ_j = 10 cm⁻¹) all theories give almost identical kinetics with long-lived oscillations due to coherences between the a602 and a603 sites, and non-oscillatory a602-610 coherence due to nonsecular transfers. Notice that this non-oscillatory component is not present in the a602-610 coherence calculated with the cmR (which is a secular theory), as shown in the right frame of the top row of figure 4. Increasing the phonon coupling to λ_j = 30 cm⁻¹ results in faster equilibration superimposed with damped oscillations. Again, the four theories predict very similar kinetics within an accuracy of some small differences. A further increase in coupling (to λ_j = 100 cm⁻¹) produces even faster equilibration with non-oscillatory decay of the coherences.

Most interesting in the example shown in figure 4 is the dynamics of transfers between the weakly coupled (M_jj' < 15.8 cm⁻¹) a602-603 and a610-611-612 clusters. Weak coupling implies intense nonsecular transfers that maintain the long-lived non-oscillatory coherences (like a602-610 and a603-610 coherences). Remember that the imaginary part of the inter-site coherences determines the coherences between the exciton eigenstates, that produce the exciton wavepacket localized at certain sites [32]. The breaking of these coherences produces exciton delocalization over the sites involved in the eigenstates. In this way the nonsecular transfers keep excitation localization on a602 and a603 sites, thus producing a relatively slow non-oscillatory population of the a610-611-612 sites. But in the secular theory the a602-610 and a603-610 coherences quickly decay (within 100 fs at λ_j = 100 cm⁻¹ and even faster for bigger couplings) thus producing quick delocalization between the a602-603 and a610-611-612 clusters. As a result, the secular cmR predicts that the population of the a610-611-612 sites is too fast (see left frames in figure 4). After we artificially break the coherent mixing between the a602-603 and a610-611-612 clusters (invoking the cmRgF approach with M_c > 15.8 cm⁻¹) the kinetics of a610-611-612 population becomes more realistic, i.e. not much different from the kinetics given by nonsecular approaches, i.e. sR and HEOM (see the third row in figure 4).

When increasing the coupling to λ_j = 300 cm⁻¹ the secular kinetics is even faster than in the λ_j = 100 cm⁻¹ case, whereas the cmRgF (as well as nonsecular sR and HEOM approaches) give slower kinetics. Comparing the kinetics at various couplings we conclude that they become gradually faster when increasing the phonon coupling (that is equal to the reorganization energy) from λ_j = 10 to 30 and 100 cm⁻¹, but further increasing it to λ_j = 300 cm⁻¹ starts to slow down the kinetics (except the secular ones). This is connected with an increase in the gap between the exciton states predominantly localized at a611-612 and at the red-most a610 site (due to an increase in the reorganization-induced shift of the lowest state localized at a610). This is mirrored by a pronounced increase in the steady-state population of the a610 site at λ_j = 300 cm⁻¹ (that can be seen in the cmR kinetics in figure 4). The observed phenomenon reflects a general rule of the electronic energy transfers: more phonons produce faster transfers, but too many phonons make excitations less mobile due to reorganization-induced red-shifting, dynamic localization and self-trapping (polaron effects). Notice that the first factor is included both in the Redfield and HEOM, the second one (dynamic localization) is present only in HEOM, and the third factor (polaron formation) is not included in the HEOM (see explanation of these terms in appendix A).

It is also interesting to note that the equilibration rates are very similar in cmRgF, HEOM and sR when reorganization is not too big (λ_j < 100 cm⁻¹). But for a realistic reorganization value of λ_j = 300 cm⁻¹ the cmRgF approach gives faster kinetics (probably due to the instantaneous character of the bath relaxation caused by the Markov approximation both in the Redfield and Förster frameworks [51], whereas sR predicts slower kinetics (most probably due to its restriction to the weak coupling limit, implying a one-phonon character of exciton relaxation). Notice also that in cmRgF kinetics there is a very fast component in the a602 and a603 kinetics (<50 fs, with an amplitude of about 20%) connected with a quick decay of the a602-603 coherences (a feature that is not present in sR or HEOM, where coherences are more long-lived due to nonsecular transfers).

Alternative versions of cmRgF, where we break more exciton couplings, do not improve the results. Thus, we tried cmRgF with M_c = 25 cm⁻¹, corresponding to breaking the coherent mixing between the a610 and a611-612 sites. This gives a localized (and more red-shifted) a610 site, which becomes more populated (compared with the a611-612 sites) in contrast with the true (HEOM) kinetics (data not shown). The cmRgF with M_c = 25 cm⁻¹ is therefore not as good as with M_c = 16 cm⁻¹. We thus conclude that breaking the couplings between 16 and 25 cm⁻¹ is not justified.

The equilibration dynamics within the whole stromal side layer (containing five Chls a and three Chls b) can be satisfactorily modeled with the same critical coupling, i.e. M_c = 16 cm⁻¹, as shown in figure S1, which is available at stacks.iop.org/JPB/50/124003/mmedia (see supplementary data).

4.2. Chls b → Chls a transfers in the luminal side (N = 6 model)

Now we switch to excitation dynamics at the luminal side. Here the a613-614 dimer is very weakly (2–3 cm⁻¹) coupled to the a604-b605-606-607 compartment. Due to weak coupling in combination with a big energy gap between the a613-614 and a604-b605-606-607 sites (the nearest in energy a604 and a614 sites are separated by 204 cm⁻¹ (see table 1)) we consider them as different clusters. We should also consider b605 separately from the a604-b606-607 sites (as will be justified below), meaning that the luminal-side pigments are split into three clusters, i.e. a613-614, b605 and a604-b606-607, corresponding to cmRgF approach with M_c = 30 cm⁻¹.

Initial coherent excitation of the b606-607 dimer is followed by the population of the b605 and a604 bottleneck sites with subsequent very slow (30 ps) population of the a613-614 dimer. Kinetics are the same at very weak coupling to phonons (λ_j = 10 cm⁻¹), but a slight increase in coupling to λ_j = 30 cm⁻¹ reveals a pronounced difference between secular (cmR) and nonsecular approaches (figure 5). Thus, the b605 population is much faster in cmR. To explain this phenomenon we note that the b605 site is in resonance with the lowest exciton level of the b606-607 dimer (predominantly localized at b607). In such a situation the transfers between the b605 and b607 sites are very sensitive to inter-site coherences. Thus, a long-lived imaginary part of the b605-b607 coherence will maintain the localization of the excitation at the initially excited b607 site and prevent quick population of the b605 site (as revealed by HEOM and sR dynamics). But in the secular approximation (in cmR) a quick decay of the b605-b607 coherence (due to the absence of nonsecular coherence transfers) produces fast delocalization with equilibration of the b605 and b607 site populations. Such a difference between the cmR and HEOM/sR kinetics becomes even more dramatic at a realistic coupling value (λ_j = 300 cm⁻¹), when the b605 and b607 populations are almost equal after 100 fs (within the cmR picture) in contrast with the nonsecular dynamics (HEOM, sR) giving slow (ps) equilibration between b605 and b607. This artifact of the cmR theory can be corrected by switching to the cmRgF approach, i.e. by considering the b605 and a604-b606-607 as separate clusters. Breaking the mixing of b605 with the b606-607 sites in cmRgF (meaning that M_c = 30 cm⁻¹) gives a result that is very close to the exact solution (HEOM).

Zoom In Zoom Out Reset image size

4.3. Equilibration within the whole subunit (N = 14 model)

Finally, we calculate the dynamics within the whole monomeric subunit of LHCII containing 14 Chls. The initial conditions are the same as those that we used in calculating the kinetics within the stromal- and luminal-side layers, i.e. we suppose broad-band excitation of all the Chls b (except red-shifted b605), including initial coherences within the b606-607 and b608-609 dimers.

Couplings between the stromal-side and luminal-side pigments are quite small (M_jj' < 12 cm⁻¹, see table 1), so the transfers between the exciton states from different layers should be calculated using the generalized Förster theory, because the modified Redfield theory can give strongly overestimated rates in this case. To check this we compare the population dynamics calculated with the cmR theory (see figure 6) with the cmRgF approach with M_c = 12 cm⁻¹, considering the stromal-side and luminal-side layers as two separate clusters, while retaining all the exciton interactions within these clusters, and, therefore, treating the dynamics within each layer with the modified Redfield theory (in figure 6 we denote this two-cluster model as cmRgF-2). In the cmRgF picture we have slow transfers from the a604-b605-606-607 sites to the a613-614 dimer on the luminal side. The a613-614 sites are also getting populated due to transfers from the stromal-side a602-603 sites. When we switch to the cmR picture and include exciton mixing between the layers, we find much faster transfers, including (i) a quick (sub-ps) rise of the a613-614 population due to transfers from a602-603 (see their reduced populations on the sub-ps scale), and (ii) faster decay of the a604-b605-606-607 populations due to transfers to the stromal-side pigments. Obviously, these fast transfers are artifacts of the modified Redfield theory, appearing when considering migration between weakly coupled and quasi-resonant sites [5]. These artifacts are connected with the secular character of the modified Redfield model [32].

Zoom In Zoom Out Reset image size

The two-cluster cmRgF-2 model gives a more correct description of the transfers between the two layers, but in order to obtain a more realistic picture of the dynamics within stromal- and luminal-side layers we have to split each of these layers into several separate clusters. As suggested in previous sections, realistic dynamics can be obtained by considering three clusters in the stromal side, i.e. a602-603, a610-611-612 and b601-608-609, and three clusters in the luminal side, i.e. a613-614, a604-b606-607 and b605. Altogether, we have a six-cluster model (denoted as cmRgF-6 in figure 6).

The main differences between the cmRgF-2 and cmRgF-6 models are related to (i) slower transfers from a602-603 to a610-611-612, and (ii) slower population of b605 in the cmRgF-6 picture (the same features have been observed in cmRgF models for the separate layers, as shown in figures 4, 5 and S1). Notice that the cmRgF-6 model is closer to the exact HEOM solution. The most remarkable differences between these two models become apparent as (i) a bigger population of the b607 site at sub-ps delays and (ii) non-monotonous kinetics of the a610 population (with a fast component on the 0.5 ps time scale) in the HEOM picture. Notice that similar features are also present in the sR model (at least within the 0.5 ps time window), so that they can be assigned to nonsecular transfers. Another factor that can influence the HEOM kinetics is the non-Markovian character of the off-diagonal phonon-induced fluctuations, that determine the (time-dependent) relaxation dynamics and a shifting of the exciton eigenstates (features not included in the standard and modified Redfield pictures). This may be the reason for the deviation of the a610 kinetics in HEOM from both sR and cmRgF models. In spite of these discrepancies the cmRgF-6 model looks very realistic. Notice that the sR approach gives the same character of the dynamics but the transfers corresponding to big energy gaps are too slow (most probably due to the one-phonon character of relaxation in the sR restricted to weak exciton–phonon coupling). Thus, the b → a transfers on the stromal side (corresponding to the biggest gaps) are most difficult to explain with the sR theory (see the dynamics of b601-608-609 depopulation in figures 6 and S1).

Notice that the damping constant (γ_j = 500 cm⁻¹) used in all our examples (figures 4–6) corresponds to fast dephasing (about 10 fs), meaning that the dynamics is close to the Markovian limit. In this respect, it is interesting to study the relation between the Markovian cmRgF model and the exact HEOM solution also for other damping constants γ_j, especially for the reduced γ_j values corresponding to, essentially, a non-Markovian regime. Modeling of the kinetics with tuning γ_j from 100 to 800 cm⁻¹ has shown (see supplementary data, figures S2 and S3) that in spite of some deviations from HEOM, the cmRgF-6 model is still capable of reproducing the equilibration within the whole monomeric subunit even at low values of γ_j. Notice, in this respect, that in real photosynthetic antennae the spectral density typically contains intense vibrations with the frequencies up to 1700 cm⁻¹ [46, 47]. Therefore, modeling of energy transfers with a low-frequency spectral density (for example, with γ_j = 100 cm⁻¹) is not realistic. On the other hand, the model with γ_j = 500–800 cm⁻¹ gives a better approximation of the spectral density and the spectral lineshapes (as shown in figure 3). So, we believe that our choice of an overdamped Brownian oscillator with γ_j = 500 cm⁻¹ yields a reasonable picture of the dynamics.

Comparing the kinetics calculated with cmRgF and HEOM, one should bear in mind that the observed differences can be determined by two factors. Firstly, the transfer rates (connecting the same pair of eigenstates) can be different when evaluated exactly or perturbatively. Secondly, the properties of eigenstates (like delocalization, phonon-induced shift, homogeneous broadening) are also different if they are treated by different approaches. These properties manifest themselves in linear spectral responses, for instance, in absorption lineshapes. In the modified Redfield model the absorption spectrum is calculated by explicitly including the diagonal exciton–phonon coupling [52], whereas off-diagonal coupling is given by a phenomenological Markovian term determined by the inverse lifetime of the exciton states (the thus-defined term is real and describes a relaxation-induced broadening) [12, 32]. In recent papers more versatile theories of absorption lineshapes have been developed [53–56], containing nonsecular and non-Markovian generalizations of the modified Redfield lineshapes. In section S3 of the supplementary data we have studied the absorption spectra of LHCII, calculated with modified Redfield theory, complex time-dependent Redfield (ctR) theory [54], and HEOM. We conclude that both ctR and HEOM give significant and nonuniform red-shifting of the exciton transitions (as compared to the modified Redfield theory). As a result, the relative energies of the exciton states are different in the cmRgF and HEOM models. This factor is responsible for the differences in kinetics shown in figure 6 (together with other factors connected with different treatment of the phonon-induced transfers between the states in cmRgF and HEOM).

5.1. The choice of critical coupling

5.2. Comparing with other models of LHCII

5.3. Comparing with models of other photosynthetic complexes

We demonstrate the applicability of the Redfield–Förster approach for the modeling of energy transfers in photosynthetic light-harvesting complexes. An adequate determination of strongly coupled compartments of the whole complex allows a quantitatively correct description of the kinetics, free from the artifacts of the pure modified Redfield theory (mostly connected with its secular character). Neglecting nonsecular transfers can give a significant overestimation of the transfer rates between weakly coupled and isoenergetic sites. To avoid this artifact we have to artificially break the exciton mixing of these weakly coupled pigments and invoke the generalized Förster theory to calculate the transfer between them. Such an approach results in a more realistic picture that is not much different from the exact solution. A critical cutoff M_c, indicating which exciton couplings should be broken, is dependent on the energy gap between the corresponding sites. Thus, our modeling of the monomeric LHCII subunit demonstrates that a formal comparison of the exciton couplings M_jj' with some fixed M_c value is not sufficient to define a suitable (optimal) compartmentalization scheme.

VIN was supported by the Russian Foundation for Basic Research (Grant No. 15-04-02136); RvG was supported by the VU University Amsterdam, the Laserlab-Europe Consortium, the TOP grant (700.58.305) from the Foundation of Chemical Sciences part of NWO, the advanced investigator grant (267333, PHOTPROT) from the European Research Council, and by the EU FP7 project PAPETS (GA 323901), and by his Academy Professor grant from the Netherlands Royal Academy of Sciences (KNAW).

We consider the molecular aggregate made of N coupled pigments. The one-exciton Hamiltonian in the site representation is:

$$\begin{eqnarray}&&{H}=\displaystyle \sum _{{j}=1}^{{N}}{\omega }_{{j}0}| {j}\rangle \langle {j}| +\displaystyle \sum _{{j}\ne {j}^{\prime} }{{M}}_{{j}{j}^{\prime} }| {j}\rangle \langle {j}^{\prime} | \end{eqnarray} \tag{ A.1 }$$

where ${\omega }_{{j}0}={{\omega }_{{j}0}}^{0}+{\lambda }_{{j}}$ is the energy of the Frank–Condon transition (including the zero-phonon energy ${{\omega }_{{j}0}}^{0}$ and reorganization energy λ_j), M_jj' is the interaction energy between the jth and j'th sites (the ω_j0 and M_jj' values are listed in table 1). A coupling of the jth site to fast nuclear modes is described by the spectral density in the form of an overdamped Brownian oscillator (with reorganization energy λ_j and damping constant γ_j):

$$\begin{eqnarray}&&{{C}}_{{j}}(\omega )=2{\lambda }_{{j}}\displaystyle \frac{\omega {\gamma }_{{j}}}{{\omega }^{2}+{\gamma }_{{j}}^{2}}.\end{eqnarray} \tag{ A.2 }$$

The scaled hierarchical equation of motion (HEOM) for the reduced density operator is [49]:

$$\begin{eqnarray}{\dot{\sigma }}^{{n}} & = & -\left({\rm{i}}\,{{H}}^{\times }+{G}+\displaystyle \sum _{{j}=1}^{{N}}{{S}}_{{j}}{{V}}_{{j}}^{\times }{{V}}_{{j}}^{\times }\right){\sigma }^{{n}}-{\rm{i}}\,\displaystyle \sum _{{j}=1}^{{N}}\displaystyle \sum _{{k}=0}^{{K}}({{L}}_{{jk}}{{V}}_{{j}}^{\times }{\sigma }^{{{n}}_{{jk}}^{+}}\\ & & +{{G}}_{{jk}}{{V}}_{{j}}{\sigma }^{{{n}}_{{jk}}^{-}}-{\sigma }^{{{n}}_{{jk}}^{-}}{{V}}_{{j}}{{G}}_{{jk}}^{* })\\ {G} & = & \displaystyle \sum _{{j}=1}^{{N}}\displaystyle \sum _{{k}=0}^{{K}}{{n}}_{{jk}}{\nu }_{{jk}};\,{{S}}_{{j}}=\displaystyle \sum _{{\rm{k}}={K}+1}^{\infty }\displaystyle \frac{{{c}}_{{jk}}}{{\nu }_{{jk}}};\\ {{L}}_{{jk}} & = & \sqrt{({{n}}_{{jk}}+1)\,| {{c}}_{{jk}}| };\,{{G}}_{{jk}}=\sqrt{\displaystyle \frac{{{n}}_{{jk}}}{| {{c}}_{{jk}}| }}\,{{c}}_{{jk}};\\ {{c}}_{{j}0} & = & {\lambda }_{{j}}{\gamma }_{{j}}\left(-{\rm{i}}+\,\cot \,\displaystyle \frac{{\gamma }_{{j}}}{2\theta }\right);\,{\nu }_{{j}0}={\gamma }_{{j}};\\ {{c}}_{{jk}} & = & \displaystyle \frac{4\theta {\lambda }_{{j}}{\gamma }_{{j}}{\nu }_{{jk}}}{{\nu }_{{jk}}^{2}-{\gamma }_{{j}}^{2}};\,{\nu }_{{jk}}=2\pi {k}\theta ;\,({k}\geqslant 1)\end{eqnarray} \tag{ A.3 }$$

where σⁿ denote auxiliary operators dependent on the set of integers describing the state of the phonon bath of the sites from j = 1 to N; n = {n₁₀,..n_1K,..n_jk,..n_N0,..n_NK}; ${{{n}}_{{jk}}}^{\pm }\,=\{\ldots \,{{n}}_{{jk}}\pm 1\ldots \};$ the index k numbers the low-temperature correction terms from k = 0 to K; θ = k_BT, where k_B is the Boltzmann constant, T is the temperature. The auxiliary operators σⁿ at n = {0,0,...0} are equal to the reduced density operator ρ. Here V^×σ denotes a commutator, i.e. V^×σ = Vσ − σV, where V_j = ∣j〉〈j∣. Notice that equation (A.3) implies that the phonons associated with different sites are uncorrelated. The dynamics of the one-exciton populations σ_ii and coherences σ_ij (in the site representation) is given by:

$$\begin{eqnarray}{\ddot{\sigma }}_{{ij}}^{{n}} & = & (-{\rm{i}}{\omega }_{{ij}}-(1-{\delta }_{{ij}})({{S}}_{{i}}+{{S}}_{{j}})-{G}){\sigma }_{{ij}}^{{n}}\\ & & -{\rm{i}}\displaystyle \sum _{{s}=1}^{{N}}({{M}}_{{is}}{\sigma }_{{sj}}^{{n}}-{\sigma }_{{is}}^{{n}}{{M}}_{{sj}})\\ & & -{\rm{i}}\,\displaystyle \sum _{{k}=0}^{{K}}({{L}}_{{ik}}{\sigma }_{{ij}}^{{{n}}_{{ik}}^{+}}-{\sigma }_{{ij}}^{{{n}}_{{jk}}^{+}}{{L}}_{{jk}}+{{G}}_{{ik}}{\sigma }_{{ij}}^{{{n}}_{{ik}}^{-}}-{\sigma }_{{ij}}^{{{n}}_{{jk}}^{-}}{{G}}_{{jk}}^{* })\end{eqnarray} \tag{ A.4 }$$

with ω_ij = ω_i0 − ω_j0., and with some initial conditions for σ_ij describing the state of the system after impulsive excitation. For example, excitation of the jth site corresponds to ${\sigma }_{{ii}}^{{n}}=1$ for n = {0,0,...0} and ${\sigma }_{{ii}}^{{n}}=0$ for other n sets. If the sum of n_jk exceeds N_c, then the terms with G, G_ik, and L_ik in equations (A.3) and (A.4) are set to zero [62]. The depth of the hierarchy N_c and the number of the Matsubara frequencies K should be chosen in order to get converging results.

We remind the reader that the Redfield approach describes the dynamics in the eigenstate (exciton) basis, where the exciton states are obtained by diagonalization of Hamiltonian (A.1) without taking into account the coupling of the diabatic electronic states ∣j〉 to the phonons. In other words, a delocalization of the exciton states does not depend on phonons (in this way, even small coupling M_jj' will produce delocalization between the jth and j'th sites). In contrast, in equation (A.3) a mixing of diabatic states is dependent on the state of the phonon bath coupled to each of the mixed sites (given by the {n₀,..n_j,..} set). If the jth and j'th sites are represented as two potential surfaces with different displacements along some nuclear coordinate, then excitation near the crossing point of the two potentials will create a delocalized excitation, but subsequent reorganization will produce relaxation of the excited-state wavepacket to the bottom of the mixed potential that is predominantly localized at one of the two states. This phenomenon can be denoted as 'dynamic localization'. Notice that in equations (A.2)–(A.3) the excited states are coupled to a continuum of phonon modes (Brownian oscillator) displaced along a single effective nuclear coordinate. In such a model dynamic localization can occur only if the states are characterized by different absolute values of the displacement (i.e. if they have different λ_j values).

In a more realistic treatment, the dynamic localization is accounted for by switching from the one-dimensional case to a multidimensional configuration space, where the coupling of individual (diabatic) excited states to a single vibrational mode are represented as potential surfaces shifted along the x_j direction (with j = 1 to N). Obviously this is impossible to do if the spectral density has the form of an overdamped Brownian oscillator, but it can be done for a limited number of underdamped vibrations. The simplest example is a dimer coupled to a single vibration, characterized by two excited states displaced along the x and y coordinates. In such a system two weakly coupled states with the same absolute values of the displacement (associated with the ground to one-exciton transition) become localized. Impulsive excitation creates a delocalized excitation located near the crossing area (corresponding to the x + y direction) with subsequent dynamic localization due to motion in the adiabatic, or 'anticorrelated' x − y direction. This model was first considered by Förster [63], who demonstrated that the rate of energy transfer between the jth and j'th sites in a weak coupling limit is given by the square of interaction energy M_jj' and by the overlap of the vibrational wavefunctions. Later this model was revisited and generalized to arbitrary strong coupling in order to describe the exciton-vibrational coherences observed in two-dimensional electronic spectroscopy [64, 65]. The dynamics in such a model can be described by Redfield theory on the exciton-vibrational basis [57] or by HEOM with appropriate choice of the displacements along nuclear coordinates x_j.

Notice that the phenomenon of dynamic localization described above is connected with reorganization dynamics within the coupled states with some fixed displacements. Generally, a change of equilibrium positions upon excitation (i.e. displacements of the excited-state potentials) can be dependent on excitation density. This means that nonuniform excitation of some molecular cluster can create a nonuniform distribution of the displacements of the diabatic potential surfaces, thus producing a local deformation of the lattice (pigment-protein matrix). This gives rise to an even more localized exciton wavefunction. This phenomenon, first described in Davydov's early work [66, 67], is usually denoted as 'self-trapping' or 'polaron' formation [68]. This additional mechanism of localization is not included in the Redfield and HEOM theories, and requires more sophisticated approaches. Thus, the polaron wavefunction can be obtained only variationally. An example of such a calculation for a photosynthetic antenna is given in [68, 69].

The coherent modified Redfield (cmR) approach [30, 31] is a non-Markovian secular master equation describing the dynamics of the one-exciton populations ρ_ss and the decay of the coherences ρ_ss' between the exciton states. Nonsecular terms, i.e. transfers between coherences and transfers between populations and coherences are not included. As in the original modified Redfield model [52] the diagonal exciton–phonon coupling is included explicitly, whereas the off-diagonal phonon-induced fluctuations (inducing relaxation between the exciton states) are treated perturbatively. The exciton wavefunctions are calculated by the diagonalization of Hamiltonian (A.1), i.e. they are independent of the phonon coordinates, so that the dynamic localization is not included. In the Markovian limit (used in our study) the s' → s population transfers are given by the tensor R_sss's' that is the same as in original modified Redfield theory [5, 52]:

$$\begin{eqnarray}{{R}}_{{\text{ss}}{s}^{\prime} {s}^{\prime} } & = & 2{\rm{Re}}\displaystyle {\int }_{0}^{\infty }{\rm{d}}{t}{{A}}_{{s}}({t})\,{{F}}_{{s}^{\prime} }^{* }({t})\,{{V}}_{{s}{s}^{\prime} }({t})\\ {{A}}_{{s}}({t}) & = & {\rm{\exp }}\{{\rm{i}}(\omega -{\omega }_{{s}}){t}-{{g}}_{{ssss}}({t})\}\\ {{F}}_{{s}}({t}) & = & {\rm{\exp }}\{{\rm{i}}(\omega -{\omega }_{{s}}){t}+2{\rm{i}}{\lambda }_{{ssss}}{t}-{{g}}_{{ssss}}^{\ast }({t})\}\\ {{V}}_{s{s}^{\prime} }({t}) & = & {\rm{\exp }}(2{{g}}_{{s}^{\prime} {s}^{\prime} {\text{ss}}}({t})+2{\rm{i}}{\lambda }_{{s}^{\prime} {s}^{\prime} {\text{ss}}}{t})\\ & & \times [{\ddot{{g}}}_{{s}{s}^{\prime} {s}^{\prime} {s}}({t})-\{{\mathop{{g}}\limits^{}}_{{s}^{\prime} {s}{s}^{\prime} {s}^{\prime} }({t})-{\mathop{{g}}\limits^{}}_{{s}^{\prime} {\text{sss}}}({t})+2{\rm{i}}{\lambda }_{{s}^{\prime} {s}{s}^{\prime} {s}^{\prime} }\}\\ & & \times \{{\mathop{{g}}\limits^{}}_{{s}^{\prime} {s}^{\prime} {s}{s}^{\prime} }({t})-{\mathop{{g}}\limits^{}}_{{\text{sss}}{s}^{\prime} }({t})+2{\rm{i}}{\lambda }_{{s}^{\prime} {s}^{\prime} {s}{s}^{\prime} }\}]\end{eqnarray} \tag{ B.1 }$$

where F(t) and A(t) are line-shape functions corresponding to the fluorescence of the donor state and the absorption of the acceptor, respectively, while V describes the interaction between donor and acceptor, and ω_s is the energy of the sth exciton state. Equation (B.1) is valid for arbitrary delocalization of the donor and acceptor states.

Other quantities are related to the exciton–phonon spectral density in the site representation C_j(ω):

$$\begin{eqnarray}{{\rm{g}}}_{{ss}^{\prime} {s}^{\prime\prime} {{s}}^{\prime\prime ^{\prime} }}({t}) & = & \displaystyle \sum _{{j}}{{c}}_{{j}}^{{s}}{{c}}_{{j}}^{{s}^{\prime} }{{c}}_{{j}}^{{s}^{\prime\prime} }{{c}}_{{j}}^{s^{\prime\prime} \,^{\prime} }{{\rm{g}}}_{{j}}({t})\\ {\lambda }_{{ss}^{\prime} {s}^{\prime\prime} {{s}}^{\prime\prime ^{\prime} }}({t}) & = & \displaystyle \sum _{{j}}{{c}}_{{j}}^{{s}}{{c}}_{{j}}^{{s}^{\prime} }{{c}}_{{j}}^{{s}^{\prime\prime} }{{c}}_{{j}}^{s^{\prime\prime} \,^{\prime} }{\lambda }_{{j}}\\ {g}_{{j}}({t}) & = & -\displaystyle {\int }_{-\infty }^{\infty }\displaystyle \frac{{\rm{d}}\omega }{2\pi {\omega }^{2}}{{C}}_{{j}}(\omega )\left[(\cos \,\omega {t}-1)\coth \,\tfrac{\omega }{2\theta }\right.\\ & & -{\rm{i}}\,({\rm{\sin }}\,\omega {t}-\omega {t})\phantom{\rule[-0000]{}{2.4ex}}]\\ {\lambda }_{{j}} & = & \displaystyle {\int }_{-\infty }^{\infty }\displaystyle \frac{{\rm{d}}\omega }{2\pi \omega }{{C}}_{{j}}(\omega )\end{eqnarray} \tag{ B.2 }$$

where g_j(t) and λ_j are the line-broadening function and reorganization energy of the jth diabatic state. Amplitudes ${{c}}_{{j}}^{{s}}$ show participation of the jth site to the sth exciton state. If the donor and acceptor states are localized at the jth and ith sites (i.e. ${{c}}_{{j}}^{{s}^{\prime} }=1$ and ${{c}}_{{i}}^{{s}}=1$) then the transfer between them is given by the Förster formula, which can be obtained from (B.1) by replacing the interaction term by

$$\begin{eqnarray}&&{{V}}_{s{s}^{\prime} }={| {{M}}_{{ij}}| }^{2}\end{eqnarray} \tag{ B.3 }$$

where M_ij is the interaction energy corresponding to a weak coupling between the localized sites i and j. Switching to the Fourier transforms of F(t) and A(t) we can rewrite the integral in a form of donor–acceptor spectral overlap [63]. The standard Förster formula can be generalized to the case of energy transfer between two weakly connected clusters [4, 5]. The rate of energy transfer from the s'th exciton state of one cluster to the sth state of the other cluster is given by (B.1) with the interaction term:

$$\begin{eqnarray}&&{{V}}_{s{s}^{\prime} }={\left|\displaystyle \sum _{{i},{j}}{{c}}_{{i}}^{{s}}{{M}}_{{i}{j}}{{c}}_{{j}}^{{s}^{\prime} }\right|}^{2}\end{eqnarray} \tag{ B.4 }$$

where i and j designate molecules belonging to different clusters. In this generalized Förster formula, the donor and acceptor states s' and s can have an arbitrary degree of delocalization (corresponding to arbitrarily strong excitonic interactions within each cluster), but the inter-cluster interactions M_ij are supposed to be weak.

In the combined Redfield–Förster approach the population transfers and the decay of the coherences within strongly coupled clusters (with intra-cluster interactions M_ij > M_c) is calculated with the coherent modified Redfield theory, whereas transfers between these clusters (with inter-cluster couplings M_ij < M_c) are modeled by the generalized Förster theory. The decay of the coherences within strongly coupled clusters is given by [31]:

$$\begin{eqnarray}{{R}}_{{s}{s}^{\prime} {s}{s}^{\prime} } & = & \tfrac{1}{2}{{R}}_{{ssss}}+\tfrac{1}{2}{{R}}_{{s}^{\prime} {s}^{\prime} {s}^{\prime} }+{{\rm{\Gamma }}}_{{s}{s}^{\prime} };\\ {{\rm{\Gamma }}}_{{s}{s}^{\prime} }({t}) & = & \displaystyle \sum _{{j}}{[{({{c}}_{{j}}^{{s}})}^{2}-{({{c}}_{{j}}^{{s}^{\prime} })}^{2}]}^{2}({\ddot{{g}}}_{{j}}({t})+{\rm{i}}{\lambda }_{{j}})\end{eqnarray} \tag{ B.5 }$$

where R_ssss is the inverse lifetime of the sth exciton state (determined by the sum of the s → s' relaxation rates R_s's'ss), Γ_ss' is pure dephasing. In the Markovian approximation the time-independent dephasing term is (bearing in mind that ${\rm{Im}}\{\dot{{g}}(\infty )\}=-\lambda$):

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{s}{s}^{\prime} }(\infty )=\displaystyle \sum _{{j}}{[{({{c}}_{{j}}^{{s}})}^{2}-{({{c}}_{{j}}^{{s}^{\prime} })}^{2}]}^{2}{\rm{Re}}\{{\ddot{{g}}}_{{j}}(\infty )\}.\end{eqnarray} \tag{ B.6 }$$