Coherence in carotenoid-to-chlorophyll energy transfer

The subtle details of the mechanism of energy flow from carotenoids to chlorophylls in biological light-harvesting complexes are still not fully understood, especially in the ultrafast regime. Here we focus on the antenna complex peridinin–chlorophyll a–protein (PCP), known for its remarkable efficiency of excitation energy transfer from carotenoids—peridinins—to chlorophylls. PCP solutions are studied by means of 2D electronic spectroscopy in different experimental conditions. Together with a global kinetic analysis and multiscale quantum chemical calculations, these data allow us to comprehensively address the contribution of the potential pathways of energy flow in PCP. These data support dominant energy transfer from peridinin S2 to chlorophyll Qy state via an ultrafast coherent mechanism. The coherent superposition of the two states is functional to drive population to the final acceptor state, adding an important piece of information in the quest for connections between coherent phenomena and biological functions.


Supplementary Methods 1: Quantum Chemical Methods
Exciton model. The excited states of the whole PCP system are described through an exciton model. The exciton HamiltonianĤ ex for a system of N pigments, each bearing n i excitations, reads: where E a i is the a-th excitation energy of pigment i (site energy), V ab ij is the electronic coupling between two transitions of different pigments. Diagonalization of the Hamiltonian matrix in Eq. S1 yields the energy levels of the exciton states E K and the exciton coefficients C K ia . The excitonic Hamiltonian of a PCP monomer includes the S 2 bright states of the eight Pers, and the Q and Soret (B) bands of the two chlorophylls (Chls), resulting in 16 total states.
Spectral densities of the exciton-phonon coupling. The spectral densities (SD) of the different pigments and transitions were modeled as a sum of one cubic continuous contribution and M discrete contributions: where i and a indicate the pigment and transition of interest, respectively, k runs on the vibrational modes of the pigment, λ c is the continuous contribution to the reorganization energy, and ω c is the cutoff frequency of the continuous contribution. Finally, ω k,i and S k,ia are the frequency and Huang-Rhys factor of mode k in pigment i / transition a. The total reorganization energy is then The discrete part of the SD was determined from a normal mode analysis [1] of the pigments in the crystal structure. The Huang-Rhys factor S k,j of each mode k was computed by projecting the gradient of the excited state j on the ground-state normal mode k [1]. The computed spectral densities are shown in Supplementary Fig. 1.
Spectra and rates. The exciton spectra of PCP were simulated employing the disordered exciton model coupled to the modified Redfield theory (mR) [2][3][4] for the description of exciton relaxation. The static disorder was modeled by allowing random Gaussian uncorrelated fluctuations of the site energies (parameters in Supplementary Table 1). For each realization of the disorder, the Hamiltonian in Eq. S1 is diagonalized, obtaining exciton energies E K and exciton coefficients C K j for every state |K . The spectrum is then calculated as: where ω K = E K / is the vertical excitation frequency of state K, τ K is the exciton lifetime as computed with mR theory, g K (t) represents the lineshape function of exciton state K, (S5) and C KKKK (ω) is the spectral density of exciton state K, which depends on the exciton coefficients and the spectral density of each state ia.
Computational details. All calculations are based on the crystal structure of main-form PCP (Ref. [5], PDB code: 1PPR). Hydrogen atoms were added using the LeAP module of AmberTools [6], considering all the titrable residues in their standard protonation state. The crystal structure was refined by relaxing the geometries of all the pigments through a QM/MM optimization in the ONIOM framework with electrostatic embedding [7], followed by normal modes and frequency calculations. At the optimized geometry, the gradient of the excited states was also computed in order to extract Huang-Rhys factors (see below). All QM/MM optimizations were performed at the B3LYP/6-31G(d) level of theory, keeping the MM part fixed. Note that B3LYP/6-31G(d) was found to give quite reasonable results for the geometry of peridinin (Per) [8].
The parameters for the protein were taken from the Amber ff99SB forcefield [9]. The charges of the cofactors were fitted to the HF/6-31G(d) electrostatic potential using the standard RESP procedure using AmberTools [6]. Selected Pers (614 and 624) were optimized in their S 2 excited state, using the same QM/MM ONIOM scheme with time-dependent DFT (TD-DFT). Additional optimizations were performed with the CAM-B3LYP functional.
The excitation energies of each pigment were computed at TD-DFT B3LYP/6-31G(d) level. The electronic couplings among all the excited states were computed with an approach based on the transition densities of the interacting pigments [10]. For all the site energy and coupling calculations, the protein environment was described through a polarizable QM/MM embedding (QM/MMPol) [10], in which the pigments of interest are described quantum mechanically, whereas the protein and the other cofactors are treated as a set of point charges and isotropic polarizabilities. The parameters for the MMPol forcefield were taken from Wang et al. [11], except for the MMPol charges of the pigments, which were fitted to the B3LYP/6-31G(d) electrostatic potential consistently with the polarization.

