Normal mode analysis of the spectral density of the Fenna-Matthews-Olson light-harvesting protein: how the protein dissipates the excess energy of excitons.

We report a method for the structure-based calculation of the spectral density of the pigment-protein coupling in light-harvesting complexes that combines normal-mode analysis with the charge density coupling (CDC) and transition charge from electrostatic potential (TrEsp) methods for the computation of site energies and excitonic couplings, respectively. The method is applied to the Fenna-Matthews-Olson (FMO) protein in order to investigate the influence of the different parts of the spectral density as well as correlations among these contributions on the energy transfer dynamics and on the temperature-dependent decay of coherences. The fluctuations and correlations in excitonic couplings as well as the correlations between coupling and site energy fluctuations are found to be 1 order of magnitude smaller in amplitude than the site energy fluctuations. Despite considerable amplitudes of that part of the spectral density which contains correlations in site energy fluctuations, the effect of these correlations on the exciton population dynamics and dephasing of coherences is negligible. The inhomogeneous charge distribution of the protein, which causes variations in local pigment-protein coupling constants of the normal modes, is responsible for this effect. It is seen thereby that the same building principle that is used by nature to create an excitation energy funnel in the FMO protein also allows for efficient dissipation of the excitons' excess energy.

Introduction

The important task in photosynthesis of collecting sunlight and transferring the excitation energy to a reaction center (RC), where it is converted into chemical energy, is performed by specialized pigment–protein complexes (PPCs) termed antennae or light-harvesting complexes. In general, the excitation energy is guided to the RC through energy funnels, which the proteins create, holding the pigments at optimal positions and varying their electronic properties. At least two different types of funnels exist.

By adjusting interpigment distances, an excitation energy funnel is created in the photosynthetic apparatus of purple bacteria.[1−3] Another type of excitation energy funnel is realized in the Fenna–Matthews–Olson (FMO) protein of green sulfur bacteria:[4−6] The local transition energies of the bacteriochlorophyll a (BChla) pigments in their binding sites, the site energies, are varied by the pigment–protein interaction in a way that pigments facing the RC are red-shifted with respect to those facing the outer antenna system. Evidence for this site energy funnel was obtained from fits of optical spectra[7,8] as well as direct structure-based quantum chemical/electrostatic computations.[9,10]

Besides this tuning of excitation energies, the protein has to dissipate the excess energy of the excitons when they relax downward to the excited states of the RC. This dissipation is achieved primarily by protein vibrations that absorb the excess energy and distribute it over many degrees of freedom. The coupling of pigment excited states to the protein vibrations is characterized by the spectral density J(ω).[3,8,11] Fortunately, information about J(ω) can be obtained from optical line-narrowing experiments.[11−13] If only a single pigment is present as in the B777 complex[14,15] derived from the light-harvesting complex 1 (LH1) of purple bacteria, the standard theory of a two-level system coupled linearly to a bath of harmonic oscillators[16,17] can be used to extract J(ω) from the experimental data.[11] If more than one pigment is present, the standard theory can still be applied, provided that the lowest exciton state is strongly localized on a single pigment.[12,13] If the lowest exciton state is delocalized, more advanced approximations have to be used to describe the line-narrowing spectra.[11] These approximations concern the theory of exciton–vibrational dynamics and the different contributions to the spectral density.

Concerning theory, the problem of equal excitonic and exciton–vibrational coupling strengths arises. In principle, this problem can be solved numerically by using one of the non-perturbative approaches, like the hierarchical equation of motion (HEOM) approach,[18−21] the density matrix renormalization/polynomial transformation approach,[22,23] or path integral techniques.[24,25] The numerical effort, however, is considerable. So far, these techniques have not been applied to analyze line-narrowing spectra of PPCs. Alternatively, the Hamiltonian of the PPC has been transformed into the basis of delocalized exciton states, with the hope to find a small parameter that can be used for perturbation theory.[11,26−32]

In this context, two non-Markovian density matrix theories, the time-local chronological ordering prescription (COP) and the time-local partial ordering prescription (POP) theory, were tested against experimental fluorescence line-narrowing (FLN) data of the dimeric B820 subunit of LH1.[11] In these calculations, the spectral density of site energy fluctuations derived from experiments on the B777 complex[15] was used. POP was found to provide more accurate results for the FLN spectrum of the B820 dimer.[11] Additionally, it was shown that COP breaks down at high temperatures in simulations of linear absorbance spectra of the water-soluble chlorophyll binding protein (WSCP), whereas POP again was in agreement with experiment.[33] The exciton–vibrational coupling constants obtained in the present work from a normal-mode analysis (NMA) provide a microscopic explanation of these results.

The spectral density of a PPC contains several contributions resulting from the fluctuation of site energies, the fluctuation of excitonic couplings as well as correlations among and between these two types of fluctuations. These contributions influence exciton relaxation and the dephasing of coherences in different ways. Therefore, it is of considerable interest to find a way to estimate the importance of the various contributions. The main goal of the present paper is to offer a tool for this purpose. In the absence of methods for a quantitative evaluation of the different parts of the spectral density, it is often presumed that the spectral density is dominated by site energy fluctuations and that the other contributions can be neglected. To model correlations among the former, an empirical correlation radius Rc has been introduced by assuming an exponential dependence of cross correlations on interpigment distance.[11,29,34] In recent years, the interest in such correlated fluctuations and their possible functional role has increased significantly.[35−39]

Rebentrost et al.[35] and Sarovar et al.[36] studied energy transfer and the trapping efficiency in the FMO protein as a function of Rc. Both studies found a slower transport and a higher degree of coherences in the system if Rc is large, in agreement with earlier results.[33] Sarovar et al.[36] raised the intriguing question of why fluctuations should be correlated, if these correlations reduce the efficiency of light harvesting. As will be shown here, the spectral density is more complicated than assumed in these phenomenological models. Abramavicius and Mukamel[37] investigated the signature of correlated site energy fluctuations in 2D photon echo spectra and confirmed their retarding influence on energy transfer. They identified the cross peak regions of the 2D spectra as particularly sensitive to the correlations.

Other groups additionally considered the case of anticorrelated site energy fluctuations, fluctuations of excitonic couplings and their correlation, as well as correlations between site and coupling fluctuations.[40−45] Whereas in the simple models, correlations prolong the quantum beating and the energy transfer times, the more general models show that, in principle, it is possible to simultaneously increase the lifetime of coherences, and to decrease the energy transfer times.[44] Thus, there is a clear need for a microscopically calculated spectral density to clarify the possible role of correlations.

Pioneering work in this direction has been performed by the Schulten group,[46] who combined molecular dynamics (MD) simulations of the LH2 complex of purple bacteria with quantum chemical (QC) calculations of the pigment’s transition energies to compute the autocorrelation functions of site energy fluctuations. The protein was included as a classical background charge distribution in the solution of the electronic Schrödinger equation. The spectral density was then obtained from the Fourier transform of the autocorrelation function. This approach has been modified by using instead of the original ab initio electronic structure method either a semiempirical method[47−51] or density functional theory (DFT).[52] The simulations are computationally demanding, since the QC site energy calculations have to be performed every 1–2 fs along the classical MD trajectories. Besides the spectral density of LH2,[46−48] those of the FMO protein[49,50,52] and the RC of purple bacteria[51] have been investigated.

