Local protein solvation drives direct down-conversion in phycobiliprotein PC645 via incoherent vibronic transport

Samuel M. Blaua,1, Doran I. G. Bennetta,b,1, Christoph Kreisbecka, Gregory D. Scholesb,c, and Ala´ n Aspuru-Guzika,b,2

aDepartment of Chemistry and Chemical Biology, Harvard University, Cambridge, MA 02138; bBio-Inspired Solar Energy Program, Canadian Institute for Advanced Research, Toronto, ON M5G 1Z8, Canada; and cDepartment of Chemistry, Princeton University, Princeton, NJ 08544

Edited by Michael L. Klein, Temple University, Philadelphia, PA, and approved March 1, 2018 (received for review January 8, 2018)

The mechanisms controlling excitation energy transport (EET) in light-harvesting complexes remain controversial. Following the observation of long-lived beats in 2D electronic spectroscopy of PC645, vibronic coherence, the delocalization of excited states between pigments supported by a resonant vibration, has been proposed to enable direct excitation transport from the highestenergy to the lowest-energy pigments, bypassing a collection of intermediate states. Here, we instead show that for phycobiliprotein PC645 an incoherent vibronic transport mechanism is at play. We quantify the solvation dynamics of individual pigments using ab initio quantum mechanics/molecular mechanics (QM/MM) nuclear dynamics. Our atomistic spectral densities reproduce experimental observations ranging from absorption and fluorescence spectra to the timescales and selectivity of down-conversion observed in transient absorption measurements. We construct a general model for vibronic dimers and establish the parameter regimes of coherent and incoherent vibronic transport. We demonstrate that direct down-conversion in PC645 proceeds incoherently, enhanced by large reorganization energies and a broad collection of high-frequency vibrations. We suggest that a similar incoherent mechanism is appropriate across phycobiliproteins and represents a potential design principle for nanoscale control of EET.

occurs by transient delocalization; otherwise incoherent excitation transport occurs by discrete hopping. Persistent beat signals in nonlinear spectroscopic measurements of multiple LHCs have been interpreted as evidence for coherent vibronic transport (7–16). Phycobiliprotein PC645, in particular, has been studied extensively and has been suggested to leverage a coherent vibronic enhancement in the down-conversion process central to its physiological function (10). Here, we instead show that
down-conversion in PC645 occurs via an incoherent vibronic mechanism.
In vivo, PC645 absorbs high-energy photons and downconverts the excitation to be resonant with the reaction center, where charge separation occurs. Transient absorption measurements on isolated PC645 reveal the process of down-conversion: Photoexcitation of the high-energy core dihydrobiliverdins (DBVs) (Fig. 1A, blue and green) is followed by transfer to the low-energy phycocyanobilin 82s (PCB82s) (Fig. 1A, red and cyan), skipping the energetically intermediate mesobiliverdins (MBVs) (Fig. 1A, brown and magenta) (17). Direct downconversion cannot be explained by interpigment electronic couplings (18), suggesting that the vibrational dynamics arising from the local protein environments determine the pathways of EET.

photosynthesis | light harvesting | excitation energy transfer |
quantum coherence | molecular dynamics
U nderstanding the solvation dynamics of multiple pigments in a heterogeneous medium is essential for being able to control excited-state dynamics in artificial materials. Natural photosynthesis demonstrates this capability in light-harvesting complexes (LHCs), which exert nanoscale control over excitation transport between pigments via the protein environment. Currently, there is no spectroscopic technique that can disentangle excitation dynamics in pigment aggregates from the many timescales of the vibrational bath, which include intramolecular pigment vibrations, fluctuations of neighboring amino acids, reorientation of proximate biological water, and reorganization of the bulk water solvent. Thus, quantitatively connecting the atomic-scale motions of the pigment–protein environment to the dynamics of excitation transport in LHCs remains a grand challenge at the interface of materials science and chemical physics.
Are there simple design principles that connect the complex atomic structure of pigment environments to the regulation of excitation energy transfer (EET) pathways in LHCs (1–5)? As demonstrated by the environmentally assisted quantum transport (ENAQT) mechanism (6), thermally occupied lowfrequency vibrations can enhance transport between energetically detuned pigments and represent the simplest design principle connecting atomistic dynamics and excitation transport. The presence of a high-frequency vibration resonant with the excitation energy gap between pigments can also enhance transport. If the vibration supports coherence between the pigment excited states, i.e., vibronic coherence, then excitation transport

Significance
To harness sunlight for growth, plants and algae rapidly convey absorbed excitation energy from antennae pigments to a reaction center, where the excitations convert to chemical energy. In deep water environments, cryptophyte algae survive by using the molecular motion of phycobiliprotein antennae to funnel excitations to low-energy pigments. The mechanism of down-conversion in phycobiliproteins remains controversial: Could specific vibrations resonant with pigment energy gaps support coherence between excited states and thereby enhance the rate of transport by transient delocalization? Here, we demonstrate that down-conversion in a specific phycobiliprotein, PC645, is both incoherent and enhanced by a broad range of high-frequency vibrations. We suggest that a similar incoherent mechanism is appropriate across phycobiliproteins.
Author contributions: S.M.B., D.I.G.B., and C.K. designed research; S.M.B., D.I.G.B., and C.K. performed research; S.M.B. and D.I.G.B. contributed new reagents/analytic tools; S.M.B., D.I.G.B., and C.K. analyzed data; and S.M.B., D.I.G.B., C.K., G.D.S., and A.A.-G. wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
This open access article is distributed under Creative Commons AttributionNonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND). 1 S.M.B. and D.I.G.B. contributed equally to this work. 2 To whom correspondence should be addressed. Email: aspuru@chemistry.harvard.edu.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10. 1073/pnas.1800370115/-/DCSupplemental.
Published online March 27, 2018.

E3342–E3350 | PNAS | vol. 115 | no. 15

www.pnas.org/cgi/doi/10.1073/pnas.1800370115

A
17,000

PCB82C

MBVB