Supplementary Methods 2: Kinetic global analysis
The total 2DES signal recorded using laser 3 was global fitted using the kinetic model reported in Figure 5a in the main text. We are interested in the cascading signal from the excitation of the Pers, thus we selected the excitation frequency range from 1550 cm −1 to 18000 cm −1 for the analysis. The method is based on the global analysis of 2D electronic spectroscopy (2DES) data reported in Ref. [12]. In this work, we focus the analysis on the population dynamics using a suitable basis set of kinetic decay components. The objective of the analysis is to determine the time constants T e e of the energy transport processes e → e , and to get spectral information on where a particular process contributes to the 2DES maps.
Kinetic model. We employed a three states kinetic model for which we can build the population transport rate matrix as where e, e = S 2 , S 1 , Q y . Several constraints can be put into this matrix, mainly some transport channels can be switched off. The specific transport rate matrix can be specified as where zero elements stand for switched off channels. Using matrix K we can calculate the population Green's functions [13], organized in matrix G(t) and we can solve the Pauli master equation Each population Green's function G e e (t) expresses the probability of the population of state e to be transferred to state e at time t. We selected only the diagonal components G ee (t) to be used in the fitting procedure as they represent a linearly independent basis set for all the population Green's functions.
In order to account for the inhomogeneity of the population transfer channel S 2 → Q y the population Green's functions were averaged for 500 realizations of the kinetics problem using a lognormal distributed time constant T Q y S2 . The choice of the lognormal probability distribution shape is guided by the satisfactorily match with the computed rates histogram reported in Supplementary Global fitting. The experimental data is represented as a three dimensional array X ijk , where i and j are the indexes of emission and excitation frequencies and k is the index the population time t. We modeled the data as a sum of separated contributions expressed as the product of complex amplitude maps, stored as pages in a three dimensional array A, and kinetic decays stored as columns of matrix C(K). Therefore, the model function M can be defined as The problem is written in term of the unconstrained minimization of the least squares of the residuals, R(A, K) = X − M(A, K), obtaining The minimization over the complex amplitudes is analytically solved employing the variable projection algorithm [12,14]. In this way the residuals are specified as a function of matrix K only as where I is the identity matrix and C + is the Moore-Penrose pseudo inverse of matrix C. The simplified minimization problem becomes After solving the minimization problem, the complex amplitudes are reconstructed as The amplitudes of the components of the basis set of the kinetic model resulting from the global fit analysis of the pure absorptive dataset recorded with laser 3 are reported in Supplementary Fig. 6a-c. The most significant feature is represented in G S2S2 map, where, on the bottom, a negative signal represents the rising behavior of (17500,14850) cm −1 cross-peak due to S 2 → Q y EET. The process is simultaneous with the decay of S 2 ESA signal, recognized as a negative peak close to the diagonal, and with S 1 ESA formation, as the positive feature. The decaying behavior of S 1 ESA signal, due to EET towards Q y , is clearly represented by the negative signal in G S1S1 map. Finally, the positive cross-peak in G QyQy map witnesses the formation of the Q y population.

S5
The evolution of the Green's function of the basis set components and of the significant off-diagonal elements of the kinetic model is shown in Supplementary  Fig. 6d and e, respectively.
As direct consequence of the EET rates determined by the fitting analysis, >67% of initially prepared S 2 population is directly transferred to Q y . Thus, according to our results, only a minor fraction of the Per624 population is funneled to Q y via the S 1 state, as shown in Figure 5b (main text).

