Trivial Excitation Energy Transfer to Carotenoids Is an Unlikely Mechanism for Non-photochemical Quenching in LHCII

Higher plants defend themselves from bursts of intense light via the mechanism of Non-Photochemical Quenching (NPQ). It involves the Photosystem II (PSII) antenna protein (LHCII) adopting a conformation that favors excitation quenching. In recent years several structural models have suggested that quenching proceeds via energy transfer to the optically forbidden and short-lived S1 states of a carotenoid. It was proposed that this pathway was controlled by subtle changes in the relative orientation of a small number of pigments. However, quantum chemical calculations of S1 properties are not trivial and therefore its energy, oscillator strength and lifetime are treated as rather loose parameters. Moreover, the models were based either on a single LHCII crystal structure or Molecular Dynamics (MD) trajectories about a single minimum. Here we try and address these limitations by parameterizing the vibronic structure and relaxation dynamics of lutein in terms of observable quantities, namely its linear absorption (LA), transient absorption (TA) and two-photon excitation (TPE) spectra. We also analyze a number of minima taken from an exhaustive meta-dynamical search of the LHCII free energy surface. We show that trivial, Coulomb-mediated energy transfer to S1 is an unlikely quenching mechanism, with pigment movements insufficiently pronounced to switch the system between quenched and unquenched states. Modulation of S1 energy level as a quenching switch is similarly unlikely. Moreover, the quenching predicted by previous models is possibly an artifact of quantum chemical over-estimation of S1 oscillator strength and the real mechanism likely involves short-range interaction and/or non-trivial inter-molecular states.