B

CHEMISTRY

BIOPHYSICS AND COMPUTATIONAL BIOLOGY

16,500

DBVC

Energy [cm-1]

16,000
15,500
15,000
D
6,000
J(ω) 4,000 [cm-1]
2,000

PCB158C MBVA

PCB158D
DBVD
PCB82D
Construct correlation
function FFT

Forces from DFT
Forces from MD

Transition density snapshot

C
18,000

TDDFT excited state

t = 10 ps

Energy gap [cm-1]

16,000

14,000

0 500 1,000 1,500 ω [cm-1]

10 20 30 40 Time [ps]

Fig. 1. PC645 and computational methodology for spectral density construction. (A) Structure and reorganized energy levels of PC645, where we have defined a bilin color scheme that remains consistent throughout. Gray lines in the energy-level diagram are drawn schematically to represent the exciton eigenstates of the core DBVs (blue and green) due to the significant excitonic coupling of 319.4 cm−1, more than three times larger than any other coupling in PC645. (B) For each bilin, we simulate 50 ps of mixed QM/MM nuclear dynamics in which that bilin is treated quantum mechanically, while the remaining protein is treated with a classical force field. (C) We construct energy-gap trajectories by calculating full TDDFT excited states for each bilin on geometries extracted at 2- fs intervals from our QM/MM trajectories. (D) We construct an energy-gap correlation function for each bilin and Fourier transform to obtain eight unique bilin spectral densities.

Quantifying the ultrafast vibrational dynamics of specific chromophores in LHCs remains a challenge for both theory and experiment. Using a combination of linear and nonlinear spectroscopic data, one can parameterize a model that describes the local vibrational environments of each pigment (19). However, multiple distinct vibrational environments quickly lead to an intractable number of free parameters. Atomistic nuclear dynamics simulations have the potential to bypass the inverse problem (20–23) but have not previously managed to quantitatively reproduce spectroscopic signals, thus restricting their ability to inform on biological function.
Here, we successfully connect unique pigment–protein vibrational environments to spectroscopic signals by incorporating pigment forces calculated with ground-state density functional theory (DFT). We extract bilin reorganization energies with a wider spread and larger average value than previously expected (10, 18, 24), allowing for quantitative reproduction of linear spectra in the absence of free parameters. Additionally, our EET simulations yield transfer pathways and timescales in good agreement with experimental transient absorption measurements. To address the underlying mechanism, we construct a representative vibronic dimer and establish the regimes of coherent and incoherent vibronic transport. We find that PC645 is far into the

incoherent parameter regime and support this conclusion by demonstrating that down-conversion in PC645 does not depend on the frequency, or even the presence, of the resonant vibrational mode previously suggested to support vibronic coherence. Our results demonstrate that transport between bilins in PC645 occurs incoherently and that Fo¨rster spectral overlaps tuned by local bilin–protein environments control EET and enable direct down-conversion. We further suggest that vibronic transport in phycobiliproteins, generally, is incoherent and enhanced by the presence of a broad collection of high-frequency, intramolecular vibrations.
Results and Discussion
Ab Initio Simulations of Protein Solvation Reproduce Linear Spectra. To understand the role of the local pigment environments in controlling EET, we first characterize the atomistic origin of the bilin solvation dynamics. The excitation energy of a pigment bound in a protein pocket, given as the energy difference between the ground and first excited state, fluctuates due to the nuclear motions of the chromophore, the surrounding protein residues, and the proximate water. Simulating excitation energy fluctuations requires performing nuclear dynamics that sample the ground-state potential energy surface (PES) and calculating

Blau et al.

PNAS | vol. 115 | no. 15 | E3343

the excitation energy at regular intervals. These fluctuations can be characterized by a spectral density that describes coupling of the electronic excited state to a continuous distribution of vibrational modes (25). Given the significant computational cost of obtaining nuclear forces quantum mechanically (QM), previous calculations have always used classical or semiempirical molecular mechanics (MM) force fields to propagate nuclear dynamics for entire LHCs (20–23, 26). However, in many cases, the MM PES sampled during nuclear dynamics differs substantially from the QM PES on which excitation energies are defined, causing the calculated spectral densities to report incorrect vibrational frequencies and coupling amplitudes (27).
We construct spectral densities using nuclear dynamics calculated on an ab initio QM/MM PES that combines bilin nuclear gradients calculated using DFT with an MM force field to treat the surrounding protein environment (Fig. 1B), thereby resolving the mismatch between nuclear dynamics trajectories and excited-state calculations. Fig. 1 provides an overview of our procedure for constructing spectral densities. For each bilin, we construct energy-gap trajectories (Fig. 1C) by calculating excitation energies with time-dependent density functional

theory (TDDFT) on 20,000 geometries extracted at 2-fs intervals from our 40-ps QM/MM production runs. Fourier transforms of the two-time correlation functions of the energy-gap trajectories define unique spectral densities that characterize individual bilin environments (Fig. 1D). We note that while TDDFT often incorrectly predicts absolute excitation energies, cancellation of error results in good descriptions of the curvature of the potential energy surface that we depend on here (28–31).
Spectroscopic signals and EET dynamics depend intimately on the solvation dynamics of individual chromophores. The solvation capacity of the local vibrational environment is quantified by the total reorganization energy (λ) of the spectral density which measures the energy dissipated as the chromophore relaxes following excitation. Previous studies have assumed that all bilins have identical spectral densities and assigned a reorganization energy of 260 cm−1 (24), 314 cm−1 (10), or 480 cm−1 (18). Our spectral densities, shown in gray in the subpanels of Fig. 2, reveal larger reorganization energies ( λ = 909 cm−1) and significant variations between different bilins. The spread in our reorganization energies (λ = 627–1,414 cm−1) is consistent with the presence of three chemically distinct bilins with different conjugation

8,000
J( ) [cm-1]

DBVC,

= 748 cm-1

4,000

Experiment HEOM

0 500 1,000

8,000
J( ) [cm-1]

DBVD,