Supplementary Note 1: Calculated parameters
Site energies. The calculated site energies of Per and Chl a are reported in Supplementary Table 2. The calculated Q y site energies of the two Chl in PCP are not identical, with Chl602 being ∼170 cm −1 blue-shifted with respect to Chl601. The Q y site energies are ∼1750 cm −1 higher than the experimental values; however, the relative site energy difference between the two Chls is sensible, and in agreement with experimental steady-state measurements [15].
The site energies of Pers differ by more than 1000 cm −1 , due to the different protein environment. Surprisingly, Per624 is strongly red-shifted by the protein environment, whereas its pseudo-symmetric Per614 is among the less red-shifted. Per614 and 624 are equivalent by position and are the most interacting with Chl. However, their binding pocket is not the same. To verify that the redshift originates from the direct effect of the protein environment, we used the internal geometry of Per614, but in the pocket of Per624, obtaining exactly the same ∼1900 cm −1 solvatochromic shift (with respect to gas phase) observed for Per624. This demonstrates that the solvatochromic shift is not a result of the different ground-state geometry, but rather a direct electrostatic/polarization effect.
The differences between the binding pockets of the two Pers are summarized in Supplementary Fig. 2. As shown in the Figure, there are significant substitutions in the charged residues of the two pockets. In particular, there are various negatively charged residues on the allene side of Per624, which are substituted by neutral residues in Per614's pocket. On its lactone side, Per614 is hydrogen bonded to Glu101, which is substituted by Lys263 in Per624's binding pocket. The particular arrangement of charged residues around Per624 generates an electrostatic field along the Per backbone, as shown by the electrostatic potential maps of Supplementary Fig. 3 and Figure 2b in the main text. The other red-shifted Pers are 611/621 and 623. These red-shifted Pers are more likely to be excited by the laser tuned on the red tail of the Per absorption.
Exciton couplings. Supplementary Table 3-5 summarize the electronic couplings calculated in the protein environment among the pigments of PCP. The magnitude of the interactions shows that the two Chls are very weakly coupled. On the contrary, Pers are strongly coupled, with some couplings reaching values larger than 350 cm −1 .

S6
Finally, turning to the Per/Chl interactions, in Supplementary Table 5 we compared the couplings between the Pers' S 2 states and the Q x/y states of Chl. While the couplings to Q x states are always <100 cm −1 , the couplings to Q y reach values around 300 cm −1 , that is same magnitude as Per-Per interactions. The differences between S 2 /Q x and S 2 /Q y couplings can be ascribed to the lower dipole strength of the Q x states compared to the Q y . Notably, Pers 614 and 624, which was identified as the ones that strongly interact with the Chl of the same cluster, show the largest coupling to the Q y state.

Supplementary Note 2: Prediction of the excitonic spectrum and rate analysis
We used the exciton Hamiltonian presented in Supplementary Tables 2-5 and we calculated the absorption spectrum of PCP in the entire VIS region using modified Redfield theory. We shifted all the site energies of the Chl Q y and B transitions by −1750 cm −1 (maintaining the same relative differences as calculated with QM/MMPol) to compute the whole spectrum since our QM level of theory yields different errors on the excitation energies of different transitions. In the same way, we shifted the Q x site energies by −2100 cm −1 , and the S 2 energies by 3000 cm −1 . This gives a reasonable agreement with the experimental spectrum (see Figure 1 in the main text).
The mR simulations were also used to analyze S 2 → Q y (and S 2 → Q x ) EET rates as a function of the Per energy, based on 10 000 realizations of the static disorder. First of all we show in Supplementary Fig. 4a and b the relationship between the Per zero-phonon energy (i.e. the adiabatic energy E K − λ K ) and the dipole strength and delocalization length, the latter defined as inverse participation ratio (IPR). The IPR of a Per state K is calculated on the basis of exciton coefficients in a particular realization: This quantity ranges from a minimum of 1, for the states localized on only one Per, to a maximum of 8, that is the total number of Pers. At the red end of the S 2 band (ZPE < 17500 cm −1 ), the static disorder localizes the exciton on only one Per ( Supplementary Fig. 4b), and the dipole strength of the excited state resembles that of a single Per (Supplementary Fig. 4a). For this reason, we can conclude that the Per states are virtually localized when exciting at the red edge of the S 2 band. Figure 4c and d shows the relationship between zero-phonon Per energies and the S 2 → Q y and S 2 → Q x EET rates. When the states are localized on only one site (c K ia > 0.9), the rates can be assigned to one specific Per. Each Per transfers the excitation energy with a different rate depending on the coupling to the Q x/y state. Notably, S 2 → Q x EET rates reach at most 3 ps −1 , whereas S 2 → Q y EET rates are much faster, between 5 and 20 ps −1 , which corresponds to transfer times of 50-200 fs. The largest EET rates occur at the red-most edge of the S 2 band, where the Per states are localized and the S 2 /Q y energy gap is minimum.
The Per → Chl EET essentially occurs from a single Per donor, which may be Per624, 611(621) or 613(623) depending on which is the red-most Per in a particular configuration of static disorder. Different Pers transfer the excitation to the Q state with different rates, depending on the coupling with the acceptor state. The fastest transfer to Q y states occurs for Per624. The energy region probed in 2DES experiments corresponds to the very red edge of Per absorption, slightly above 17000 cm −1 . In this region, different Pers contribute to the EET toward the Q y state, but only Per624 can transfer energy to the Q y state of Chl with an EET time shorter than 100 fs, as observed in experiments; this Per was already proposed as the most coupled to Chl.