INTRODUCTION
Non-photochemical quenching (NPQ) in higher plants is a regulatory response to a sudden increase in light intensity (Horton et al., 2000;Niyogi, 2000;Müller et al., 2001;Ruban et al., 2012). It is a (mostly Malnoë et al., 2018) reversible down-regulation of the quantum efficiency of the Photosystem II (PSII) lightharvesting antenna (LHCII) with the purpose of defending the saturated reaction centers from over-excitation and photoinhibition (Powles, 1984;Aro et al., 1993).
Essentially, it is due to the creation of exciton-quenching species within LHCII which trap and dissipate chlorophyll excitation before it can accumulate in PSII and damage the reaction centers. While the fine molecular details of the mechanism are still unclear, a general consensus has emerged over the basic scheme. The primary trigger of NPQ is an acidification of the thylakoid lumen ( pH) due to a high rate of electron transport (Strand and Kramer, 2014), in large part arising from cyclic electron flow about PSI (Sato et al., 2014). The pH activates three components of NPQ: the PSII antenna sub-unit PsbS (Li et al., 2004), the enzyme violaxanthin deepoxidase (VDE) (Jahns et al., 2009), and the LHCII antenna proteins themselves (Walters et al., 1994;Liu et al., 2008;Belgio et al., 2013). VDE converts the violaxanthin pool to zeaxanthin which may lead to violaxanthin-zeaxanthin exchange in the loose, peripheral xanthophyll-binding site of LHCII (Xu et al., 2015). It has been shown that the presence of zeaxanthin affects the kinetics and amplitude of NPQ but is not a strict requirement for it (Ruban and Horton, 1999;Nicol and Croce, 2021). Quenching can similarly be achieved in the absence of PsbS if pH is driven to non-physiological levels . Either way (for further information, the reader is directed to a comprehensive review of this complex and ongoing topic Ruban, 2016) the combined effect is to induce an in-membrane aggregation or clustering of LHCII (Horton et al., 1991) and some subtle internal conformational changes (Ilioaia et al., 2011). These somehow modulate the pigment-pigment and pigment-protein couplings to create a quenching species, although the nature of the quencher and molecular dynamics of the conformational "switch" are still unclear.
Recently, it has become broadly (though by no means universally) accepted that the quencher is or involves one of the LHCII carotenoid (Cart) pigments (Ma et al., 2003). These are attractive candidates as they are, in a sense, intrinsically quenched, possessing a very short (≈10 ps) excitation lifetime relative to chlorophyll (Chl -≈4-6 ns). Various mechanisms have been suggested, such as excitation energy transfer (EET) to the Cart which quenches simply by virtue of its short lifetime (Ruban et al., 2007); mixing of the Chl and Cart lifetimes brought on by excitonic resonance (Bode et al., 2009;Holleboom and Walla, 2014); and formation of fast-relaxing Chl-Cart CT states (Holt et al., 2005;Ahn et al., 2008). The lutein (Lut) bound to the L1 binding site (Wei et al., 2016) of the LHCII trimer is often cited (Ruban et al., 2007) as being the particular carotenoid involved in quenching, but zeaxanthin at an equivalent site in one of the minor PSII antennae has also been proposed (Ahn et al., 2008). We also note that Holzwarth and co-workers present a Cart-independent quencher model that involves Chl-Chl CT states (Müller et al., 2010;Ostroumov et al., 2020). The differences between these models often come down to specific interpretations of highly-congested timeresolved spectral measurements on these complexes. Moreover, any involvement of the Carts is obscured by the fact that their lowest singlet excitation, S 1 , is optically forbidden and decays very quickly (Polívka and Sundström, 2004).
The X-ray structure of LHCII (Liu et al., 2004) can provide some insight into the quencher, particularly since it was found to correspond to a highly dissipative configuration, meaning it could serve as a model structure for the quenched state (Pascal et al., 2005). Several detailed models of this structure very accurately predicted the steady state and time-resolved spectra of LHCII (Novoderezhkin et al., 2004(Novoderezhkin et al., , 2011Müh et al., 2010) but they did not capture the dissipative character (in fairness that was never their goal). One possible reason for this was their neglect of the Carts, due to the fact that they contribute nothing to the spectrum in the red region and that there are no truly reliable methods for calculating the excitation energy and one-electron transition density of the S 1 state. The latter is due to the strong electron correlations giving it a complex multi-electron character (Tavan and Schulten, 1987;Andreussi et al., 1993). Beginning in 2013 Duffy and co-workers used a semi-empirical quantum chemistry method to estimate the S 1 transition density and its potential effect on the excitation lifetime of the LHCII crystal structure Chmeliov et al., 2015;Fox et al., 2017Fox et al., , 2018. These models suggested that quenching was due to EET from the Chl Q y band to the S 1 state of the centrally bound Luts (L1 and L2), followed by fast decay of S 1 . This EET was mediated by weak resonance couplings between Q y and S 1 (due to the latter's lack of oscillator strength) and was therefore assumed to be incoherent (Förster transfer) and slow (20-50 ps) relative to excitation equilibration across the Chls (≈1-2 ps). This is essentially the mechanism proposed by Ruban et al. based on global target analysis of transient absorption (TA) measurements on LHCII aggregate, although they propose L1 as the sole quencher (Ruban et al., 2007). Of course these models are all based on a single, timeaveraged structural configuration, and a highly artificial one at that. It therefore tells us nothing about how such quenching is switched on and off and can only very tentatively be applied to the actual in vivo quenching mechanism. More recently, the model was extended to a molecular dynamics simulation of the LHCII trimer within a lipid bilayer (Balevičius et al., 2017). Although a stable, unquenched conformation was not identified, it predicted that the Q y − S 1 coupling was highly sensitive to very small changes in inter-pigment orientations, suggesting that the lifetime could be modulated by very subtle conformational changes. Unfortunately, this appears to have been incorrect for two reasons. Firstly, the coupling sensitivity appears to have been an artifact of the semi-empirical Hamiltonian used to calculate S 1 . Khokhlov and Belov showed that this sensitivity disappears when more rigorous methods are used (Khokhlov and Belov, 2019). Moreover, by simulating the near-identical CP29, Lapillo et al. showed that even with the semi-empirical method, large lifetime fluctuations are significantly dampened if one accounts for the excitonic structure of the Chl manifold in the complex (the previous model assumed a Chl-Lut dimer embedded in some coarse-grained, iso-energetic Chl pool) (Lapillo et al., 2020). In addition to these errors, the model has a series of weaknesses that here we attempt to address: • The S 1 excitation energy is neither easy to measure directly or calculate. Transient absorption in near-IR gives a phononless excitation energy of 14, 050 ± 300cm −1 for Lut in recombinant LHCII (Polívka et al., 2002), while two-photon excitation (TPE) in native LHCII gave < 15, 300cm −1 (Walla et al., 2000). The latter value is likely the first vibronic peak which is ≈ 1100cm −1 higher than the phononless peak, meaning the two values agree reasonably well. Nevertheless, it is often treated as a free parameter and large changes to its value have been proposed as a part of the NPQ switch (Holleboom and Walla, 2014;Lapillo et al., 2020). • The vibronic structure and relaxation dynamics of S 1 were not properly considered. It was treated as a single optical transition with a line-broadening function chosen to provide a convincing visual fit to the TPE spectrum which implied very large reorganization energies. It was assumed that reorganization on S 1 was instantaneous and that internal conversion (IC) to the ground state (S 0 ) occurred with a single rate constant of ≈ 10-20ps (Polívka et al., 2002). The end result is a picture of S 1 as an deep, irreversible trap. In reality, S 1 is composed of several vibronic transitions that could couple differently to Q y , relax on finite timescales and undergo IC at different rates. • Limited sampling of the LHCII potential energy surface (PES) means we might not be probing biologically relevant conformations. Unsteered MD simulations start from a quenched minimum close to the crystal structure. Single molecule spectroscopy has shown that LHCII trimers will spontaneously switch between quenched and unquenched states but the typical dwell time in each is of the order 1-10 s (Krüger et al., 2010), meaning prohibitively long unsteered simulations may be needed to capture this switching.
Here we attempt to correct these errors in several ways. We obtain a detailed picture of the S 1 energy gap, vibronic structure and relaxation kinetics by fitting a detailed model to the TA kinetics of Lut in pyridine. These parameters (along with a secular Redfield model of the Chl manifold Malý et al., 2016) are then used to model energy relaxation in LHCII. The LHCII model structures that we use come from an exhaustive steered search of the LHCII PES which was previously published (Daskalakis et al., 2020). The motivation is to determine whether NPQ can realistically be switched on and off simply by altering the relative distance/orientation of Lut and its neighboring chlorophylls.

Steady-State Spectra of the Chlorophyll Excitonic Manifold
The Chl-Chl relaxation dynamics are modeled according to the method in Malý et al. (2016) and briefly recapped in the Methods. For a given LHCII monomer trajectory we take a set of uncorrelated snapshots and for each calculate the population relaxation. The snapshots sample disorder in the inter-pigment excitonic couplings and the different minima in the original steered MD may reveal differences in the average couplings. We do not calculate the Chl excitation (site) energies in situ but simply take the average values reported in Müh et al. (2010). The reason for this is partly to spare computational expense and partly because these fluctuations have almost no effect on quenching (Balevicius and Duffy, 2020

Relaxation Kinetics of Lutein in Pyridine
For Lut we adopt the Vibrational Energy Relaxation Approach (VERA) (Balevičius et al., 2018;Balevicius et al., 2019) to reproduce several independent spectral measurements. The details are discussed in Methods (and section C of the Supporting Information) but essentially the four singlet electronic states (|S 0 , |S 1 , |S 2 , |S n ) are replaced by sets of vibronic states, (|i a 1 a 2 = |i |a 1 i |a 2 i ), where i is the electronic index and a 1 and a 2 are the vibrational quantum numbers associated with the high-frequency, optically-coupled C − C and C = C modes, respectively (Balevičius et al., 2018). For example, the state |0 00 corresponds to the absolute ground state and |1 10 to the first excited vibrational state of the C-C stretching mode on the first excited electronic state |S 1 . The LA is given by the sum of all vibronic transitions belonging to |S 0 → |S 2 (weighted by the Franck-Condon overlaps and the initial populations on |0 a 1 ,a 2 ) and is shown in Figure 1C alongside the experimental profile for Lut in pyridine. The fit is very good up to the blue edge of the first vibronic peak after which there is a deviation due to contributions from different geometrical conformers that are not accounted for in our model (Lukeš et al., 2011). The rise above 27, 500cm −1 is a solvent artifact.
The static properties of S 1 and all dynamical properties were obtained by fitting the VERA model to the TA of Lut in pyridine. Figures 2A,B show the calculated and experimental difference spectra at intermediate (1 − 20ps) and long (10 − 42ps) delay times, respectively. The sub-picosecond kinetics are not shown as they are less relevant to the final quenching model. The S 1 Excited State Absorption (ESA, positive feature around 18, 000cm −1 ) is well fit but there is some discrepancy for the Ground State Bleach (GSB, negative feature) at earlier times. While the fit can be improved by adjusting the S 2 parameters, this disrupts the original LA fit. This could be linked to GSBdistorting local heating effects which have previously been reported  or simply an artifact. Either way it is the S 1 parameters and kinetics that we are primarily interested in. All fitting parameters are reported in section A of the Supporting Information but there are a few key quantities: • The phononless S 1 energy, ε S 1 = 14, 050 cm −1 : We assumed the previously reported value during the fit to reduce the number of free parameters. Varying ε S 1 naturally ruins the fit but it can be recovered by adjusting other parameters (mainly ε S n and the dimensionless displacements between S 1 and S 0 . As an independent check we calculated the S 0 → S 1 lineshape and compared it to the TPE of Lut in octanol (Walla et al., 2002). This is shown in Figure 1D and apart from a 150 cm −1 blue shift to account for the different solvent there is a reasonable visual agreement. However, we must note that the fit (nor the data, really) does not match the mirror image of the S 1 FL line-shapes observed for Carts such as neurosporene and spheroidene Fujii et al. (1998). These have a much less defined vibronic structure and deconvolution suggests that the largest peak is the 0 − 2 line (|0 00 → |1 01 in our model) rather than 0 − 1. We found it impossible to reproduce such a lineshape while retaining any kind of fit to the TA and TPE data. This may be a limit of the displaced oscillator model but it was later suggested that S 1 FL measurements may be distorted by the presence of cis-isomers (Christensen et al., 2007). • The S 1 lifetime, τ S 1 →S 0 ≈ 14ps: Figure 2C shows the simulated evolution of total population on S 2 , S 1 and S 0 . S 2 → S 1 internal conversion (IC) occurs on the ≈ 100fs timescale while S 1 undergoes near mono-exponential decay in within the 10 − 20 ps range usually quoted for xanthophylls. • Vibrational relaxation on S 1 , τ vib−S 1 ≈ 1ps: Figure 2D shows the simulated population evolution of the S 1 vibronic levels |1 00 , |1 10 and |1 01 . While it is difficult to assign a single lifetime to a multi-component process it is clear that vibrational relaxation on S 1 is an order of magnitude faster than S 1 IC but far from instantaneous. • Vibrational relaxation on S 0 , τ vib−S 0 < 14ps: There is a very small transient population on |0 01 and |0 10 which reaches a peak at ≈ 9 ps (not shown) and makes a very small contribution to the blue shoulder on the S 1 ESA. However, unlike Carts such as canthaxanthin and rhodoxanthin , the S 1 IC is too slow to generate a vibrational population inversion on S 0 and hence there is no S * -type signal ).

Excitonic States and Intermolecular Couplings
As a baseline we first simulated relaxation in the LHCII crystal structure (Liu et al., 2004) following protonation and minimization. The Chl-Chl couplings are essentially as previously reported (Novoderezhkin et al., 2004;Müh et al., 2010) and diagonalization leads to a set of exciton states that have already been described elsewhere (Novoderezhkin et al., 2004). is localized on the Chl a610-a611-a612 domain closely associated with Lut1, which we label |TE − . There is also an analogous 'anti-bonding' state, |TE + , at around 15, 120cm −1 . This is shown diagrammatically in Figure 3. We then consider the purely electronic Chl-Lut1 couplings, J 0 n,Lut using Lut transition charges from Khokhlov and Belov (2019). In the site basis we have the same picture as previously reported (Chmeliov et al., 2015), weak (10 − 20cm −1 ) couplings to the terminal emitter Chls and negligible couplings otherwise. Excitonic mixing among the Chls naturally mixes these couplings, the strongest (|J 0 i,Lut | ≈ 7cm −1 ) being between Lut1 and |TE − and |TE + . These small couplings justify (Balevicius and Duffy, 2020) our mixed kinetic model in which the Chl excitonic and Lut1 vibronic subsystems can exchange energy incoherently. When we model the relaxation kinetics the presence of Lut1 results in a decreased excitation lifetime of τ ex ≈ 500 ps, compared to the unquenched value of 4ns (Pascal et al., 2005). The pathway is two-fold, involving fairly-reversible transfer from |TE + to the nearresonant |S 10 1 = |1 10 level and steep down-hill transfer from |TE − to |S 00 1 = |1 00 (see Figure 4A). Recent femtosecond stimulated Raman spectroscopy (Artes Vivancos et al., 2020) has also suggested that vibrationally excited states on Lut 1 can participate in light-harvesting by internal conversion from S2 and then energy transfer to the chlorophylls, which agrees well with our observation that the |1 10 state is resonant with |TE+ . The transfer is typically slow, however. For example, the rate constant for transfer from |TE − to |S 00 1 is k −1 S 00 1 ,TE − ≈ 300 ps. There are, however, several pathways that contribute, involving other exciton states and the |S 10 1 vibronic level, resulting is a net timescale of ≈ 100 ps. This is too slow for any transient accumulation of population on |S 00 1 (see Figure 4B). While these results are essentially identical to those previously reported (Chmeliov et al., 2015), it is important to realize that the absolute value of τ ex may not be quantitatively accurate, due to the fact that the Chl-Lut couplings are derived from unscaled S 1 transition charges from quantum chemical calculations (Khokhlov and Belov, 2019), but the relative changes in lifetime between minima are meaningful.

Exploring the LHCII Potential Energy Surface
We calculated the average, relative mean excitation times for several minima identified by a previous steered MD study (Daskalakis et al., 2020). It was reported that different monomers within the same LHCII trimer could access different conformational states and so we consider the monomer in our calculations. The minima are broadly classified into 'low pH' and 'neutral pH' depending on the protonation state of several lumen-exposed residues. In all minima there were fluctuations in the snapshot lifetimes, 300 < τ ex < 1000ps, but the average value varies little within the range 500 < τ ex < 600ps (see Figures 4C,D).
It is premature to say that 'all of these minima are quenched' but we can state that there is no evidence of a simple, purelygeometric switch between states with significantly different lifetimes. We find (as previously noted Fox et al., 2017) that τ ex is correlated with Lut1-Chl a612 coupling but the coupling is not sufficiently sensitive to the small movements of the pigments to produce any kind of functional transition.

S 1 Excitation Energy and Asymmetry Between Lut1 and Lut2
If we trust the TA-derived value of ε S 1 = 14, 050cm −1 then the relative arrangement of excitonic and vibronic levels is notable (see Figure 4A). |TE + and |S 10 1 are near-resonant but since |TE + acquires little exciton population at room temperature (see Figure 4B) this is not a very effective pathway for quenching. The terminal emitter state, |TE − , lies almost precisely in the middle of |S 00 1 and |S 10 1 meaning any reasonable shift in the relative energy actually increases the quenching. This is shown in Figure 5A where we alter ε S 1 to bring either |S 00 1 (ε S 1 = 14, 750cm −1 ) or |S 10 1 (ε S 1 = 13, 600cm −1 ) into resonance with |TE − . In both cases τ ex drops by around 50% to roughly 300 ps. Within the smaller range we find that changes in the energy of ε S 1 = 14, 050 ± 300cm −1 , i.e., within the error bar of the reported value, the largest decrease is by about 25%. For ε S 1 > 18, 000cm −1 the quenching disappears completely ( τ ex → Ŵ Chl = 4 ns), as has been previously reported (Lapillo et al., 2020). This is simply because there is no energetic overlap between the two sub-systems and energy transfer between them is impossible by construction.
Although we initially excluded Lut2 we then put it back into the model, assuming the same transition charges and, initially, the same excitation energy as Lut1. The binding pocket of Lut2 (superficially) mirrors that of Lut1 with weak couplings to Chls a603 and a604 which participate in several excitonic states between |TE − and |TE + . This leads to a 40% decrease in lifetime relative to the Lut1-only model. Figure 5B shows the excitation population on the Lut ground state at time t = τ ex and we see that when ε Lut2 S 1 = ε Lut1 S 1 = 14, 050cm −1 Lut2 is almost as effective a quencher as Lut1. This is contrary to the observed features of NPQ and the known properties of Lut2. The initial TA measurements that lead to the proposal of a Lut-mediated NPQ identified Lut1 as the sole quencher (Ruban et al., 2007). This is likely because Lut2 has a significantly distorted electronic structure relative to that in solution. The S 2 excitation energy of Lut 2, ε Lut2 S 2 , is significantly lower than ε Lut1 S 2 (Son et al., 2020) and if this shift is caused by twisting of the backbone then it is likely accompanied by a concomitant upward-shift in ε Lut1 S 1 Artes Vivancos et al., 2020). In fact, recent ultra-broadband 2D measurements on LHCII identified a dark state (termed S X ), lying above the Chl a Q y band, which belongs exclusively to Lut2 (Son et al., 2019). This is most likely a strongly-distorted S 1 . Figure 5B shows that quenching by Lut2 can be completely abolished if we introduce some energetic asymmetry between Lut1 and Lut2. The shifts are not actually that large with ε Lut1 S 1 = 13, 600cm −1 being almost within the error bars of the measured value of 14, 050 ± 300cm −1 and ε Lut2 S 1 = 15, 000 − 15, 500cm −1 being roughly in the region of the proposed S X state. It is important to note that this is not a rigorous analysis, which would require independent, in situ parameterization of Lut1 and Lut2. However, it points to an energetic sensitivity in the quenching pathway(s) that was absent in previous models.

DISCUSSION
The essence of the quenching mechanism investigated here (and previously proposed Balevičius et al., 2017) is that trivial geometric modulation of Chl-Lut coupling is sufficient to drive the system between quenched and unquenched states. This appears to be incorrect, as the complex simply does not seem to possess the conformational flexibility to induce significant changes in the coupling. We are not saying that the different minima do not represent different functional states or that carotenoids are not involved in quenching, merely that our model does not capture its key features. There are several possible quenching scenarios that can be discussed.

NPQ May Involve Modulation of the Properties of S 1
The point of this study was to try and cast this problem in terms of experimental parameters, specifically the S 1 energy, vibronic structure and relaxation dynamics. While the TA fits seem reliable, the data is for Lut in pyridine and obviously there is a question of whether this can be applied to Lut in protein.
Actually, the default value of ε S 1 = 14, 050cm −1 was taken from NIR measurements on Lut in LHCII, although the error bars are quite big (±300cm −1 ) and the study did not compare quenched and unquenched configurations (Polívka et al., 2002). An earlier NPQ model proposed that quenching was induced by bringing the Chl Q y band and S 1 into resonance (Holleboom and Walla, 2014), which would require either ε 00 S 1 or ε 10 S 1 ≈ 15, 000cm −1 . Balevičius Jr. and Duffy recently provided a very general physical argument as to why fine tuning of this energy gap cannot modulate quenching (Balevicius and Duffy, 2020) and showed that significant quenching is possible for large energy gaps, even if the quenching state lies above the donor state. S 1 has to be quite far above the Q y band to abolish quenching, as was recently proposed by Lapillo et al. (2020). They reported a sharp dependence of the EET rate (and overall quenching) on the energy gap when ε S 1 ≈ 18, 000cm −1 . This is simply because at this point the Q y band coincides with the steep red edge of the S 1 lineshape. However, as they point out, ε S 1 is not a free parameter and a protein-induced blue-shift of 14, 050 → 18, 000cm −1 (712 → 555 nm) would be quite large. Saccon et al. recently performed TA measurements on quenched LHCII immobilized in polyacrylamide gel (a model for NPQ)  and found that linear excitation of Lut (i.e., via S 2 ) produces the usual S 1 → S n ESA at ε S n − ε S 1 ≈ 18, 500cm −1 (540 nm). This is reasonably close to the value in pyridine, ε S n − ε S 1 ≈ 17, 900cm −1 (558 nm). Of course this is an indirect measurement and a massive shift in ε S 1 could be hidden by a correlated shift in ε S n . However, this would have a significant affect on the ESA formation and decay kinetics which is not observed.

NPQ May Involve Non-coulomb Interactions and/or Non-adiabatic Inter-molecular States
When estimating Chl-Chl couplings the Q y transition density/charges are typically re-scaled to reproduce an experimentally determined oscillator strength (Knox and Spring, 2003;Müh et al., 2010). This is clearly not an option for S 1 (for which a non-zero oscillator strength has never been measured) and therefore the absolute values of S 1 transition charges and the Chl-Lut couplings are difficult to estimate. For the planar geometry of Lut the published transition atomic charges (Khokhlov and Belov, 2019) yield a dipole moment of |µ S 1 | ≈ 0.1 − 0.2D which, although very small, is nonzero (|µ Q y | ≈ 3 − 5D for comparison Knox and Spring, 2003). It is possible that the amplitude of the S 1 transition density is overestimated and therefore so are the Coulomb couplings, Q y → S 1 transfer rates, and overall level of quenching. In fact, given that there appears to be insufficient conformational flexibility in LHCII to switch this Coulomb-mediated quenching off, it may be an artifact. The reason that it was initially considered promising was that it qualitatively matched the NPQ scheme proposed by Ruban et al. in 2007, based on TA measurements of quenched LHCII aggregates (Ruban et al., 2007). The role of S 1 was implied by global target analysis of the kinetics rather than any visually detectable S 1 signal and so must be treated with caution. Direct observation of S 1 -mediated quenching was later reported for the cyanobacterial High light inducible proteins (Hlips) which are ancestors of LHCII (Staleva et al., 2015). Hlips are perpetually quenched by ≈ 2 ps (hence observable) EET from a small pool of Chl a to β-carotene in one of the central binding pockets that are analogous to L1 and L2 in LHCII. More recently, sub-picosecond EET to Lut1 was directly observed in LHCII via ultra-fast 2D spectroscopy (Son et al., 2020). In both the EET is much faster than predicted by this or any of the previous models and it is difficult to see how such fast transfer could be Coulomb-mediated and be in any way switchable or involve an optically-forbidden transition. Cignoni et al. provide a possible answer via a detailed QM/MM study of CP29 in which short-range interactions (exchange, overlap, etc.) were found to make large contributions to the Chl-Cart couplings (Cignoni et al., 2021). These are naturally far more sensitive to minor conformational changes than the long-range Coulomb interactions. The picture gets even more complicated when one considers quenched LHCII in gel. It was recently shown that excitation of Q y results in the immediate appearance of a large-amplitude positive peak at 19, 417cm −1 (515 nm) which we'll label A 515 . This is not merely a shifted S 1 as direct excitation of Lut gave the usual S 1 → S n ESA at 18, 500cm −1 (540 nm), although A 515 is detectable at later times and may simply be initially hidden by S 1 . This suggests that S 2 → S 1 and S 2 → A 515 are competing pathways. A 515 is in the region of the S * signal which some people suggest is either a distorted S 1 or a dipole-forbidden singlet electronic state lying below S 1 (Mascoli et al., 2019). That argument aside, since the Chl ESA is typically flat and featureless, it seems reasonable to assume that A 515 is associated with the Cars, although it is difficult to assign it solely to Lut1. A 515 is independent of whether it is Chl a or Chl b that is excited and the GSB bands in the S 2 region (< 500 nm) looks very different to the classic S 2 GSB. This all suggests some type of delocalized quenching pathway that involves several Carts and possibly even some non-adiabatic intermolecular states not accessible simply by exciting S 2 . This is exactly the picture emerging from the elaborate QM/MM models of CP29 being reported by Mennucci et al. Lapillo et al., 2020).

Quenching Requires Hydrophobic Mismatch and Aggregation
It is possible that the conformational switch cannot be revealed by simulating a single LHCII monomer/trimer in a lipid bilayer. In vitro quenching is induced by low detergent concentration which in solution leads to aggregation. LHCII aggregates are the original model system for studying NPQ (Horton et al., 1991) and there is compelling evidence that some form of aggregation or clustering in the membrane is part of the in vivo mechanism . Key to this is PsbS, with overexpression observed to enhance LHCII clustering and its absence frustrating it (Goral et al., 2012). Recent simulations show that PsbS's lumen-exposed side is covered in titratable residues with protonation causing an unfolding of a specific region implicated in protein-protein interactions . Other studies have shown that it is capable of interacting with the minor PSII antenna complexes (Daskalakis, 2018), possibly helping LHCII to partially detach from the reaction center complex and form the quenching clusters. It has also been suggested that active PsbS has an affinity for certain lipids, altering local membrane composition and causing the hydrophobic mismatch that drives aggregation/clustering (Daskalakis et al., 2019;Ruban and Wilson, 2020).
The role of external interactions on LHCII conformation was (at least partially) considered in the metadynamical simulations used in this work (Daskalakis et al., 2020). They were performed by first calculating the free energy surface (FES) of LHCII while considering a wide range of external stimuli, such as pH and interactions with PsbS. The next step involved steering MD simulations of an LHCII trimer around this FES in order to reach the minima. Equilibrium MD trajectories were then performed on an isolated trimer once each minimum was reached. Due to the absence of external stimuli in these equilibrium trajectories, it is possible that the configuration space even at these minima still favors the quenched conformation over the unquenched.

TA Measurements of Lutein in Pyridine
All transient absorption data were measured with a spectrometer described in detail in Saccon et al. (2020). Lutein (Sigma Aldrich) was dissolved in spectroscopic grade pyridine to yield an optical density of ≈ 0.2 mm −1 at the absorption maximum. The sample was placed in a 2 mm path-length quartz cuvette equipped with a micro-stirrer to avoid sample degradation during measurement. The mutual polarization of the excitation and probe beams was set to the magic angle (54.7 • ) and excitation intensity was kept below 1014 photons pulse −1 cm −2 .

The Chlorophyll Exciton Manifold
Modeling of energy relaxation within the chlorophylls is carried out according to previous work (Novoderezhkin et al., 2004;Malý et al., 2016) and is described in detail in section B of the Supporting Information. Briefly, for a single uncorrelated MD snapshot (at time t k ) the relevant system of Chl excited (Q y ) states is determined by the usual spin-boson Hamiltonian, where ε = ε r ε 0 = 2ε 0 . Both {E n } and q α are taken from Müh et al. (Madjet et al., 2006;Renger et al., 2007). Equation (1) is then diagonalised to give the exciton basis, where c i n are the participation coefficients of each pigment state, |n for a given exicton state, |i . Site energies, oscillator strengths and couplings to Cart vibronic levels are also mixed. The exciton states are initially populated according to their oscillator strengths and relaxation is modeled using the approach outlined in Novoderezhkin et al. (2004). The population relaxation rates are given by, where C ′′ n (ω) is the spectral density of bath-induced site energy fluctuations and ω ij is the gap between the zero-phonon lines of excitons i and j. The ansatz spectral density (Novoderezhkin et al., 2004) is assumed throughout. For a single uncorrelated snapshot along a trajectory (at time t k ) the instantaneous LA and FL spectra are given by and,

The Carotenoid Vibronic Subsystem
The full VERA (Balevičius et al., 2018;Balevicius et al., 2019) Hamiltonian is, ε a 1 a 2 i |i a 1 ,a 2 i a 1 ,a 2 | (7b) c iκ x κ √ a 1 + 1 |i a 1 ,a 2 i a 1 +1,a 2 | + |i a 1 +1,a 2 i a 1 ,a 2 | (7d) c iκ x κ √ a 2 + 1 |i a 1 ,a 2 i a 1 ,a 2 +1 | + |i a 1 ,a 2 +1 i a 1 ,a 2 | (7b) is the Hamiltonian of the system (H S ) of uncoupled vibronic levels, {|i a 1 a 2 }, where ε a 1 a 2 i = ǫ i + ǫ a 1 ,a 2 = ǫ i + a 1 + 1 2 h ω 1 + a 2 + 1 2 h ω 2 (8) is the sum of the electronic, ǫ i and vibrational energies. H B is the bath Hamiltonian which is composed of a large set of harmonic oscillators representing the non-optical modes of the Cart itself plus librations, solvent modes, etc. We split this into two parts, H IVR SB and H IC SB . H IVR SB describes the bath-induced couplings between adjacent vibrational levels of the optical modes and are therefore responsible for vibrational relaxation on the electronic states. {c iκ } are the coupling constants and {x κ the bath mode displacements. Energy (mostly) relaxes into the non-optical modes of the Cart and therefore reflects Intramolecular Vibrational Redistribution (IVR). Note that there is no population transfer between the two optically-coupled modes. H IC SB couples different electronic states and is therefore responsible for Internal Conversion (IC). It is characterized by coupling constants, f iκ , and the Franck-Condon (FC) overlaps, If we assume that the frequencies of the optical modes (ω α ) are independent of the electronic state (i.e., no Duschinsky rotations) then {F α i aα ,j bα } are entirely determine by their relative dimensionless displacements, {d ij α }. The relaxation dynamics are obtained by a second-order perturbative treatment of H IVR SB and H IC SB . The resulting equations of motion are rather complicated and are listed in section C of the Supporting Information. The various IVR, k ± α , and IC, k ij a 1 a 2 ,b 1 b 2 , rate constants are defined in terms of Drude-type spectral density functions C ′′ i α (ω) and C ′′ i α ,j α (ω). We therefore have a large set of fitting parameters including electronic, {ǫ i } and vibrational, ǫ a 1 a 2 , energies, modes frequencies, ω α , mode displacements, d ij α , and the reorganization energies, λ i α , λ i α ,j α , and dephasing frequencies, γ i α , γ i α ,j α . Solving the dynamics yields a set of vibronic populations, n i a 1 a 2 (t), which are used to calculation the TA difference spectrum as a combination of ESA, GSB and stimulated emission (SE) components, which are given by, A X (ω, t) = i,a 1 ,a 2 n i a 1 a 2 (t)I a 1 a 2 ,X (ω) where I a 1 a 2 ,X (ω) are FC-weighted Gaussian/Lorentzian lineshape functions that account for line-broadening.

Energy Transfer Between the Chlorophyll and Lutein Subsystems
Having parameterized the subsystems separately we can now couple them via the calculated resonance couplings, J 0 n,Lut (t k ). We make two assumptions. Firstly, since the inter-pigment couplings between the Chls and Lut is an order of magnitude smaller than the nearest-neighbor chlorophyll couplings (there is essentially no coupling between Lut1 and Lut2), we treat the Chl-Lut system as two weakly-interacting subsystems and assume that energy transfer proceeds incoherently (Balevicius and Duffy, 2020). Secondly, since there is almost no accumulation of vibronic population on the ground state ('hot' ground state), we do not explicitly include the Chl or Lut ground states in the dynamics. S 1 can decay to higher vibrational levels on S 0 but excitation proceeds from |S 00 0 = |0 00 . Essentially, we are assuming instantaneous IVR on the ground state. The couplings in the exciton basis are given by, where {J 0 n,Lut (t k )} are the purely electronic Chl-Lut couplings. The rate of transfer from Chl exciton state |i to Lut vibronic level |S a 1 a 2 1 = |1 b 1 b 2 is given by the Fermi Golden Rule, where ij b α a α = (ǫ i + ǫ b 1 ,b 2 ) − (ǫ j + ǫ a 1 ,a 2 ) is the vibronic energy gap,χ ′ i is the normalized excitonic fluorescence lineshape and σ (ω; 10 b α 0 α , ω 10 ) is the normalized vibronic Gaussian lineshape of width ω 10 = 1070cm −1 determined by the TA fit. The backward rate is similarly defined and Boltzmann factors are added to uphill rates to enforce the detailed balance condition.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
CD, CG, and TW devised the project. TP provided TA data and assisted with data fitting. CD performed TA data fitting. VD provided steered MD data. TW contributed pilot MD data for development of model. All calculations and overall model development by CG. CG and CD wrote the draft, all authors took part in editing.
Foundation (19-28323X) for financial support. VD would like to thank European Regional Development Fund and the Republic of Cyprus through the Research and Innovation Foundation (POST-DOC/0916/0049) for supporting the metadynamics studies of LHCII trimers.