Despite considerable differences in the shape of the reported spectral densities of site energy fluctuations, all the QC/MD approaches discussed above (except for the ones by Jing et al.[51] and Kim et al.,[53] see below) have in common that they severely overestimate the contribution from intramolecular pigment vibrations (occurring at high frequencies) when compared with the spectral density obtained from line-narrowing spectra.[12,54,55] Furthermore, significant deviations between the QC/MD results and the spectral densities derived from experiment are also observed in the low-frequency region (ω < 500 cm1) that is particularly important for light-harvesting. As pointed out by Aspuru-Guzik and co-workers,[52] one possible reason for the overestimation of the high-frequency part could be the classical treatment of these modes because of the quantum mechanic freezing-out effect, which suppresses the motion of high-frequency degrees of freedom. Alternatively, Kim et al.[53] noted that the missing zero point energy of the vibrations in a classical description may be responsible for the spurious appearance of high-frequency components in the spectral density.

Probably the main reason, as reported by Jing et al.,[51] lies in the mismatch between the ground state geometries of the pigments created by the classical force field in the MD run and the optimal geometry obtained from a QC calculation. During the development of computational schemes for site energy shifts,[9,56] we observed a similar behavior: QC calculations based directly on heavy-atom coordinates from a crystal structure (supplemented with hydrogen atom positions from a classical force field optimization) usually resulted in erroneous charge distributions as judged from a too large difference in the permanent dipole moment between ground and excited states (compared to values suggested by Stark effect measurements). This discrepancy was one reason to use a two-step procedure: In the first step, QC calculations are performed on pigments fully geometry-optimized in vacuo, ultimately resulting in atomic partial charges that describe the permanent charge distributions of ground and excited states and the transition densities. In the second step, these partial charges are used in a crystal-structure-based all-atom electrostatic calculation including the whole PPC.[9,10,56,57] A second reason for the use of a two-step procedure is the possible occurrence of electron leakage and overpolarization resulting from an artificial distortion of the molecular electron density due to classical background charges.[58,59] It remains to be explored in detail to what extent these effects contribute to the error of QC/MD approaches in the calculation of spectral densities.

Another subtlety of the calculation of site energies concerns the dependence of the protein charge distribution on the protonation states of titratable amino acid residues. In the framework of an all-atom electrostatic calculation, the protonation states can be handled by solving the linearized Poisson–Boltzmann equation and performing a Monte Carlo average.[60,61] The electrostatic computation of site energies is readily incorporated into this scheme,[9,56] resulting in a method that has been baptized the Poisson–Boltzmann/quantum chemical (PB/QC) approach. In a simplified version thereof, termed the charge density coupling (CDC) method, the protonation states are first determined by the Poisson–Boltzmann methods, then fixed, and, finally, the site energies are computed by approximating the polarization effects through an effective dielectric constant εeff.[10,56,62] In addition, εeff takes into account uncertainties about the exact values of the QC partial charges of the ground and excited states of the pigments, also present in the PB/QC approach. εeff is essentially a fit parameter, and has been found to be εeff = 3 for the FMO protein.[10] The QC/PB and CDC methods have been successfully applied to predict low-energy excited states in various PPCs, including the FMO protein,[9,10,56] the major light-harvesting complex II of plants,[63] and photosystem I.[62]

Jing et al.[51] reported for the first time a combination of the CDC method with MD simulations to infer the spectral density of the pigment–protein coupling. The coupling to intramolecular modes was obtained from a QC-based NMA of the pigments in vacuo. They also tested a variant of the QC/MD method, in which the intramolecular part of the spectral density was calculated separately from a QC/MD simulation of the pigments in vacuo and then subtracted from the QC/MD results of the PPC. The qualitative agreement of the latter difference with the CDC/MD results demonstrates that the mismatch between classical force field and QC calculation may be compensated in this way. Further, it shows that the fluctuations of the Coulomb interaction between pigments and protein are the main contributors to the intermolecular spectral density. However, the trajectory length limiting the accuracy of the spectral density in the low-frequency range remains a problem in MD-based approaches. Jing et al.[51] implicitly took care of this limitation by assuming a temperature-dependent inhomogeneous broadening. In the present work, this problem is circumvented by applying, instead of a MD simulation, a NMA of the PPC. Another comparison between a CDC/MD and a QC/MD approach was recently reported by Kim et al.,[53] who concluded that these two methods provide limiting cases for the experimental spectral density. Whereas with CDC/MD a similar shape of the spectral density was obtained as measured in the experiment, but the calculated amplitude was too small, with QC/MD the amplitude of the spectral density at low frequencies seems to be in better agreement with experimental data, but a comparison is difficult because of the overwhelming high-frequency part that is not observed experimentally (Figures 4 and 5 of their paper[53]). Consequently, the absorbance spectrum obtained with the CDC/MD spectral density fitted the experimental absorbance spectrum much better than the one calculated with QC/MD (Figure 3 of their paper[53]).

The development of 2D electronic spectroscopy and its application to photosynthetic light-harvesting complexes pioneered by the Fleming group has enormously revived the interest in coherent exciton motion in these systems.[64−67] Quantum beats have been observed that last for a couple of hundred femtoseconds even at room temperature[68] and for weakly coupled systems.[69] First signatures of excitonic quantum beats were found in the experimental[70] and theoretical[29] anisotropy of pump–probe spectra of the FMO protein, 15 years ago. An interesting question[71] that awaits to be answered by structure-based simulations is whether the protein environment in light-harvesting systems acts as a “quantum protector” of coherences by a correlated modulation of site energies. So far, neither the QC/MD[48,50,52] nor the CDC/MD[51] methods could find correlations in site energy fluctuations. The QC/MD simulations by Olbrich et al.,[50] however, found signatures of correlations in fluctuations of excitonic couplings and between the fluctuations of couplings and site energies. A quantification of the related spectral densities was, however, not possible so far and will be provided in the present work, together with that of the correlation in site energy fluctuations.

The remaining part of this article is organized in the following way: We first introduce the CDC/TrEsp/NMA method for the calculation of the spectral density and show how this spectral density enters the calculation of optical spectra and the evolution of coherences and populations that are generated by ultrashort optical pulses. Next, the method is applied to the monomeric subunit of the FMO protein, revealing all details of the intermolecular spectral density including the correlations. The relative importance of the different parts of the spectral density is investigated in calculations of the temperature-dependent dephasing of excitonic wavepackets and the population transfer. Finally, the negligible influence of the correlations is explained and implications of the present results for the theory of excitation energy transfer and the mechanism of light harvesting in PPCs are discussed.

Theory and Computational Methods

Spectral Density and Normal Mode Analysis

We consider a standard time-dependent Hamiltonian of the PPCwhere |m⟩ denotes a state where pigment m is excited whereas all other pigments are in their electronic ground state. The diagonal element H(t) is the electronic transition energy between the ground state and the local excited state |m⟩, that is, the site energy of pigment m, whereas the off-diagonal element H(t) with m ≠ n is the excitonic coupling between states |m⟩ and |n⟩. The excitonic couplings are responsible for the transfer of excitation energy and for the delocalization of excited states. The protein environment is described by a set of harmonic oscillators obtained from a NMA, where the dimensionless coordinates Qξ and momenta Pξ are related to creation and annihilation operators of vibrational quanta Cξ† and Cξ, respectively, of normal mode ξ by Qξ = Cξ† + Cξ and Pξ = i(Cξ† – Cξ) and to mass-weighted normal coordinates qξ and momenta pξ by[34,72]The time dependence of H(t) is induced by the nuclear dynamics via the pigment–protein coupling. The matrix elements H(t) are expanded into a Taylor series with respect to small displacements R(t) of atoms j = 1, ..., Natom of the PPC from their equilibrium positions R(0). Including terms up to first order in the displacements giveswhere H(0) and ∇H|0 are the values of H and of its gradient, taken with respect to the three cartesian coordinates of the jth atom, respectively, at the equilibrium position of nuclei in the electronic ground state, i.e. for R = R(0), k = 1...Natom. The mass-weighted normal coordinates qξ(t) are related to the displacements (R(t) – R(0)) by[34]where M is the mass of atom j and A(ξ) contains the contributions of this atom to the eigenvector of normal mode ξ. From eqs 3 and 4, we obtainwhere a dimensionless coupling constant gξ(m,n) was introduced asThis coupling constant will enter the spectral density below (eq 10). In order to evaluate the r.h.s. of eq 6, we need to know how the matrix elements H(t) = H({R(t)}) depend on the nuclear coordinates R. These dependencies are revealed by the TrEsp and CDC methods.