= 672 cm-1

4,000

1,500 [cm-1]

Fluorescence

0 500 1,000

8,000
J( ) [cm-1]

MBVA,

= 1414 cm-1

4,000

0

1,500

[cm-1] 8,000
J( ) [cm-1]

PCB158C,

14,000
= 974 cm-1

4,000

Absorption

16,000
8,000
J( ) [cm-1]

PCB82C,

18,000
= 677 cm-1

4,000

ω [cm-1]

0 500 1,000

8,000
J( ) [cm-1]

MBVB,

= 1104 cm-1

4,000

1,500

0 [cm-1]
8,000
J( ) [cm-1]

500 1,000 PCB158D, = 1026 cm-1

1,500

0 [cm-1]
8,000
J( ) [cm-1]

500 1,000 PCB82D, = 627 cm-1

4,000

4,000

1,500 [cm-1]

000

500

1,000

1,500 [cm-1]

500

1,000

1,500 [cm-1]

500

1,000

1,500 [cm-1]

Fig. 2. Comparison between absorption and fluorescence simulated with atomistic spectral densities and experimental spectra. Unabridged (gray) and abridged (colored) bilin spectral densities are shown in the subpanels, where the latter are used for spectroscopic simulations. Each spectral density panel is labeled with the bilin name and the calculated reorganization energy λ. The main panel compares the calculated absorption and fluorescence spectra (solid black lines) with experimental spectra (open gray circles).

E3344 | www.pnas.org/cgi/doi/10.1073/pnas.1800370115

Blau et al.

CHEMISTRY

lengths and either one or two covalent protein linkages. The variability of the reorganization energies between chemically identical bilins bound in distinct protein environments (e.g., MBVA and MBVB) demonstrates the importance of the protein scaffold in determining the local solvation dynamics. On the other hand, all eight bilin spectral densities show a peak near 1,650 cm−1 which is consistent with the assignment of a long-lived 1,580cm−1 mode as intramolecular C = C or C = N vibrations in broad-band transient absorption measurements (10).
Incorporating the ab initio QM/MM spectral densities into absorption and fluorescence simulations (solid black lines) results in excellent agreement with experimental spectra (open gray circles), as shown in Fig. 2. We note that the negative features observed in simulated spectra are unphysical but do not impact the quality of the lineshapes (SI Appendix, Fig. S4). To account for both the large reorganization energies and the multiple timescales of vibrational relaxation encoded in the structure of our spectral densities, we use numerically exact hierarchical equations of motion (HEOM) (32, 33), as implemented in QMaster (34). HEOM has been shown to yield realistic simulations of exciton dynamics in LHCs (35–37), but its computational complexity limits our spectral densities to including a total of 24 Drude–Lorentz peaks for full-system simulations, each representing a distributions of bath modes. We construct four classes of abridged spectral densities for each bilin, starting with class 1 (eight or nine peaks per bilin) and successively coarse graining to incorporate fewer peaks until we reach class 4 (two peaks per bilin). Peak parameters for all spectral densities are given in SI Appendix, Tables S1–S4. Abridged spectral densities used for spectroscopic calculations are shown with our consistent bilin color scheme in the subpanels of Fig. 2. We justify our abridged spectral densities by comparing simulated monomer absorption and fluorescence spectra, as described in SI Appendix, section 2. SI Appendix, section 2 discusses the system Hamiltonian while SI Appendix, section 3 covers ab initio transition dipole moments, details of the fluorescence simulations, and inhomogeneous broadening.
We validate our spectral densities by comparing them to two features of the experimental spectra: the Stokes shift, defined as the frequency difference between the highest-energy fluorescence peak and the lowest-energy absorption peak, and the overall width of the absorption spectrum. Together, these observables are quite sensitive to the distribution of reorganization energies between the different bilins of PC645. Using the lowest-energy bilin spectral density for all pigments, such as would be extracted from a fluorescence line-narrowing experiment, results in an absorption spectrum with a full-width at half maximum (FWHM) that is 25% too small (SI Appendix, Fig. S5). The absorption spectrum contracts in this case because we underestimate the reorganization energies of the PCB158s and MBVs. If we instead use the average of the eight bilin spectral densities for all pigments, we obtain an absorption spectrum with a Stokes shift that is 55% too large (SI Appendix, Fig. S6). The Stokes shift expands because we substantially overestimate the reorganization energy of the PCB82s that dominate the fluorescence and low-energy absorption features. Thus, 18% error between the experimental and simulated Stokes shift and very good agreement with the overall absorption linewidth in the absence of free parameters represent a strong validation of the solvation dynamics encoded in our atomistic spectral densities.
Protein Solvation Drives Down-Conversion. The spectral density connects the atomistic dynamics of the bilin vibrational environment to both spectroscopic lineshapes and EET pertinent to light harvesting. Having validated our atomistic QM/MM spectral densities against absorption and fluorescence spectra, we now confirm that they also describe the experimental observables associated with direct down-conversion in PC645. Experimen-