Supplementary Note 3: Bilinear time-frequency analysis
Signals from the S 2 /Q y cross-peaks were investigated using time-frequency methods which can provide simultaneously frequency and time resolution, unveiling the dynamics of the relevant beating components. In particular, four signals extracted from the 2DES rephasing maps obtained with laser 3 at excitation frequencies 17200, 17600, 17800, 18000 cm −1 and emission frequency 14850 cm −1 were analyzed with the smoothed-pseudo-Wigner-Ville distribution (SPWV) [16]. This time-frequency decomposition is derived from the Wigner-Ville distribution, two smoothing time windows, h(t) and g(t), are added in order to mitigate the presence of interference terms. SPWV is defined as (S16) where s(t) is the investigated signal. The time windows are chosen with Gaussian shape and standard deviations σ h = 75 fs and σ g = 100 fs. The resulting distributions are shown in Supplementary Fig. 7. Two kinds of contributions are identified: long lasting features with time constants of hundreds of femtoseconds ascribable to Per vibrational coherences and rapidly decaying features with time constants of tens of femtoseconds assigned to electronic coherences between Q y and S 2 states. The most striking difference between the two kinds of features is their shape in the time-frequency plane which is determined by the timefrequency uncertainty principle. Long lasting signals have a narrow frequency footprint, whereas rapidly decaying signals have a large frequency bandwidth. In the overlap between the two contributions, negative interference artifacts are present.
To clarify how the different signal contributions affect the SPWV distribution, a detailed model analysis was performed on the signal extracted at (17800, 14850) cm −1 . The signal was fitted with a model function composed by four oscillating components. The model function f (t 2 ) is specified as where a i are the complex amplitudes, τ i are the damping times, ω i are the oscillating frequencies and σ 4 is the standard deviation for the Gaussian damping employed for the electronic coherence. The parameters obtained from the fitting procedure are reported in Supplementary Table 6, whereas the plot of the model function, superimposed to the real data, is reported in Supplementary Fig. 8. The model fitting allowed us to calculate separately the SPWV distribution of each independent contribution to the signal. As shown in Supplementary Fig. 9, the contributions overlapping in the time-frequency plane were thus disentangled. The model fitting allowed us to calculate separately the SPWV distribution of each independent contribution to the signal and thus to disentangle the contributions overlapping in the time-frequency plane. The time-frequency analysis has also been performed to the above diagonal Q y /S 2 cross-peak. The interpretation of the results at those coordinates is however a bit more challenging because the overwhelming contribution of the vibrational modes of Per almost completely hide the presence of the 1900 cm −1 component. We report below Supplementary Fig. 10  Representation of the differences in the binding pockets of Per614 (a) and Per624 (b). Only the polar/charged residues that differ in the two pockets are shown, and they are colored according to the residue name. Negatively and positively charged residues are highlighted in red and blue, respectively.        Fourier transform (a) and bilinear time-frequency transform (b) of the oscillating part of the rephasing signal recorded with laser 3 and extracted at coordinates corresponding to the Qy/S2 above diagonal cross-peak.   , that match quite well the frequency of the three vibrational modes of Per more strongly coupled to the S0 → S2 transition. The intensity of Per vibrational modes increases as the laser spectrum is shifted to higher energies, i.e. S2 transition is progressively brought into resonance. Moreover, power spectrum recorded with laser 1 reveals several peaks below 1000 cm −1 due to Chl vibrational modes. S19  Amplitude of the three main oscillating contributions as a function of excitation and emission frequency (Fourier maps) of the 2DES rephasing (upper row) and non-rephasing (lower row) datasets. Experimental data (a) and simulated data (b) obtained with laser 1. Experimental and simulated rephasing maps of all the vibrational modes differ mainly in the strong cross-peak above diagonal.  Supplementary Figure 18: Fourier maps for laser 2. Same as Supplementary  Fig. 17 with laser 2. The intensity of the cross-peak above diagonal not predicted by simulations increases as the exciting pulses progressively bring S2 transition into resonance.  Supplementary Figure 19: Fourier maps for laser 3. Same as Supplementary  Fig. 17 with laser 3. The intensity of the cross-peak above diagonal not predicted by simulations increases as the exciting pulses progressively bring S2 transition into resonance.