In the TrEsp method,[73] the excitonic coupling H is calculated from the Coulombic coupling of the transition charges of two pigments. The transition charges are determined from a fit of the electrostatic potential (ESP) of the ab initio transition density obtained with time-dependent density functional theory. Because of uncertainties in the magnitude of the QC transition charges and the neglect of local field and screening effects induced by the dielectric environment, the TrEsp method needed to be refined. For this purpose, the Poisson-TrEsp method has been developed,[8,56,74] which takes into account the optical polarizability of the environment represented by the refractive index n. In this method, a Poisson equation is solved for the potential of rescaled transition charges that are placed into molecule-shaped cavities with n2 = 1 inside and n2 = 2 outside. The environmental refractive index has been estimated from comparison of the oscillator strength of protein-bound and solvent-extracted pigments.[75] Rescaling of the transition charges is done by a constant factor such that the correct vacuum value of the transition dipole moment results, as inferred by Knox and Spring[76] from an empty cavity analysis of the oscillator strength of BChla in different solvents. Due to this calibration and the explicit inclusion of the dielectric environment, there is no free parameter left. A comparison of excitonic couplings of the FMO protein obtained with Poisson-TrEsp to those obtained with TrEsp shows that the influence of the optical polarizability on the excitonic couplings may be well approximated by a constant screening/local field correction factor of f = 0.8.[3,8,56] Hence, the excitonic coupling between pigments m and n is given asHere, q((0, 1) is the transition charge of atom k of pigment m (rescaled as described above) and R is the position of this atom, for example. We note in passing that the exponential distance dependence of the screening factor assumed in ref (50) cannot be justified by Poisson-TrEsp.[74]

In order to simplify the calculation of site energies H and to obtain an explicit coordinate dependence, our original PB/QC approach[9,63] has been modified to result in the CDC approach,[10,56,62] as discussed in the Introduction. In the CDC method, the site energy is given aswhere E0 is a constant that is independent of m, q(bg) denotes the ground state partial charge of the background atom i, and q((1, 1) and q((0, 0) are the excited state and ground state partial charges of atom k at the mth pigment, respectively. Note that i runs over all atoms of the PPC, except for those of the macrocycle of pigment m, and k runs only over the latter. The most probable protonation state of the titratable residues of the protein was determined by performing Poisson–Boltzmann type calculations as described in ref (10).

The effective dielectric constant εeff was introduced to describe screening and local field effects of the Coulomb coupling, in the same spirit as f was introduced in the calculation of excitonic couplings in eq 7. Please note, however, that excitonic couplings are influenced only by the electronic polarizability, since the nuclei of the environment have no time to move during an electronic transition, whereas the site energies experience both types of polarizabilities.[9] For the calculation of static-structure site energies in the FMO protein, we estimated εeff = 3 from a comparison of the resulting optical spectra with experimental data.[10] The site energies obtained with CDC[10] are in good agreement with earlier PB/QC values,[9] showing that screening/local field effects can indeed be approximated by using an effective dielectric constant.

Since in the present work the nuclear polarizability is modeled explicitly by the NMA, we are just left with the effect of electronic polarizability that can be described by a factor of f = 0.8, as discussed above for the excitonic couplings. Hence, we use εeff = 1/f = 1.25 in the NMA of site energy fluctuations.

By calculating the gradient in eq 6, using eqs 7 and 8, the coupling constants gξ(m,n) finally are obtained aswhere the equilibrium vectors R(0) are obtained from energy minimization of one FMO monomer taken from the crystal structure with hydrogen atoms addded by molecular modeling. The A( are obtained from the eigenvectors of the NMA. The dimensionless coupling constants gξ(m,n) are the weighting factors for the density of vibrational states N(ω) = ∑ξ δ(ω – ωξ), revealing the spectral density J(ω) of the exciton vibrational couplingHere, J(ω) describes the fluctuation of the site energy of pigment m and J(ω) that of the excitonic coupling between pigments m and n. The correlation of site energy fluctuations between pigments m and n is contained in J(ω), that of coupling fluctuation between m and n and the site energy fluctuation of k in J(ω), and the correlation between fluctuations of excitonic couplings m ↔ n and k ↔ l in J(ω).

Computational Details

The sets of atomic partial charges used in the CDC and TrEsp methods were taken from earlier work[73] and are based on QC calculations on geometry-optimized BChla in vacuo. For the computation of the ESP of the charge density of ground and first excited states as well as of the transition density between gound and first excited state, (TD)DFT was employed with the B3LYP exchange correlation functional. We note that, in the CDC and TrEsp methods, no transition energies from QC but only the charge sets are used. In this context, a critical evaluation of the performance of various QC methods is only possible by comparison of simulated with measured optical spectra. In the case of the FMO protein, it was found that the charge sets inferred from B3LYP allow for a quantitative description of optical spectra.[9,10] Therefore, these charge sets are used here as well.

The atomic partial charges representing the ground state charge density of the protein were taken from the CHARMM force field.[77,78] All calculations are based on the 1.3 Å resolution crystal structure of the FMO protein of Prosthecochloris aestuarii by Tronrud et al.[6] For simplicity, we restricted the analysis to one monomeric subunit of this complex (Figure 1) and did not include the eighth BChla (apo-form). The use of the apo-form is justified, as the eighth BChla is likely to be absent in most samples used for spectroscopy.[10] Hydrogen atoms were added by using CHARMM.[77,78] A three-step geometry optimization of the whole complex was performed with steepest descent, conjugate gradient and the Newton–Raphson methods using CHARMM. The force field parameters of BChla were adopted from Marchi and co-workers,[79] who based the parametrization on DFT computations.[80] The published force field of BChla[79] was intended for use with AMBER94[81] and was adapted for the current application to the format and functional forms of CHARMM.[78] Finally, the NMA was performed with CHARMM, revealing 19 392 normal-mode frequencies and corresponding eigenvectors. The frequencies of the first six normal modes, which correspond to translation and rotation of the whole FMO protein, were zero. All other eigenvalues of the NMA were positive, indicating a stable minimum of the optimized structure.

Figure 1

Monomeric subunit of the FMO protein of Prosthecochloris aestuarii(6) containing seven BChla pigments (apo-form).

Theory of Optical Spectra and Exciton Dynamics

The Hamiltonian of the PPC in eq 1 is divided into three partsthe exciton Hamiltonian Hex, the vibrational Hamiltonian Hvib, and the Hamiltonian Hex–vib containing the exciton–vibrational coupling. The exciton Hamiltonian readsThe exciton matrix H(0) contains in the diagonal the site energies and in the off-diagonal the excitonic couplings, obtained from the CDC and Poisson-TrEsp methods, respectively. By diagonalizing this matrix, delocalized exciton statesare defined, where the coefficients c( are obtained from the eigenvectors, and the eigenvalues ℏω are the excitation energies of the delocalized states.

The exciton–vibrational Hamiltonian reads, using eq 5,where the exciton–vibrational coupling constants gξ(M,N) of the delocalized (exciton) states were introduced aswith the coefficients c( of exciton states and the local coupling constants gξ(m,n) of fluctuations of site energies (m = n) and excitonic couplings (m ≠ n) obtained here from a NMA (eq 9).

Finally, the vibrational Hamiltonian Hvib is that of uncoupled harmonic oscillatorswith the normal-mode frequencies ωξ. Hence, all parts of the Hamiltonian are fully parametrized by structure-based simulations.

Time-local partial ordering prescription (POP) non-Markovian density matrix theory[11] is used to describe linear absorbance and exciton dynamics following an ultrashort pulse excitation. This theory contains an exact treatment of the diagonal parts of the exciton–vibrational coupling in the basis of delocalized states, characterized by the coupling constants gξ(M,M) (eq 15), and treats the off-diagonal parts gξ(M,N) with M ≠ N by perturbation theory. A microscopic justification for this approximation based on the present NMA will be provided below.

Linear Absorbance

The absorbance spectrum α(ω) is obtained from the homogeneous absorbance spectrum α(hom)(ω) viawhere ⟨...⟩dis denotes an average over static disorder in site energies. An independent variation of the site energies and Gaussian distribution functions, centered around the static-structure site energies, is assumed, and the disorder average is performed numerically by a Monte Carlo method, as usual.

The homogeneous absorption spectrum readswhere the transition dipole moment μ⃗ of the optical transition between the ground state and the Mth exciton state is given asHere, μ⃗ is the transition dipole moment of the local optical transition at site m.

The line shape function D(ω) is obtained from the non-Markovian time-local POP density matrix theory within secular approximation and Markov approximation for the off-diagonal elements of the exciton–vibrational coupling[11]The time-dependent function G(t) describes the vibrational sideband of the exciton transition and is related to the spectral density J(ω), defined in eq 10, via (M = N)where n(ω) is the Bose–Einstein distribution function of vibrational quantaThe dephasing time constant τ in eq 20 contains the lifetime broadening due to exciton relaxation between exciton state |M⟩ and the other exciton states |N⟩with the Redfield rate constant

The frequency ω̃ in eq 20 is the transition frequency between the ground state |0⟩ and the exciton state |M⟩ that contains a renormalization due to the diagonal and off-diagonal parts of the exciton–vibrational coupling[11,30]where Eλ( is the reorganization energy of the Mth exciton state defined asand denotes the principal part of the integral. Note that J(ω) = 0 for ω < 0.

Excitation by an Ultrashort Optical Pulse

We consider excitation by a delta-shaped pulse E(t) = Ae⃗δ(t) with amplitude A and polarization vector e⃗. Using second-order perturbation theory for the coupling of the pulse with the system, and neglecting any nuclear dynamics during the action of the pulse, it is seen that the δ-pulse prepares the system in a statewhere we have included an average over the orientation of the complex with respect to the external field, and α is the angle between transition dipole moments μ⃗ and μ⃗.

Within the secular approximation, the equations of motion for the off-diagonal elements of the density matrix are obtained as[11]with the time-dependent function F(t)The correlation function C(t) is related to the spectral density J(ω) in eq 10 byA Markov approximation is applied to the off-diagonal parts (L ≠ M and L ≠ N) in eq 29 by setting the integration limit t = ∞ for those terms. Hence, F(t) becomeswhere denotes the real part, and the half-sided Fourier transform C̃(ω) of C(t) was introduced asThe real part C̃(Re)(ω), up to a factor of 1/2, equals the Redfield rate constant k for exciton relaxation introduced in eq 24and the imaginary part C̃(Im)(ω) is obtained from the real part by a principal part integral[30]By introducing the above equations into eq 29, the resulting expression into eq 28, and performing the time integrations, we obtainwithwhere the function G(t) was introduced in eq 21, and ω̃ = ω̃ – ω̃, with the latter being defined in eq 25.

In secular approximation, the dynamics of diagonal elements of the density matrix is decoupled from that of the off-diagonal elements reading simplywith the rate constants defined in eq 24.

Results

Spectral Density

In Figure 2, the diagonal elements of the spectral density, i.e., the J(ω) describing the fluctuation of site energies of the seven BChla (m = 1, ..., 7) in the monomeric subunit of the FMO protein, as obtained from the NMA are compared to the spectral density J(ω) as extracted earlier[11] from FLN spectra of the B777 complex. The latter J(ω) is similar in shape to the vibrational sideband measured in FLN spectroscopy on the FMO protein[12] (see below and the discussion in ref (8)) and has been successfully applied to a large number of PPCs.[3] Its amplitude is determined from a fit of the temperature dependence of linear absorbance spectra (see below), giving a Huang–Rhys factor S = ∫0∞ dωJ(ω) = 0.42. Similar Huang–Rhys factorsare obtained for the seven pigments from the NMA with variations between S5 = 0.19 and S2 = 0.54 and an average value of 0.39. The J(ω) are similar in shape for the different pigments, but there are systematic deviations from the experimental J(ω). At small frequencies, the NMA spectral densities are larger, whereas at larger frequencies they are somewhat below the experimental values. In Figure 3, the average NMA spectral density of site energy fluctuationsis compared with the experimental J(ω) of the B777 complex[11] and that of the FMO protein, estimated from the vibrational sideband measured in FLN.[12] For better comparison, the experimental data are rescaled to the same Huang–Rhys factor (area) as the J̅(ω) obtained from the NMA. As noted above, the NMA overestimates the low-frequency contributions to the spectral density. On the basis of the comparison in Figure 3, we will define further below a factor f(ω) that shall correct for the missing anharmonic terms of the force field in the NMA in an effective way. Later, results will be compared obtained with the original and the corrected spectral densities.

Figure 2

Spectral densities J(ω) describing the fluctuation of site energies of the pigments m = 1···7. The J(ω) are represented as histograms, where the ω-axis has been discretized in steps of 5 cm–1. For comparison, we show J(ω), the shape of which was obtained from a fit of fluorescence line narrowing spectra of a one-pigment system, the B777 complex (solid black line),[11] and the integral coupling strength, i.e., the Huang–Rhys factor S = 0.42, was obtained from a fit of the temperature dependence of the linear absorbance spectrum of the FMO protein (Figure 6). The individual Huang–Rhys factors S (eq 38) of the pigments, obtained from the NMA, are shown as well.

Figure 3

Comparison of the average spectral density J̅diag(ω) (eq 39) of site energy fluctuations of the seven pigments, shown separately in Figure 2, with the experimental spectral density of the B777 complex[11] and the FMO protein.[12] The latter two have been rescaled to the Huang–Rhys factor S̅ = 0.39 of J̅diag(ω) (eq 39), resulting from the NMA, for better comparison.

Figure 6

Linear absorbance spectrum at three different temperatures. The experimental data[12] (left) are compared to simulations obtained with the directly calculated spectral density J(ω) (middle) and with the corrected spectral density Jc(ω) (right).

The correlation in site energy fluctuations between BChla at sites m and n is described by the J(ω) shown in Figure 4 for those pigment pairs with the largest correlations. To quantify correlations, we introduce a generalized Huang–Rhys factorHere, we use the absolute magnitude of the coupling constants to avoid a compensation of positive and negative contributions resulting from correlated and anticorrelated fluctuations, respectively. Note that, for m = n = k = l, the generalized Huang–Rhys factor reduces to the ordinary Huang–Rhys factor S = ∑ξ (gξ(m, m))2. Interestingly, the S for the correlation in site energy fluctuations (Figure 4) are in the same order of magnitude as the S of the site energy fluctuations (Figure 2), but the amplitudes are somewhat weaker at larger frequencies.

Figure 4

Spectral densities J(ω) describing the correlations in site energy fluctuations of pigments m and n. The pigment pairs with the largest correlation strengths, characterized by the generalized Huang–Rhys factors S (eq 40), are shown.

The fluctuations of excitonic couplings are characterized by the J(ω) shown in Figure 5 for those pigment pairs with the largest exciton–vibrational coupling strength. The latter has been evaluated according to the generalized Huang–Rhys factorThe largest Huang–Rhys factor S34 = 0.031 is still about 1 order of magnitude smaller than those of the site energy fluctuations in Figure 2 and their correlations in Figure 4. The spectral densities of the correlations that involve fluctuations of excitonic couplings (site-coupling J (k ≠ l) and coupling-coupling J (m ≠ n, k ≠ l)) are in the same order of magnitude as the spectral density of the coupling fluctuations J (m ≠ n), as shown in the Supporting Information. The J (k ≠ l) is largest if m equals k or l (Supporting Information, Figure S1), and J (m ≠ n, k ≠ l)) is largest for those combinations of pigment pairs, which share one pigment (Supporting Information, Figure S2). From the J (m ≠ n, m ≠ l), one can get a rough picture of the characteristics of the low-frequency normal modes. For example, the fluctuations of excitonic couplings H12(0) and H16(0) are anticorrelated throughout the spectrum of normal modes, indicating an overall oscillating motion of pigment 1 between pigments 2 and 6.

Figure 5

Spectral densities J(ω) characterizing the fluctuations of excitonic couplings between pigments m and n are shown for those pigment pairs with the largest fluctuations, characterized by the Huang–Rhys factors S (eq 41).

The spectral density J(ω) determined above from a NMA is used to calculate linear absorbance spectra of the FMO protein as a function of temperature. The site energies H(0) obtained recently by using the CDC method and a refinement fit (ref (10), Table 1, column 2 (apo)) are used as well as the excitonic couplings (Supporting Information of ref (10), Table 5, please note the misprint concerning the coupling H12(0) = −94 cm–1, where the minus sign was missing in the Supporting Information) obtained from a point-dipole approximation with an effective dipole strength of 29.8 D2 for BChla. The point-dipole approximation has been verified by comparison with the TrEsp method, and the effective dipole strength was determined by comparison with Poisson-TrEsp values.[8,56] The direction of the transition dipole of BChla is taken along the NB–ND axis (e.g., ref (73)). We use the inhomogeneous broadenings of the BChla pigments from our earlier publications.[9,10] The width of the inhomogeneous distribution function is 60 cm–1 for pigments 1, 3, and 4, 100 cm–1 for pigment 2, and 120 cm–1 for pigments 5–7. This choice has been motivated by the presence of mobile water molecules in the neighborhood of the pigments with a larger inhomogeneous broadening.

The nonbonded interactions of molecular mechanics force fields describe the soft degrees of freedom of the protein, which are responsible for the conformational flexibility of the macromolecule. Therefore, it can be expected that the harmonic approximation of a NMA is particularly critical for the latter degrees of freedom. We will assume in the following that anharmonic effects can be taken into account by a frequency-dependent factor f(ω) resulting in the corrected spectral density Jc(ω) that is obtained from the NMA spectral density J(ω) asThe correction factor f(ω), which in a simple picture takes into account a change in the (effective) density of vibrational states, is obtained by comparing the average diagonal part of the NMA spectral density J̅diag(ω) defined in eq 39 with the experimental spectral density Jexp(ω) (Figure 3) asJexp(ω) is obtained, up to a constant factor, from the vibrational sideband of the FLN spectra of the FMO protein (the red solid line in Figure 3). The constant factor, that is the experimental Huang–Rhys factor, is obtained from an analysis of the temperature dependence of linear absorbance (Figure 6). The rationale behind this definition of f(ω) is that the lowest exciton state of the FMO protein has a dominant contribution from one pigment, and therefore, the vibrational sideband observed in fluorescence line narrowing spectra at low T is mostly influenced by the diagonal part of the spectral density. The average over all pigments in J̅diag(ω) (eq 39) is taken to make the corrections factor less dependent on site-specific details that might be not accurate enough.

With both spectral densities, original and corrected, the change of the absorbance spectrum with temperature between 4 and 300 K can be qualitatively understood (Figure 6). At higher temperatures, the increased dephasing introduced by exciton relaxation and vibrational excitations leads to a broadening of the absorbance peaks. Overall, the corrected spectral density provides a slightly better description of the temperature dependence. From this temperature dependence, a Huang–Rhys factor of S = 0.42 was determined for the average diagonal part of Jc, that is (1/7)∫0∞ dω∑ Jc = 0.42, in good agreement with earlier studies, which gave S = 0.45[12] and S = 0.5.[8] We note that the NMA includes only the intermolecular part of the spectral density, the Huang–Rhys factor of which has been estimated to S = 0.3 from the FLN spectra on the FMO protein of P. aestuarii.[12] The same intermolecular Huang–Rhys factor was measured for the FMO protein of C. tepidum.[82] The average Huang–Rhys factor S̅ = ∑17S/7 = 0.39 of the pigments obtained without any correction directly from the NMA is in good agreement with this value.

Population Transfer and Decay of Coherences

In the next step, the two spectral densities are used to simulate exciton relaxation at 77 K, where the initial population of exciton states is assumed to be prepared by an ultrashort (δ-shaped) optical pulse acting at t = 0. In this calculation and also in the calculation of dephasing below, the homogeneous behavior is studied; i.e., no inhomogeneous broadening is included. As seen in Figure 7, the excitons equilibrate faster when the corrected spectral density Jc(ω) is used. The time scale obtained for the latter is in agreement with experimental data,[83] whereas the exciton dynamics obtained for the original spectral density J(ω) contains a slow component that is not observed experimentally.

Figure 7

Population of exciton states after δ-pulse excitation at t = 0, calculated with the original spectral density J(ω) (left part) and the corrected one Jc(ω) (right part) at T = 77 K. The solid curves show results obtained by taking into account the full spectral density J(c)(ω). The dashed and dotted curves show simulations, where only uncorrelated site energy fluctuations (J(c)(ω)) and uncorrelated site energy and coupling fluctuations (J(c)(ω)), respectively, were included. Please note that the sum probability of excited state populations is a constant that depends on the amplitude of the external field, which was chosen small enough to justify the second-order perturbation theory used for the initial populations (eq 27).

In addition, it was investigated which parts of the spectral density J(c)(ω) are important for exciton relaxation by setting certain other parts to zero. Whereas the solid lines in Figure 7 were obtained for the full spectral densities, the dashed lines represent the case for which only the diagonal parts J(c)(ω) were included. The dotted lines were obtained by including only spectral densities J(c)(ω), which for m = n describe the fluctuation of site energies and for m ≠ n that of excitonic couplings. Note that in the latter two cases all correlations are neglected. It is clearly seen that exciton relaxation is determined by the diagonal part J(ω) of the spectral density and that neither the fluctuations in excitonic couplings nor the correlations between the different types of the spectral density play any significant role for the energy transfer.

Next, we study the decay of coherences created by the ultrashort pulse. As before, we consider the original and the corrected spectral densities. In Figure 8, the decay of the density matrix element ρ12 describing the coherence between the two lowest exciton states in the FMO protein is investigated for four different temperatures between 77 and 277 K as in the experimental study.[68] With increasing temperature, the coherence decays faster, as expected. The decay, in particular at low T, is faster for the original spectral density than for the corrected one. The latter, however, describes the experimental data better (see Figure 3 in ref (68)). There is a second set of red dotted lines present in Figure 8 that were calculated without taking into account correlations. However, the agreement with the full calculations is so complete that the two sets of curves lie practically on top of each other. We find that for both spectral densities the influence of correlations is negligible.

Figure 8

Dephasing of coherences ρ12, created by a δ-pulse acting at t = 0, between the two lowest exciton states in dependence on temperature, obtained for the original (left part) and the corrected (right part) spectral density, J(ω) and Jc(ω), respectively. The real part of ρ12 is shown. The solid black lines show calculations obtained using all parts of the spectral densities, and for the red-dotted lines, all correlations were neglected by including only J(c)(ω) and setting all other elements to zero. The solid and dotted lines lie practically on top of each other.

In an attempt to find out which properties the spectral density of correlation in site energy fluctuations Jc(ω) (m ≠ n) should have to allow the protein to protect coherences between different exciton states, we replaced Jc(ω) by the average of the diagonal parts of the spectral density of the two pigments, scaled by a constant αwhere the spectral densities on the r.h.s. were taken from the NMA (corrected as described above). As shown in the upper part of Figure 9, for negative α values, the coherences are even further suppressed, whereas positive values of α prolong their lifetimes. For α = 1, the coherences live longer than 10 ps, even at 277 K, whereas, for absolute magnitudes of α smaller than 0.3, there is almost no influence of the correlations in site energy fluctuations. In the lower part of Figure 9, the corresponding population ρ11 of the lowest exciton state is shown for the different α-values. As seen there, large positive α-values, which lead to long coherence times, strongly suppress exciton relaxation.

Figure 9

Dephasing of coherences ρ12 (upper part) and population of the lowest exciton state ρ11 (lower part) at T = 277 K, following excitation by a δ-pulse acting at t = 0, for different forms of the site energy correlation part of the spectral density Jprotect(ω), created artificially from the site energy fluctuation parts Jc(ω) and Jc(ω) of the NMA, as described in detail in the text (eq 44).

Finally, we investigate how far away the actual spectral density Jc(ω) obtained from the NMA is from the case of ideal protection of coherences as defined in eq 44 for α = 1. For this purpose, we introduce the following dephasing coefficient:which is essentially 1 – α and does not depend on the correction from J(ω) to Jc(ω) (i.e., the factor f(ω) in eq 42 cancels out). If there is no correlation in site energy fluctuations, we have Jc(ω) = 0 for m ≠ n and the resulting d(ω) = 1 indicates a strong dephasing. For anticorrelated fluctuation in site energies, it holds that J(ω) < 0, and the resulting d(ω) > 1 shows even stronger dephasing (corresponding to negative α values). If, on the other hand, it holds that 2Jc(ω) = Jc(ω) + Jc(ω), we have d(ω) = 0 corresponding to α = 1, that is, weak dephasing.

In Figure 10, the dephasing coefficients obtained from the normal mode spectral density are shown. The coefficients fluctuate around 1 at all frequencies, with α = 1 – d values much smaller than 0.3 for almost all normal modes, indicating a strong dephasing. In many cases, the correlations at some frequency seem to be compensated by anticorrelations at another frequency.

Figure 10

Dephasing coefficients d(ω), defined in the text (eq 45), for those pigment pairs m and n with the largest correlations in site energy fluctuations.

Discussion

Structure-Based Modeling of Pigment–Protein Complexes

The present study represents one further step toward an entirely structure-based modeling of excitation energy transfer and optical sepctra of PPCs. By combining the previously developed CDC and TrEsp methods for the computation of site energies and excitonic couplings, respectively, with NMA, information about the spectral density can be obtained. It is important to note that, by construction, this approach allows only to calculate intermolecular exciton–vibrational couplings, i.e., couplings of a transition on pigment m with the vibrations of the protein and other pigments n ≠ m. The coupling between a transition on pigment m with its own vibrations cannot be obtained. While this appears to be a drawback at first glance, it actually turned out to be an advantage. This seemingly counterintuitive result has its roots in one major problem of MD-based computations of the spectral density. The latter approach involves the computation of time correlation functions that in turn are based on transition energies computed directly with QC on snapshots of the MD trajectory. As a consequence, the correlation functions suffer from systematic errors due to the inherent limitations of QC calculations of excited states (depending on the electronic structure methods employed) and suboptimal pigment geometries, as mentioned in the Introduction. Apparently, these errors result in a severe overestimation of the contributions from intrapigment exciton–vibrational couplings to the spectral density predominantly (but not exclusively) in the high-frequency range, masking the actually more important contributions from intermolecular couplings being predominant at lower frequencies. By inevitably omitting the intrapigment part, the CDC/TrEsp/NMA approach opens up the view on the intermolecular couplings. Nonetheless, the omission of intramolecular vibrational excitations remains a problem that we discuss below. Another problem of the MD-based methods is also related to the computation of time correlation functions. In order to resolve the low-frequency part of the spectral density along with eventual correlations of site energy fluctuations, long simulation times are required. This problem does not show up in the NMA-based approach, as the spectral density is obtained directly. For the same reason, artifacts due to the classical description of high-frequency vibrational motions occurring in MD simulations are not an issue in the NMA method.

The present analysis provides a qualitatively correct description of the shape and the amplitude of the spectral density, as the comparison of the diagonal parts with spectral densities extracted from experimental data suggests (Figures 2 and 3). However, at very low frequencies, the normal mode spectral density is systematically too large, and at higher frequencies, it is lower than the experimental data. This deviation can have several reasons. First of all, many vibrational degrees of freedom of the PPC contributing to the low-frequency region are soft and anharmonic, so that a NMA might be an oversimplification for some of them. Further, the force field parameters are not necessarily optimal for the analysis of vibrations and so may contain uncertainties that lead to a systematic deviation of the calculated spectral density from the experimental data.

In addition, the neglect of intramolecular vibrational excitations certainly contributes to the underestimation of the high-frequency part of the spectral density. However, most of the pigment modes have very small Huang–Rhys factors[51,84] and the high-frequency range is less important for excitation energy transfer. Therefore, we do not consider this approximation to be critical. We note that, in principle, these contributions could be included by performing QC-based NMA of the excited and ground states of BChla in vacuo, as shown by Jing et al.[51] and for a different system (a pheophorbid a molecule dissolved in ethanol) by Megow et al.[85] Jing et al.[51] performed QC calculations on geometry-optimized BChla in vacuo, and therefore, the Huang–Rhys factors obtained do not suffer from the geometry-distortion problem. Megow et al.[85] performed QC calculations on the geometry-optimized chromophore in a homogeneous dielectric, representing the solvent, revealing the Hessians of the ground and excited states. These calculations were combined with MD simulations of the chromophore in an explicitly described solvent. Using a suitable mapping procedure, the deviations of the chromophore geometry from the equilibrium structure were determined along the MD trajectory, and, using the Hessians, the fluctuation of the chromophore’s excitation energy calculated. The good agreement of the resulting absorbance spectrum with experimental data suggests that the geometry–distortion problem cannot be severe. Therefore, the latter approach might be a good compromise between a simple QC-based NMA in vacuo and a very involved QM/MM approach of the PPC, where the force constants are obtained directly from ab initio quantum chemical calculations, to avoid the geometry-mismatch problem.

In the present work, we decided to correct for the missing anharmonicities, the uncertainties of force field parameters, and the missing contributions from intramolecular modes by introducing a scaling function f(ω) that is determined on the basis of a comparison of the average diagonal part of the spectral density with experimental data. It was demonstrated before that, in the spirit of the central limit theorem and a second-order cumulant expansion, even a strongly anharmonic system might still be described by an effective spectral density of harmonic oscillators.[86,87] Our scaling function f(ω) may be seen as one way to create such an effective spectral density. As the calculations of linear absorbance (Figure 6), exciton relaxation (Figure 7), and dephasing of coherences (Figure 8) show, in all three cases, the corrected spectral density provided an improved description of the experiment. Note that the correction of the spectral density according to eq 42 does not affect the relative amplitudes of the different parts of the spectral density.

At present, we believe that the most serious approximation of our NMA is the harmonic treatment of the nonbonded interactions in the molecular mechanics force field. We are currently performing CDC/TrEsp/MD calculations to investigate this point in detail.

Concerning the decay of quantum beats measured in 2D spectroscopy[68] and our simple calculation of the decay of the exciton coherence ρ12 between the two lowest exciton states, which reproduces the observed temperature dependence of the experimental decay, we have to add a word of caution. Recent results obtained by Christensson et al.[88] and Caycedo-Soler et al.[89] suggest that the static disorder present in this system, if it is uncorrelated, can lead to significant additional decay of the quantum beating observed in the 2D spectra. This result implies that the evolution of off-diagonal density matrix elements cannot be directly compared with the oscillations in the 2D spectra. On the other hand, Kreisbeck and Kramer[21] reported a large similarity between these two quantities (as seen from a comparison of Figures 3 and 4 of their paper).

It will be interesting to use the present model for a direct calculation of the 2D signal, including also inhomogeneous broadening.[21,89,88] An open question in that respect is whether there are correlations in the static disorder of the site energies of different pigments. The correlation functions in Figure 4 have a large amplitude at small frequencies. The lowest frequency (2.7 cm–1) corresponds to an oscillation period of 2 ps, which is already in the range of the time scale of exciton equilibration. A NMA on the whole FMO trimer would allow us to reach even lower frequencies (longer time scales). In this way, it will be possible to check also the assumption of uncorrelated static disorder in site energies, commonly used in the theory of optical spectra and excitation energy transfer, with a NMA. Considering the fact that the present NMA of the FMO monomer reveals the largest correlations in site energy fluctuations for low frequencies, it seems likely that the correlations in static disorder are significant. If so, these correlations have to be included in the calculation of the 2D signal and can be expected to lead to long dephasing times.[88,89] We note that Christensson et al.[88] found that, in the absence of static site energy correlations, intramolecular vibronic transitions of the pigments are able to provide the means for correlated energy level fluctuations and thereby long-lived coherences.

To check this idea, Engel and co-workers[90] recently performed 2D experiments on BChla in two different solvents but did not find the characteristic beating pattern that they observed before in the FMO protein. This result seems to suggest that the beating indeed involves coherent motion of excitonic wavepackets, rather than vibrational wavepackets. On the other hand, the coupling between intra- and intermolecular vibrational modes could be much stronger in the solvent than in the FMO protein. Butkus et al.[91] recently reported model calculations suggesting that vibrational quantum beats in 2D spectra can be suppressed by a strong coupling to environmental modes.

Finally, we note that Kreisbeck and Kramer[21] found that a small slope of the spectral density at small frequencies is required for long-lived coherences to occur. The faster dephasing of coherences in the left part of Figure 8 obtained for the original spectral density, that has a larger amplitude in the low-frequency region than the corrected one, is in agreement with this result. Exciton relaxation is faster for the corrected spectral density than for the original (Figure 7). This result demonstrates that it is not the exciton-relaxation induced dephasing, described by the factor e– in eq 35, but rather the pure dephasing term e– in this equation that dominates the overall dephasing of coherences. The latter corresponds to vibrational excitations in excitonic potential energy surfaces, which is stronger if the spectral density has large low-frequency components, because the latter can be most easily excited thermally.

Role of Correlations

A main goal of the present work was to quantify the correlations among the different contributions to the spectral density and to investigate their role in energy transfer. Despite the rather strong amplitude of the correlation in site energy fluctuations present in the spectral density (Figure 4), their function in energy transfer and dephasing of coherences is practically zero (Figures 7 and 8, respectively). There exists an optimal form of the correlation part of the site energy fluctuations that would allow the protein to protect excitonic coherences even at 277 K (upper part of Figure 9). However, the analysis of dephasing coefficients, introduced in eq 45, of the coherences in the FMO protein clearly shows that the protein has realized a correlation regime that is far from such a protection scenario. Two questions arise: (i) Is there an advantage for the light-harvesting function of the protein, if coherences between exciton states dephase fast? (ii) What is the microscopic origin of the fast dephasing? The answer to the first question is yes. A protection state of coherences would considerably slow down exciton relaxation and thereby energy transfer (lower part of Figure 9), as found earlier also by using phenomenological models for the spectral density.[35,36] On the other hand, we know that exciton delocalization may persist over long times and is used, e.g., in the photosystem of purple bacteria to create an excitation energy funnel. Therefore, we have to conclude that even an uncorrelated fluctuation of site energies appears to be too weak to destroy a delocalization of excitons over a certain number of pigments. Upon optical excitation, these delocalized states are created and excitation energy relaxes within the domains of delocalized states and is transferred also between different such domains. A fast relaxation in the domains helps to direct the excitation energy to the domain with the lowest energy. The same mechanism that would protect exciton coherences created by an ultrashort pulse would inhibit exciton relaxation in the domains of strongly coupled pigments. Therefore, it is not surprising that the protein does not protect these coherences. It rather protects interpigment coherences (i.e., between localized excited states), as will be discussed further below.

Concerning the second question about the microscopic origin of the fast dephasing, we believe that the answer is connected to the inhomogeneous charge distribution of the protein and the long-range Coulomb interaction. To make this point clear, let us consider a single normal mode ξ. Our ideal protection scenario would require gξ(m,m)gξ(n,n) = (gξ(m,m)gξ(m,m)+gξ(n,n)gξ(n,n))/2 or (gξ(m,m) – gξ(n,n))2 = 0, which is fulfilled, if gξ(m,m) = gξ(n,n). Hence, in order to maximally protect coherences, the protein has to create a situation, in which the pigments m and n are coupled to every normal mode in exactly the same way. A strong dephasing, therefore, can be accomplished in two ways: (i) by localizing the normal modes such that every pigment feels its own independent heat bath and (ii) for delocalized normal modes ξ by strongly varying the local coupling constants from site to site and excitonic coupling to excitonic coupling. Since low-frequency modes in general involve the motion of larger parts of the protein, one might expect to find smaller dephasing coefficients for small frequencies ω. If at all, this expectation is only fulfilled for d12 and d34 (Figure 10) and only for very low frequencies. The strikingly small variation of dephasing coefficients with frequency suggests that (ii) is the dominating mechanism, responsible for fast dephasing of coherences and fast energy transfer. As a consequence of the inhomogeneous distribution of charges and vibrating atoms, every pigment feels a different Coulomb field and hence the local coupling constants will not be identical. It is interesting that the inhomogeneous charge distribution is also responsible for the creation of an excitation energy funnel in the FMO protein.[9,10] Thus, the same physical mechanism that directs excitation energy flow also causes fast dephasing of coherences and thereby exciton relaxation.

As pointed out by Huo and Coker,[44] a strong correlation in excitonic coupling fluctuations could change the simple picture of long-lasting coherences and fast energy transfer excluding each other. We have, therefore, estimated the maximum contribution of the part of the spectral density containing the coupling fluctuations in the following way. Assuming that the NMA cannot describe the polarizability of the protein, we would have to use εeff = 3, determined from static-structure site energy calculations,[10] instead of εeff = 1.25, used so far in the CDC part of the NMA. This change would increase the relative strength of the coupling fluctuations by a factor of (3/1.25)2 = 5.8, bringing it closer to that of the site energy fluctuations but leaving it still smaller by an average factor of about 5. In this case, the exciton relaxation is influenced somewhat stronger by the coupling fluctuations than in Figure 7, but the influence of correlations is still negligible. This result shows that the present system is far away from the scenario that Huo and Coker identified for a more complicated relation between energy transfer and dephasing. In other words, the dominating influence of site energy fluctuations, found here, implies that a fast decay of exciton coherences is needed for a fast energy transfer.

Implications for the Theory of Optical Spectra and Excitation Energy Transfer

The Spectral Density

In the absence of any microscopic analysis of the spectral density, the most reasonable and simple approximation, used in many previous theories, was to neglect all parts of the spectral density J(ω) except for the diagonal part J(ω) describing the uncorrelated fluctuation of site energies. The present work shows that the fluctuations of excitonic couplings and the related correlations are small compared to the site energy fluctuations and, therefore, can indeed be neglected. The correlations in site energy fluctuations are of the same order of magnitude as the site energy fluctuations. However, since their influence on the exciton dynamics is negligible due to the inhomogeneous charge distribution of the protein, as explained in detail above, this part of the spectral density may also safely be neglected. Although only shown for the FMO protein here, the inhomogeneous charge distribution of other PPCs suggests that they should not behave in a qualitatively different way. Therefore, we conclude that the present work provides a quite general microscopic justification for the above standard assumptions.

Another assumption of theories was to assume the same coupling strength for all pigments, i.e., to neglect the dependence of J(ω) on m. The present approach, in contrast, shows considerable variations of the related Huang–Rhys factors varying by up to 50% (Figure 2). Since in line-narrowing spectroscopy only the low-energy exciton states can be investigated, it is an important contribution of microscopic theory to provide information about the exciton–vibrational coupling of those pigments contributing to high-energy exciton states. At the present level of agreement between directly calculated spectral density and experimental data, it is too early to draw definite conclusions about quantitative differences of local exciton–vibrational coupling strengths of particular pigments. We leave this question for future studies that include anharmonicities of the force field like the CDC/TrEsp/MD method.

Is There a Small Parameter for Perturbation Theory?

As noted before, when modeling optical spectra and excitation energy transfer in PPCs, one is confronted with the problem that the exciton–vibrational (pigment–protein) as well as the excitonic (pigment–pigment) coupling are in the same order of magnitude. As an example, for the present system, the average reorganization energy of local fluctuations in site energies E̅λ = ∫0∞ dωJ̅(ω) amounts to 46 cm–1 and the absolute magnitudes of nearest neighbor excitonic couplings vary between 30 cm–1 (pigments 2 and 3) and 95 cm–1 (pigments 1 and 2). Such scenarios can be treated by non-perturbative methods (see the Introduction); however, these methods are rather involved. In other, numerically less costly approaches, like polaron and non-Markovian density matrix theory, one tries to find a representation, in which a small parameter appears that can be used for a perturbative treatment. In the time-local (POP) non-Markovian density matrix theory, used in the present work, this small parameter is the off-diagonal part of the exciton–vibrational coupling in the representation of delocalized exciton states, characterized by the coupling constants gξ(M,N) (M ≠ N) in eq 15.

Thus, by using a representation of delocalized excited states, we transform a problem, where we had competing excitonic couplings V = H(0) and reorganization energies Eλ = ∑ξ  ℏωξgξ2(m,m) of the local exciton–vibrational coupling, into a problem with only exciton–vibrational coupling constants gξ(M,N), which are still difficult to treat non-perturbatively. However, as Figure 11 shows, the off-diagonal elements gξ(M,N) (M ≠ N) are found to be smaller than the diagonal elements gξ(M,M). We take this difference as a justification of using a perturbation theory for the off-diagonal parts of the exciton–vibrational coupling (including a Markov approximation), while using an exact summation of the diagonal parts as provided by the time-local non-Markovian (POP) density matrix theory. It now becomes clear why the alternative COP theory (see the Introduction) gave less accurate results than POP,[11,33] since COP contains a partial summation of both the diagonal as well as the off-diagonal parts of the exciton–vibrational coupling. A comparison of the more accurate POP with a numerically exact theory would be helpful to judge about the conversion of the present perturbation theory of the off-diagonal exciton–vibrational coupling. A close inspection of the red wing of the absorbance spectrum calculated with the corrected spectral density in the right part of Figure 6 reveals a somewhat too strong broadening with increasing temperature when compared with the experimental data shown in the left part of this figure. This effect was first noted by Novoderezhkin et al. and was termed Redfield artifact.[92] The origin could lie in the perturbation theory used for the off-diagonal parts of the exciton–vibrational coupling.

Figure 11

Coupling constants gξ(M,N) of exciton–vibrational coupling of delocalized exciton states, obtained from eq 15, using the microscopic coupling constants gξ(m,n) from the NMA as a function of normal mode index ξ for the first 4000 normal modes. The black solid line shows the corresponding normal-mode frequencies ω = ωξ.

How Can the Protein Protect Electronic Coherences?

The dominance of the diagonal parts of the exciton–vibrational coupling in the basis of delocalized states found above (Figure 11) shows that exciton states stay delocalized over a certain number of pigments, as determined by their differences in site energies and static disorder. A dynamic localization of excitons would be caused by the off-diagonal parts which are, however, small. By keeping the pigments at the right distance, the protein is able to allow for a certain delocalization of excited states and thereby for prevailing interpigment coherences, that is, density matrix elements ρ that are off-diagonal in the basis of localized excited states. Even when excitons are relaxed, where we have approximately ρ ∝ δ exp{−ℏω/kBT}, interpigment coherences are present, since ρ = ∑  c(c(ρ. Note that the dominance of the diagonal parts gξ(M, M) of the exciton–vibrational coupling in the delocalized basis justifies the assumption of a Boltzmann equilibrium in that basis. As we have seen in the present calculations, the correlations in site energy fluctuations have no influence on the relaxation dynamics of excitons (described by ρ) and the decay of coherences (ρ) between different delocalized states and, therefore, also not on the interpigment coherences ρ. An interesting question is: How could the protein change the delocalization of excited states? A dynamic localization is observed, if two pigments come so close that they can exchange, in addition to excitations, also electrons.[93] In this case, the mixing of exciton states with charge transfer states[94,95] leads to very strong exciton–vibrational coupling that can localize the excited states and also allows for nonradiative transitions to the ground state. Hence, we may conclude that the protein protects interpigment coherences by locating pigments at close enough distances such that the excitonic coupling is comparable with the difference in site energies and the amount of static disorder and, on the other hand, at large enough distances to prevent dynamic localization effects.

Summary and Conclusions

A method for the microscopic calculation of the spectral density of the pigment–protein coupling in light-harvesting systems was established that combines NMA with the CDC and TrEsp methods. The quality of the calculations was checked by comparison with spectral densities extracted from line narrowing spectroscopy, revealing good qualitative agreement concerning shape and amplitude. The correlations in site energy fluctuations obtained from the NMA are of similar magnitude as the site energy fluctuations, whereas the fluctuations of excitonic couplings and the related correlations are 1 order of magnitude smaller. Therefore, a microscopic justification was obtained for the standard assumption of neglecting this part of the spectral density. Concerning the correlations in site energy fluctuations, it was found that their influence on the relaxation dynamics of excitons and the dephasing of coherences between exciton states is practically zero. A detailed analysis of the dephasing process and of possibilities of its inhibition was performed. It reveals that the same mechanism that creates an excitation energy sink in this system, namely, the inhomogeneous charge distribution of the protein, is responsible for the fast dephasing of coherences and its insensitivity to correlations in site energy fluctuations. In this way, a directed energy transfer and a fast dissipation of the excitons’ excess energy become possible. The large amplitude correlations in site energy fluctuations found at low frequencies could be important for the interpretation of 2D spectra.