tal transient absorption measurements, combined with global kinetic analysis, indicate that initial excitation of the DBV core is followed by direct transport to the low-energy PCB82 pigments with a rate of 1.7 ps−1 (10, 17). Neither the rates nor selectivity of direct down-conversion have been successfully reproduced using previously estimated unstructured spectral densities (24). Here, we show that numerically exact HEOM simulations combined with our spectral densities predict pathways and rates of down-conversion that reproduce experimental results.
The large reorganization energies of our QM/MM spectral densities are essential for correctly predicting the selectivity of transport from the DBV core to the low-energy PCB82s. Simulated EET pathways following photoexcitation of the DBV core are shown in Fig. 3 A and B. Arrow thickness in Fig. 3 is proportional to exciton flux, the net amount of excitation that is transferred from the DBVs to an acceptor pigment in a 1.0-ps interval following photoexcitation. We perform exciton flux calculations using class 2 spectral densities for the MBVs and class 4 spectral densities for the remaining pigments (SI Appendix, section 2). In Fig. 3A we rescale the spectral densities to match the average magnitude of reorganization energy used in previous simulations ( λ = 260 cm−1, labeled 30%). In the presence of a smaller average reorganization energy, the more weakly solvated DBV core primarily transports excitation to the energetically adjacent MBVB, qualitatively reproducing previous results (24). EET simulations using the full reorganization energy (Fig. 3B, labeled 100%), however, show enhanced direct transfer to the low-energy PCB82s, in reasonable agreement with global kinetic analysis of transient absorption measurements.
In addition to influencing selectivity, the large reorganization energies also increase the rate of direct down-conversion. HEOM is a non-Markovian theory, and as a result, the rates of transport vary as a function of time in response to changes in the vibrational energy distribution. We extract a best-fit rate of down-conversion from HEOM simulations using class 1 spectral densities with a four-site model containing only the core DBVs and the low-energy PCB82s (SI Appendix, section 4). Class 1 spectral densities (e.g., SI Appendix, Fig. S1A) explicitly incorporate the high-frequency mode (∼1,650 cm−1) previously assigned to an intramolecular bilin vibration. The rescaled spectral densities (Fig. 3A, 30%) result in a transport rate of 0.8 ps−1, substantially slower than the experimentally observed 1.7 ps−1. However, using the full reorganization energy of our spectral densities (Fig. 3B, 100%), we find the simulated rate of transport to be 1.6 ps−1, within 8% of the experimental value.
Down-Conversion Occurs via an Incoherent Vibronic Mechanism. We have demonstrated that HEOM simulations using our atomistic spectral densities reproduce both the linear spectra and the rate of direct down-conversion in PC645. However, the question remains: What are the mechanistic principles connecting bilin vibrational environments to regulation of EET pathways? In particular, the relative importance of shortlived, low-frequency, intermolecular motions vs. long-lived, highfrequency, intramolecular vibrations remains controversial (1, 15). In this section, we assess how the mechanism of EET in PC645 arises from the interplay of electronic couplings between pigments with structured spectral densities. First, we introduce the simplest representative model system that captures the essential features of vibronic transport in the presence of a long-lived high-frequency vibration and distinguish between the coherent and incoherent transport regimes. Second, we demonstrate that PC645 experiences an incoherent vibronic mechanism and that the long-lived near-resonant vibration assigned to support vibronic coherence is not essential for efficient transport. Finally, we show that generalized Fo¨rster theory (38–41) explains the mechanism of transport in PC645, consistent with recent work on other PPC aggregates (34, 42–45), and that local protein

BIOPHYSICS AND COMPUTATIONAL BIOLOGY

Blau et al.

PNAS | vol. 115 | no. 15 | E3345

A 30% of reorganization energy

PCB82C

MBVB

B 100% of reorganization energy

PCB82C

MBVB

Arrow thickness = exciton flux

MBVA

PCB82D

DBV to PCB82 rate: 0.8 ps-1

MBVA

PCB82D

DBV to PCB82 rate: 1.6 ps-1

Fig. 3. EET pathways and down-conversion rates simulated with HEOM. Arrow thickness is proportional to the exciton flux leaving the DBV core integrated over 1 ps after an initial photoexcitation with (A) 30% or (B) 100% of calculated reorganization energy. The corresponding rates from the DBV core to both PCB82s are reported below each panel.

solvation enables direct down-conversion by controlling Fo¨rster
spectral overlaps.
The transition from coherent to incoherent vibronic transport
is governed by the ratio of the vibronic coupling (Vvib) to the rapidly relaxing component of the pigment reorganization energy (λdeph). Fig. 4A depicts an energetically detuned dimer with a long-lived resonant high-frequency vibration on the low-energy acceptor pigment (further details in SI Appendix, section 5). In
brief, we define the basis as a direct product of three indexes: the electronic state of the acceptor (|ga , |ea ), the electronic state of the donor (|gd , |ed ), and the vibrational state of the acceptor on either the ground- or excited-state harmonic oscillator (|0g , |0e , |1e ). Assuming the energy gap between the donor (|ga, ed, 0g ) and the acceptor (|ea, gd, 0e ) is large compared with λdeph, rapid excitation transport cannot occur (Fig. 4B, purple line) unless a vibrationally excited acceptor state (|ea, gd, 1e ) is generated by a bridging high-frequency mode (Fig. 4B, black line). Vibronic
coupling (Vvib, Eq. 1) between the nearly degenerate states relevant for vibronic enhancement (Fig. 4B, gray area) is described
by the product of the electronic coupling (V) with the Frank–
Condon factor (46):

Vvib = V · 1e|0g ∼ V ·

λvib . Ωvib

[1]

The vibronic coupling facilitates delocalization between the
donor and the vibrationally excited acceptor. In contrast, the
component of the pigment reorganization energy that relaxes
faster than the timescale of transport (λdeph) drives excitations to localize on individual pigments (47–49). Therefore, when
Vvib λdeph, transport is coherent (i.e., aided by delocalization) and initial dynamics rise and oscillate (Fig. 4C, blue/green lines) on a timescale determined by Vvib (Fig. 4C, black dashed line). However, when λdeph ≈ Vvib (Fig. 4C, orange line), dynamic localization inhibits the initial coherent rise. Finally, when λdeph Vvib (Fig. 4C, red line), complete localization results in incoherent hopping that can be described by Fo¨rster theory,

Ka←d

=

|V|2 π2

Re

∞
Aa(t)F∗d (t)dt ,
0

[2]

where Aa(t) is the absorption lineshape of the acceptor and Fd(t) is the fluorescence lineshape of the donor. In the incoherent regime (λdeph Vvib) HEOM and Fo¨rster theory predict identical transport rates (Fig. 4B, red dashed line) and dynamics (Fig. 4C, gray dashed line), assuming vibrational relaxation is not too slow. Thus, the presence of a resonant, intramolecular vibration can enhance incoherent transport by providing a vibronic sideband that increases the absorption/fluorescence overlap (14, 46) (Fig. 4C). We note that the general conclusions drawn from this model hold whether the vibration is placed on the acceptor, the donor, or both (SI Appendix, section 5).
We find that direct down-conversion in PC645 occurs in the incoherent regime of vibronic transport. Limited delocalization between the DBVs due to their strong electronic coupling and the many timescales of vibrational relaxation complicates the assignment of λdeph. The short-timescale (<200 fs) dynamic Stokes shift of the DBV core provides a lower bound on λdeph of 125 cm−1, as described in SI Appendix, section 6. The electronic coupling between the delocalized DBV state and the PCB82 pigments is estimated to be 50 cm−1. Combined with the Frank– Condon factor of ∼0.25 for the high-frequency vibration, the pertinent vibronic coupling is 12 cm−1. We, therefore, find that the ratio of λdeph to the vibronic coupling in PC645 has a lower bound of 10, indicating the dominance of an incoherent mechanism. Thus, since the energy gap between the DBV fluorescence and PCB82 absorption maxima is larger than their peak widths (SI Appendix, section 6), we find that direct down-conversion in PC645 occurs in the incoherent vibronic regime. Given the order of magnitude difference between the lowest estimate of λdeph and the vibronic coupling, our assignment of incoherent vibronic transport is robust to even substantial variation in Hamiltonian parameters (SI Appendix, section 6).
The presence of long-lived oscillations in nonlinear spectroscopic measurements led to the suggestion of functionally

E3346 | www.pnas.org/cgi/doi/10.1073/pnas.1800370115

Blau et al.

ABC

CHEMISTRY

BIOPHYSICS AND COMPUTATIONAL BIOLOGY

D

Fig. 4. Regimes of coherent and incoherent vibronic transport. (A) Energy levels of a detuned model dimer with a vibration on the acceptor pigment that is resonant with the site energy gap. The yellow shading represents the possible delocalization between the donor and the vibrationally excited acceptor states. (B) Rate of transport from donor to acceptor as a function of the reorganization energy of the low-frequency background (λdeph), where ED–EA = 250 cm−1, V = 2.75 cm−1, λvib = 24.7 cm−1, Ωvib = 274.7 cm−1, and γvib = 0.125 cm−1. For the model dimer, our calculations use a vibrational frequency much lower than that found in PC645 to ensure that no negative populations are observed in the coherent regime (we check the validity of PC645 calculations in SI Appendix, Fig. S4). Calculated rates of transport are shown for HEOM with a high-frequency vibration (HEOM with vib, solid black line), HEOM without vibration (HEOM no vib, solid purple line), and Fo¨ rster with vibration (Fo¨ rster with vib, dashed red line). The solid gray area represents the vibronic enhancement arising from the presence of the high-frequency vibration. The green, orange, and red boxes represent the population dynamics shown in C. (C) The acceptor population dynamics following the donor excitation are plotted for λdeph = 0.01, 0.1, 1, and 10 in blue, green, orange, and red, respectively. The first oscillation of purely coherent dynamics between the donor and vibrationally excited acceptor is plotted as a dashed black line. The population dynamics calculated using Fo¨ rster theory for λdeph = 10 are plotted as a dashed gray line. (D) Mechanism of incoherent vibronic enhancement of transport, in which a nearly resonant vibration can generate a vibronic sideband in the absorption (black) that can enhance overlap with the fluorescence. Absorption and fluorescence lineshapes in the absence of a high-frequency vibration are shown in purple.

relevant vibronic coherence in PC645 (10). However, recent work has demonstrated that oscillatory vibronic signatures in 2D electronic spectra do not necessarily imply a role of vibronic coherence in excitation energy transfer (50). We directly probe the role of vibronic coherence in PC645 by manipulating the high-frequency vibration previously assigned to support vibronic signatures and demonstrate that it is not required for efficient transport. Coherent vibronic transport is known to exhibit a sharp resonance condition as a function of the vibrational frequency (7, 10). The resonance condition arises because vibronic delocalization between the donor and the vibrationally excited acceptor occurs only when the energy gap is smaller than, or the same order of magnitude as, the vibronic coupling (i.e., Vvib ∼ 12 cm−1). In the case of PC645, partial delocalization in the DBVs means that there are two possible resonance conditions to consider, depending on whether the DBV core eigenstates or pigment states act as the donor. To explore the possible resonance conditions, we examine population dynamics calculated in a four-site model containing only the DBVs and PCB82s. We use class 1 spectral densities with four possible positions of the high-frequency mode (Fig. 5A): (i) the original positions (∼1,650 cm−1), (ii) the the energy difference between the DBVD and PCB82C sites, (iii) the energy difference between the low energy DBV exciton and the PCB82C site, and (iv) the complete removal of the high-frequency peak. We observe minimal differences in the resulting population dynamics (Fig. 5B). An additional exhaustive scan of the peak position between 1,100 cm−1 and 1,800 cm−1 shows no greater variations in popula-

tion dynamics. We have therefore demonstrated the absence of a sharp resonance condition, further supporting our assignment of incoherent vibronic transport.
Consistent with the incoherent transport regime, Fo¨rster theory captures the dominant contributions to the rate of downconversion in PC645. To describe incoherent transport in the presence of strong coupling between DBVs, we use a generalization of Fo¨rster theory (38) that treats excitations within the DBV core as delocalized but assumes transport out of the core can be described as an incoherent hop. Due to rapid exciton relaxation within the DBV core, the rate of transport out of the core is determined by the overlap between the lowestenergy DBV exciton fluorescence and the absorption spectra of the remaining pigments. We simulate absorption spectra for the non-DBV pigments using Kubo lineshapes, which are exact for local excitations. We simulate the low-energy DBV exciton fluorescence using an extension of Kubo theory (38, 42, 43), which neglects the localization of excitons resulting from vibrational fluctuations (SI Appendix, Fig. S10). Generalized Fo¨rster theory predicts a rate of 1.4 ps−1, within 15% of the HEOM result (Fig. 5C).
Using generalized Fo¨rster theory, we can explain the relative mechanistic role of low-frequency, intermolecular vibrations and high-frequency, intramolecular vibrations in controlling direct down-conversion in PC645. To understand the impact of bilin vibrational environments on transport rates and pathways, we examine the absorption/fluorescence overlap between the DBV core (Fig. 5 D and E, gray lines) and either PCB82C (Fig. 5E,

Blau et al.

PNAS | vol. 115 | no. 15 | E3347

J(ω) [cm-1]

A
3,000
2,000

Actual peak position Site basis resonance Eigenstate resonance

1,000

B
1.0
0.5

PCB82 population

Eigenstate resonance No peak
Site basis resonance Actual peak position

Energy [cm-1]

0 500
C 17,000
16,500 16,000 15,500 15,000

1,000 ω [cm-1]

1,500

DBV to PCB rate
rster 1.4 ps-1
HEOM 1.6 ps-1
Exp. 1.7 ps-1

2,000

0

D
MBVB abs. Overlap V2
DBV fl.

E
PCB82C abs. Overlap V2
DBV fl.

0.5 Time [ps]

1

14,500

14,500 15,000 15,500 16,000 16,500 17,000 Energy [cm-1]

Fig. 5. Direct down-conversion in PC645 is driven by incoherent vibronic transport. (A) Four different PCB82C spectral densities that we use to examine the impact of the high-frequency vibration on direct down-conversion. While cyan represents the case where the high-frequency vibration has been removed
entirely, it is also implicitly added to the following three cases in which the high-frequency vibration is shown by itself for clarity: the vibration in its original position in our QM/MM spectral densities (black), the vibration shifted to be in resonance with the DBVD – PCB82C site energy difference (1,550 cm−1, dashed gray line), and the vibration shifted to be in resonance with the low-energy DBV exciton – PCB82C energy difference (1,220 cm−1, black open circles). Note that all four sites present in our population dynamics simulations use their individual spectral densities, but the changes to the high-frequency peak are
equivalently applied to the other three sites. (B) Population dynamics of a four-site system containing the DBVs and PCB82s. Lines representing the sum PCB82 population are labeled by their representation in A. (C) An energetic depiction of down-conversion where site energies are fully reorganized and arrows represent the combined flux from the DBV core to the MBVs or PCB82s. The flux arrow from DBVs to PCB82s is labeled with the effective transport
rate obtained from theory and experiment. Note that excitons are drawn schematically to represent DBV delocalization. (D) Kubo monomer absorption for MBVB (magenta), Kubo DBV core fluorescence (gray), and absorption/fluorescence overlap weighted by the coupling to the lowest-energy DBV exciton (outlined magenta). (E) Kubo monomer absorption for PCB82C (cyan), Kubo DBV core fluorescence (gray), and absorption/fluorescence overlap weighted by the coupling to the lowest-energy DBV exciton (outlined cyan).

cyan) or MBVB (Fig. 5D, magenta). The Stokes shift of the DBV core (∼280 cm−1) enhances the energetic overlap with the lower-energy pigments and dramatically increases the overall rates of transport. Despite almost perfect alignment of the DBV fluorescence with the MBV absorption spectra, excitation is preferentially transported to PCB82s due to the presence of broad vibronic sidebands. In contrast to the narrow resonance condition associated with coherent vibronic transport (width on the order of |V 1e|0g |), the existence of wide absorption and fluorescence features allows for incoherent vibronic enhancement with a broad resonance condition (width on the order of the homogeneous linewidth). Finally, we note that given the

importance of the Stokes shift in determining the rate and pathways of excitation transport, the demonstration of modified local solvation dynamics for chemically identical bilins within PC645 represents a potential design principle for nanoscale control of EET.
Concluding Remarks
Here, we have demonstrated how specific features of the pigment vibrational environments in PC645 exert nanoscale control of EET pathways and enable direct down-conversion from the DBV core to the low-energy PCB82s. We accessed the atomistic origin of local solvation for each bilin, using ab initio QM/MM nuclear

E3348 | www.pnas.org/cgi/doi/10.1073/pnas.1800370115

Blau et al.

CHEMISTRY

dynamics. In contrast with the previous assumption of identical baths, our simulations identified substantial disparities between bilin vibrational environments which were essential to reproducing the experimental absorption and fluorescence spectra. Further, we have determined that down-conversion in PC645 proceeds via an incoherent vibronic transport mechanism where (i) excitations are localized on individual bilins (except for the DBV core) and transport occurs via incoherent hops and (ii) direct down-conversion is enhanced by large reorganization energies and the presence of a wide collection of high-frequency vibrations that induce broad vibronic sidebands. Given similar bilin flexibility and weak electronic couplings generally reported for biliproteins, our findings suggest that vibronic transport in cryptophyte algal antenna complexes proceeds incoherently.
The incoherent vibronic mechanism assigned here to PC645 is far more robust to imperfections than its coherent counterpart and could act as a blueprint for the design of artificial excitonic materials. Recent advances in biomimetic light-harvesting technologies have demonstrated novel architectures, e.g., metalorganic frameworks (MOFs) (51) and DNA origami (52), that can precisely place pigments within a hierarchically organized assembly. To fully realize nanoscale control of EET, we must iteratively simulate and design the local vibrational environments of pigment-scaffold architectures. However, the complexity of our current procedure is infeasible for high-throughput application, pointing toward the need for new and more efficient computational strategies and approximations.
Materials and Methods
We perform eight QM/MM nuclear dynamics simulations where the forces
on one bilin are constructed from ground-state DFT calculations for each
trajectory. We begin each QM/MM trajectory with 10 ps of equilibration

followed by 40 ps of production run. We construct energy gap trajectories for each bilin, using TDDFT calculations on geometries sampled every 2 fs, for a total of 20,000 geometries per bilin. We use the energy gap trajectories to construct two-time correlation functions and then Fourier transform to obtain the spectral density for each bilin. Overall, the eight 50-ps QM/MM trajectories took over 9 mo to run and cost more than 2 million central processing unit (CPU) hours. The details of the initial geometry preparation, separation of the quantum from classical regions, excited-state energy calculations, and the construction of the spectral densities are described in SI Appendix, section 1.
Exciton dynamics simulations were performed with the QMaster software
package (34) that provides a high-performance implementation of HEOM,
which runs flexibly on both graphics processing unit (GPU) and CPU architec-
tures (13, 34, 53). HEOM fluorescence calculations used a development ver-
sion of QMaster. All HEOM results presented were run at a hierarchy depth
of six except for rate calculations, which were run at a hierarchy depth of
five (SI Appendix, section 4).
ACKNOWLEDGMENTS. We gratefully acknowledge detailed discussions with and extensive commentary on the manuscript by Kapil Amarnath, as well as helpful conversations with Jacob Dean. We thank Thomas Markovich for help obtaining correlation functions and spectral densities. We acknowledge support from the Center for Excitonics, an Energy Frontier Research Center funded by the US Department of Energy (DOE), Office of Science and Office of Basic Energy Sciences, under Award DE-SC0001088. S.M.B. acknowledges support from the US DOE through the Computational Sciences Graduate Fellowship. D.I.G.B. and A.A.-G. acknowledge the John Templeton Foundation (Grant 60469). D.I.G.B., A.A.-G., and G.D.S. acknowledge the Canadian Institute for Advanced Research for support through the Bio-Inspired Solar Energy program. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the US DOE under Contract DE-AC02-05CH11231. We thank Nvidia for support via the Harvard CUDA Center of Excellence. This research used computational time on the Odyssey cluster, supported by the Faculty of Arts and Sciences Research Computing Group at Harvard University.

1. Scholes GD, et al. (2017) Using coherence to enhance function in chemical and biophysical systems. Nature 543:647–656.
2. Scholes GD, Fleming GR, Olaya-Castro A, van Grondelle R (2011) Lessons from nature about solar light harvesting. Nat Chem 3:763–774.
3. Caram JR, Lewis NHC, Fidler AF, Engel GS (2012) Signatures of correlated excitonic dynamics in two-dimensional spectroscopy of the Fenna-Matthew-Olson photosynthetic complex. J Chem Phys 136:104505.
4. Kruger TPJ, van Grondelle R (2016) Design principles of natural light-harvesting as revealed by single molecule spectroscopy. Phys B Condens Matter 480:7–13.
5. Muh F, et al. (2007) α-Helices direct excitation energy flow in the Fenna-MatthewsOlson protein. Proc Natl Acad Sci USA 104:16862–16876.
6. Rebentrost P, Mohseni M, Kassal I, Lloyd S, Aspuru-Guzik A (2009) Environmentassisted quantum transport. New J Phys 11:033003.
7. Novoderezhkin VI, Romero E, Prior J, Grondelle R (2017) Exciton-vibrational resonance and dynamics of charge separation in the photosystem II reaction center. Phys Chem Chem Phys 19:5195–5208.
8. Fuller FD, et al. (2014) Vibronic coherence in oxygenic photosynthesis. Nat Chem 6:706–711.
9. Romero E, et al. (2014) Quantum coherence in photosynthesis for efficient solarenergy conversion. Nat Phys 10:676–682.
10. Dean JC, Tihana M, Toa ZSD, Oblinsky DG, Scholes GD (2016) Vibronic enhancement of algae light harvesting. Chem 1:858–872.
11. Christensson N, Kauffmann HF, Pullerits T, Mancal T (2012) Origin of long-lived coherences in light-harvesting complexes. J Phys Chem B 116:7449–7454.
12. Chin AW, et al. (2013) The role of non-equilibrium vibrational structures in electronic coherence and recoherence in pigment–protein complexes. Nat Phys 9:113–118.
13. Kreisbeck C, Kramer T (2012) Long-lived electronic coherence in dissipative exciton dynamics of light-harvesting complexes. J Phys Chem Lett 3:2828–2833.
14. Kolli A, O’Reilly EJ, Scholes GD, Olaya-Castro A (2012) The fundamental role of quantized vibrations in coherent light harvesting by cryptophyte algae. J Chem Phys 137:174109.
15. Perlik V, et al. (2015) Vibronic coupling explains the ultrafast carotenoid-tobacteriochlorophyll energy transfer in natural and artificial light harvesters. J Chem Phys 142:212434.
16. Novelli F, et al. (2015) Vibronic resonances facilitate excited-state coherence in lightharvesting proteins at room temperature. J Phys Chem Lett 6:4573–4580.
17. Marin A, et al. (2011) Flow of excitation energy in the cryptophyte light-harvesting antenna phycocyanin 645. Biophys J 101:1004–1013.
18. Mirkovic T, et al. (2007) Ultrafast light harvesting dynamics in the cryptophyte phycocyanin 645. Photochem Photobiol Sci 6:964–975.
19. Novoderezhkin VI, Palacios MA, Amerongen H, van Grondelle R (2004) Energytransfer dynamics in the LHCII complex of higher plants: Modified Redfield approach. J Phys Chem B 108:10363–10375.

20. Shim S, Rebentrost P, Valleau S, Aspuru-Guzik A (2012) Atomistic study of the longlived quantum coherences in the Fenna-Matthews-Olson complex. Biophys J 102:649– 660.
21. Viani L, et al. (2014) Molecular basis of the exciton-phonon interactions in the PE545 light-harvesting complex. Phys Chem Chem Phys 16:16302.
22. Olbrich C, Strumpfer J, Schulten K, Kleinekathofer U (2011) Theory and simulation of the environmental effects on FMO electronic transitions. J Phys Chem Lett 2:1771– 1776.
23. Rosnik AM, Curutchet C (2015) Theoretical characterization of the spectral density of the water-soluble chlorophyll-binding protein from combined quantum mechanics/molecular mechanics molecular dynamics simulations. J Chem Theor Comput 11:5826–5837.
24. Pengfei H, Coker DF (2011) Theoretical study of coherent excitation energy transfer in cryptophyte phycocyanin 645 at physiological temperature. J Phys Chem Lett 2:825– 833.
25. Fleming GR, Cho M (1996) Chromophore-solvent dynamics. Annu Rev Phys Chem 47:109–134.
26. Aghtar M, Kleinekathofer U, Curutchet C, Mennucci B (2017) Impact of electronic fluctuations and their description on the exciton dynamics in the light-harvesting complex PE545. J Phys Chem B 121:1330–1339.
27. Lee MK, Coker DF (2016) Modeling electronic-nuclear interactions for excitation energy transfer processes in light-harvesting complexes. J Phys Chem Lett 7:3171– 3178.
28. Sinnokrot MO, Sherrill CD (2001) Density functional theory predictions of anharmonicity and spectroscopic constants for diatomic molecules. J Chem Phys 115:2439– 2448.
29. Falzon CT, Chong DP, Wang F (2005) Prediction of spectroscopic constants for diatomic molecules in the ground and excited states using time-dependent density functional theory. J Comput Chem 27:163–173.
30. Laurent AD, Adamo C, Jacquemin D (2014) Dye chemistry with time-dependent density functional theory. Phys Chem Chem Phys 16:14334–14356.
31. Send R, Kuhn M, Furche F (2011) Assessing excited state methods by adiabatic excitation energies. J Chem Theor Comput 7:2376–2386.
32. Tanimura Y, Kubo R (1989) Time evolution of a quantum system in contact with a nearly Gaussian-Markoffian noise bath. J Phys Soc Jpn 58:101–114.
33. Tanimura Y (2012) Reduced hierarchy equations of motion approach with Drude plus Brownian spectral distribution: Probing electron transfer processes by means of twodimensional correlation spectroscopy. J Chem Phys 137:22A550.
34. Kreisbeck C, Kramer T, Aspuru-Guzik A (2014) Scalable high-performance algorithm for the simulation of exciton dynamics. Application to the light-harvesting complex II in the presence of resonant vibrational modes. J Chem Theor Comput 10:4045–4054.
35. Kreisbeck C, Aspuru-Guzik A (2016) Efficiency of energy funneling in the photosystem II supercomplex of higher plants. Chem Sci 7:4174–4183.

BIOPHYSICS AND COMPUTATIONAL BIOLOGY

Blau et al.

PNAS | vol. 115 | no. 15 | E3349

36. Strumpfer J, Schulten K (2009) Light harvesting complex II B850 excitation dynamics. J Chem Phys 131:225101.
37. Ishizaki A, Fleming GR (2009) Theoretical examination of quantum coherence in a photosynthetic system at physiological temperature. Proc Natl Acad Sci USA 106:17255–17260.
38. Sumi H (1999) Theory on rates of excitation-energy transfer between molecular aggregates through distributed transition dipoles with application to the antenna system in bacterial photosynthesis. J Phys Chem B 103:252–260.
39. Scholes GD, Fleming GR (2000) On the mechanism of light harvesting in photosynthetic purple bacteria: B800 to B850 energy transfer. J Phys Chem B 104:1854–1868.
40. Jang S, Newton MD, Silbey RJ (2004) Multichromophoric Forster resonance energy transfer. Phys Rev Lett 92:218301.
41. Mirkovic T, et al. (2017) Light absorption and energy transfer in the antenna complexes of photosynthetic organisms. Chem Rev 117:249–293.
42. Novoderezhkin V, Marin A, van Grondelle R (2011) Intra- and inter-monomeric transfers in the light harvesting LHCII complex: The Redfield-Forster picture. Phys Chem Chem Phys 13:17093.
43. Bennett DIG, Amarnath K, Fleming GR (2013) A structure-based model of energy transfer reveals the principles of light harvesting in photosystem II supercomplexes. J Am Chem Soc 135:9164–9173.
44. Raszewski G, Renger T (2008) Light harvesting in photosystem II core complexes is limited by the transfer to the trap: Can the core complex turn into a photoprotective mode? J Am Chem Soc 130:4431–4446.

45. Amarnath K, Bennett DIG, Schneider AR, Fleming GR (2016) Multiscale model of light harvesting by photosystem II in plants. Proc Natl Acad Sci USA 113:1156–1161.
46. Forster T (1965) Delocalized excitation and excitation transfer. Modern Quantum Chemistry. Istanbul Lectures. Part III: Action of Light and Organic Crystals, ed Sinanoglu O (Academic, New York), pp 93–137.
47. Fleming GR, Cho M (1996) Chromophore-solvent dynamics. Ann Rev Phys Chem 47:109–134.
48. Jordanides XJ, Lang MJ, Song X, Fleming GR (1999) Solvation dynamics in protein environments studied by photon echo spectroscopy. J Phys Chem B 103:7995–8005.
49. Ishizaki A, Calhoun TR, Schlau-Cohen GS, Fleming GR (2010) Quantum coherence and its interplay with protein environments in photosynthetic electronic energy transfer. Phys Chem Chem Phys 12:7319–7337.
50. Fujihashi Y, Fleming GR, Ishizaki A (2015) Impact of environmentally induced fluctuations on quantum mechanically mixed electronic and vibrational pigment states in photosynthetic energy transfer and 2D electronic spectra. J Chem Phys 142:212403.
51. So MC, Wiederrecht GP, Mondloch JE, Hupp JT, Farha OK (2015) Metal-organic framework materials for light-harvesting and energy transfer. Chem Commun 51:3501– 3510.
52. Hemmig EA, et al. (2016) Programming light-harvesting efficiency using DNA origami. Nano Lett 16:2369–2374.
53. Kreisbeck C, Kramer T, Rodriguez M, Hein B (2011) High-performance solution of hierarchical equations of motion for studying energy transfer in light-harvesting complexes. J Chem Theor Comput 7:2166–2174.

E3350 | www.pnas.org/cgi/doi/10.1073/pnas.1800370115

Blau et al.