Theoretical Characterization of the Spectral Density of the Water-Soluble Chlorophyll-Binding Protein from Combined Quantum Mechanics/Molecular Mechanics Molecular Dynamics Simulations.

Over the past decade, both experimentalists and theorists have worked to develop methods to describe pigment-protein coupling in photosynthetic light-harvesting complexes in order to understand the molecular basis of quantum coherence effects observed in photosynthesis. Here we present an improved strategy based on the combination of quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations and excited-state calculations to predict the spectral density of electronic-vibrational coupling. We study the water-soluble chlorophyll-binding protein (WSCP) reconstituted with Chl a or Chl b pigments as the system of interest and compare our work with data obtained by Pieper and co-workers from differential fluorescence line-narrowing spectra (Pieper et al. J. Phys. Chem. B 2011, 115 (14), 4042-4052). Our results demonstrate that the use of QM/MM MD simulations where the nuclear positions are still propagated at the classical level leads to a striking improvement of the predicted spectral densities in the middle- and high-frequency regions, where they nearly reach quantitative accuracy. This demonstrates that the so-called "geometry mismatch" problem related to the use of low-quality structures in QM calculations, not the quantum features of pigments high-frequency motions, causes the failure of previous studies relying on similar protocols. Thus, this work paves the way toward quantitative predictions of pigment-protein coupling and the comprehension of quantum coherence effects in photosynthesis.


Introduction
Photosynthesis is the process by which sunlight is captured and converted into chemical energy, ultimately sustaining the life of essentially all organisms on earth. Despite its importance to life, the actual mechanism by which photons are captured and their energy is passed among pigment-protein complexes toward reaction centers is still heavily debated. In particular, the recent discovery that electronic energy transfer (EET) in photosynthesis involves quantum coherence, or correlation between quantum states in multiple objects, has challenged the traditionally assumed Förster incoherent hopping mechanism. [1][2][3] Understanding the molecular basis sustaining quantum coherence could lead to the design of more robust, affordable means to obtain renewable energy, e.g., a new generation of solar cells. 4 Even if the impact of quantum coherence on light harvesting efficiency is still controversial, 3 such knowledge could also be useful for designing quantum computing technologies, where a long-standing goal is to maintain coherences among qubits at room temperature. 5 Researchers previously considered that quantum coherence could not be involved in photosynthetic light harvesting due to the thermal disorder present in biological systems.
However, when coherences were observed to survive for hundreds of femtoseconds, [6][7][8] it became evident that the protein environment and its structural dynamics play a critical role in the light harvesting process. 2,9,10 In particular, quantum coherence seems to be intimately related to the structured nature of the protein environment, as opposed to the behavior of a random thermal noise that should favor decoherence. One hypothesis supposes that correlated fluctuations in the individual pigments electronic energy levels, or site energies, could give rise to long-lived coherences. 7 Nonetheless, theoretical studies combining quantum chemistry and classical molecular dynamics (MD) simulations predict such correlations to be negligible. 11,12 Another popular hypothesis emphasizes the coupling of protein vibrations to pigments excitations, which could support coherence in cases where the vibrational energy commensurates the energy gap among exciton states. 13,14 The emerging consensus suggests that EET in photosynthesis often occurs in an intermediate coupling regime. This implies that neither pigment-pigment nor pigment-protein interactions can be treated perturbatively when solving the dynamics of EET, as is done in the popular Förster or Redfield theories, because the magnitudes of these interaction strengths are similar. Much progress has been made to develop theoretical models simulating EET dynamics in the intermediate regime. 3 That said, the outcome of these models critically depends on the tools used to parameterize the photosynthetic system at hand -namely, the site energies, the electronic couplings, and the spectral density of exciton-phonon coupling that quantifies the pigment-protein interactions. Whereas recent research has made substantial progress towards accurately determining model Hamiltonian parameters (site energies and couplings) both from theory and from experiment, [15][16][17] the determination of the spectral density of individual pigments in a pigment-protein complex -the quantity needed to understand the molecular basis of quantum coherence -still presents a great challenge to the field, primarily since a large number of protein vibrations couple to many pigment electronic transitions over multiple timescales.
Consequently, it is quite difficult to construct experiments in which the correlation between pigments electronic excited states and protein vibrational states can be discerned. Site-selective spectroscopies provide an appealing means to determine spectral densities; however, these techniques often only probe the lowest-energy fluorescing pigments in the complex. 18 An intriguing alternative to experiment is the theoretical determination of spectral densities. Several groups have contributed important advances in this direction. [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34] The typical technique involves the extraction of the spectral density from the Fourier transform of the autocorrelation function of site energy fluctuations, the latter estimated using mixed quantum mechanical/molecular mechanical (QM/MM) calculations performed along a classical MD trajectory. 19,22,24 This strategy has been applied to several pigment-protein complexes by performing semi-empirical Zerner's intermediate neglect of differential overlap (ZINDO) or time-dependent density functional theory (TD-DFT) calculations in the presence of the protein atomic charges, which mimic the electrostatic influence of the protein on the pigments' excited states. Recently, we have explicitly accounted for protein polarization effects on the pigmentprotein coupling using a polarizable force field. 31 Nevertheless, the spectral densities predicted in this way systematically show strongly overestimated peaks in the high frequency region of the spectrum, a problem often ascribed to the classical nature of the MD simulations performed, which prevents the capture of the quantum mechanical features of such high-frequency modes.
Another limitation of this strategy arises from the "geometry-mismatch" problem, as low-quality geometries extracted from the classical simulation are used for QM excited state calculations.
In this publication, we present an improved strategy based on the combination of semiempirical QM/MM MD simulations 35 with polarizable QM/MMpol calculation of the excited state properties. 36,37 We apply this novel approach to the prediction of spectral densities of Chl a and Chl b in class IIb water-soluble chlorophyll-binding protein (WSCP). Our results are critically compared with predictions from classical MD simulations as well as from detailed differential fluorescence line-narrowing data. [38][39][40] They reveal that when the MD trajectory is obtained along a QM/MM potential energy surface, theoretical predictions can capture the intramolecular features characterizing the spectral densities of the chlorophylls in WSCP to a surprisingly high degree of accuracy.
The term "WSCP" actually implies two different classes of proteins, of which there are multiple types. The WSCP complex studied here corresponds to class IIb WSCP from Lepidium virginicum (Virginia pepperweed), whose crystal structure was solved by Horigome et al. 41 This structure is organized as a homotetramer, where each monomer binds one Chl a or Chl b molecule, as illustrated in Figure 1. This leads to a dimer-of-dimers organization of the chlorophylls. WSCP is not a photosynthetic protein but rather a chlorophyll carrier protein found in higher plants; its proposed function is to scavenge for and hoard chlorophyll when the plant experiences stress. 41 The tractable size, dimer-of-dimer symmetry, presence of multiple pigments, and proximity of the chlorophylls within a dimer greatly simplify complexities that are observed in photosynthetic complexes, making WSCP a simple yet effective model to investigate in detail the role of pigment-pigment and pigment-protein interactions in photosynthetic systems.
Most theoretical and spectroscopic studies, including steady-state spectra, transient absorption, difference fluorescence line-narrowing, and hole-burning experiments, have so far addressed the properties of class IIa WSCP, 38,39,[42][43][44][45] which binds up to two chlorophylls per complex. Only recently, 2D electronic spectroscopy experiments have also been presented for class IIb WSCP. 46 Some studies have shown that the crystal structure of class IIb WSCP can be used to describe the optical properties of class IIa WSCP, suggesting a similar Chl binding motif in both cases. 47,48 Following this idea, we adopt the class IIb crystal structure of WSCP to simulate the spectral density determined from differential fluorescence line narrowing by Pieper and co-workers 38 on class IIa WSCP reconstituted with either Chl a or Chl b from cauliflower (Brassica oleracea var. botrytis). We then perform QM/MM as well as classical MD simulations to explore the groundstate potential energy surface of the system and later process the trajectory with ZINDO excitedstate calculations based on the QM/MMpol approach we have recently developed, 36,37 which accounts for mutual pigment-protein polarization effects in a self-consistent way, in order to extract the corresponding spectral densities of Chl a and Chl b in the complex. The paper is organized in the following way. First, the methods and computational details -the extraction of the spectral density from the autocorrelation function and the details of the MD and the QM/MMpol excited-state calculations -are presented. Next, the results for both WSCP with Chl a and WSCP with Chl b are shown and discussed. Finally, we provide the conclusions of our study and some perspectives on future developments for quantitative predictions of pigment-protein coupling in photosynthetic complexes.

Spectral densities
The spectral density, J(ω), describes the frequency-dependent coupling strength among electronic transitions and matrix vibrations, including both intramolecular pigment and intermolecular environmental modes. For practical purposes, it is often decomposed into two contributions: where J 0 (ω) describes the coupling of the pigments' excitations to a continuum of low-frequency damping modes due to the solvent and the protein environment, and J vib (ω) accounts for the coupling with high-frequency modes, mostly coming from intramolecular pigment vibrations.
The spectral density defines the total reorganization energy The expression in Eq. 3 has a classical pre-factor introduced to negate the temperature dependence of the classical correlation function, thus ensuring a temperature-independent spectral density. 24 In the following sections, we describe the classical and PM6 QM/MM MD simulations as well were multiplied by a Gaussian envelope of variance ! = 3.6 × 10 ! fs 2 to ensure they converge to zero. 24 The high-frequency offset in the spectral densities was corrected with negligible effects in the low-frequency part by shifting the Fourier transform before its multiplication by the prefactor. 28,29 Here our theoretical estimates are compared to the experimental differential fluorescence linenarrowing (ΔFLN) data measured by Pieper and co-workers. 38,39 In these experiments, both the continuous J 0 (ω) and the discrete J vib (ω) components of the spectral density were estimated. In order to model J 0 (ω), however, the log-normal expression and the corresponding parameters introduced by Kell and co-workers to describe the Chl a and Chl b WSCP ΔFLN spectra (parameters reported in Table 1 of Ref. 40 ) are used; such an approach has been shown to describe better the phonon spectral density of different photosynthetic complexes compared to other expressions commonly used. On the other hand, the discrete J vib (ω) component is modelled as a collection of damped quantum harmonic oscillators using the vibrational frequencies and Huang-Rhys factors derived from the ΔFLN spectra (reported in Table 1 of Ref. 38). The complete experimental spectral densities were thus modelled using the following expressions: where ! is the cutoff frequency, σ ! the standard deviation, and S ! the Huang-Rhys factor for the phonon contribution i. The discrete oscillators have frequencies ! with associated damping ! and reorganization energy ! = ! ! . The associated damping ! was set to 16.5 cm -1 in order to reproduce the overall Huang-Rhys factors 0.80 and 0.74 measured for the vibrational profile when Eq. 2 is applied to the set of discrete vibrations of Chl a and Chl b WSCP, respectively 38 .

Molecular dynamics simulations
The Chl a and Chl b WSCP simulation systems were based on the X-ray crystal structure of class IIb WSCP homotretramer reported at 2.0 Å resolution by Horigome et al. (PDB code 2DRE). 41 The four chlorophylls in the asymmetric binding sites of the crystal were assigned as Chl a molecules (Chl a WSCP), so for the Chl b simulations we substituted those with four Chl b molecules. The missing residues in the PDB file were carefully added by selecting appropriate rotamers for the side chain groups in order to avoid any steric clash with the neighboring amino acids. The structure is organized as a dimer of dimers, despite the sequence identity of the four subunits. Protonation states of all residues were determined by computing their corresponding pKa's at neutral pH using the PROPKA3 server, 49,50 indicating a standard protonation state for all residues. The system was solvated in a TIP3P water box (a truncated octahedron with a buffer zone of 12 Å) using the Leap module of the Amber 12 suite of programs. 51 The protein was described using the parm99SB Amber force field. 52 In the classical MD simulations, Chl a was described using the parameters developed by Zhang and co-workers 53 based on a previous parameterization of bacteriochlorophyll performed by Ceccarelli et al., 54 with slight modifications. First, in Zhang's force field the four nitrogen atom types correspond to one ns and three nmh. In order to prevent a distortion of the in-plane position of the Mg atom occurring due to the nmh-mgc-nmh and nmh-mgc-ns angle bendings, a new nitrogen atom type was introduced (note that originally such bending angles had an equilibrium value of 90º, inappropriate when the two nitrogens involved are in opposite sides of the Mg atom, where this value should rather be 180º). Second, a few missing dihedral terms for Chl a were added, and new atom types and parameters needed to describe Chl b were assigned based on the general Amber force field (gaff). 55 Finally, the RESP charges for both Chls were estimated via the usual Amber RESP approach at the HF/6-31G(d) level using a geometry optimized at the B3LYP/6-311G(d,p) level of theory. The final system, which contained ~90,000 atoms, first underwent 50 ns of MD simulation at constant pressure (1 atm) and temperature (298 K) using standard coupling schemes using the Amber12 code. 51  For the Chl a-WSCP and Chl b-WSCP systems, 100 ps QM/MM MD simulations were run using the built-in semi-empirical QM/MM code available in Amber 12. 35 These runs were started from the same point as the 100 ps classical trajectory, i.e., after 50 ns of classical MD simulation. The same conditions as for the classical MD were used, except the fact that Chl-1 and Chl-2 were described at the semi-empirical PM6 QM level of theory, 56 whereas Chl-3 and Chl-4 molecules were described through the classical force field to make the simulations feasible. Dispersion and hydrogen bond corrections with PM6 were not implemented because there were no available parameters for Mg. Finally, a QM-MM cutoff (qmcut) of 8 Å was utilized. In simulations that use periodic boundary conditions, the QM/MM cutoff only applies to the division of interactions between "direct" and "reciprocal" parts, and 8 Å is the recommended value. 51 As for the classical trajectory, 20,000 structures were extracted every 5 fs for the excited-state QM/MMpol calculations.

Polarizable QM/MM calculations
Excited state calculations along the structures extracted from the MD simulations were performed using the polarizable QM/MM strategy developed to study excited states and energy transfer in heterogeneous environments. 36,37,57 In particular, excited state calculations were performed for each chlorophyll, and the corresponding electronic couplings were calculated between the lowest-energy π→π * Q y states of the chlorophylls at the ZINDO level of theory. 58 For structures extracted from the QM/MM MD simulation, the excited state properties of Chl-1 and Chl-2 were computed, but no electronic couplings were estimated because the other chlorophylls were treated classically. The QM/MMpol model employed here couples a polarizable classical description of the protein and solvent environment based on the induced dipole model with a QM description of the pigments excited states. This requires assignment of atomic charges and polarizabilities to all MM atoms, and mutual polarization effects among the QM and MM regions are solved during the self-consistent-field process and the excited-state environment response. A detailed description of the methodology can be found in previous publications. 36,57 Electronic couplings between the Q y states also were computed using the structures extracted from the classical MD trajectory. In particular, the couplings are estimated from the transition densities derived from the ZINDO QM/MMpol excited state calculations, using the following expression: 36,37 where !/! ! indicates the transition densities of the interacting pigments, and ! !"# the induced dipoles describing the MM polarization response of the environment to a given transition density. The !"#$ term describes the Coulomb interaction between the transition densities, thus representing an extension of the Förster dipole-dipole term that accounts for the shape of the interacting molecules. The !"# term describes the environment-mediated contribution, which typically counteracts the !"#$ term and therefore leads to an overall screening of the coupling.
We can thus define an effective screening factor s as: 37 In all QM/MMpol calculations the AL parameters of the recent Amber pol12 polarizable force field to describe the protein MM region were adopted. For the water solvent, the POL3 model charges were used. 59 To speed up the calculations, the water molecular polarizability (1.44 Å 3 ) located on its center of mass was used instead of a three-point distributed polarizability model.

Spectral densities of electron-vibrational coupling
In Figures 2 and 3 we report the spectral densities derived from QM/MM calculations pigment. In addition, this finding suggests that the 100 ps window of the MD trajectory suffices to converge the primary characteristics of the spectral densities, except for very low-frequency modes, which arise from slow protein motions characterized by longer timescales, which would require longer simulations to resolve these features.
In Figures 2 and 3, our predictions are plotted along with the spectral density measured by Pieper and co-workers from differential fluorescence line-narrowing spectra (ΔFLN) at 4.5 K for class IIa WSCP from cauliflower. [38][39][40] Compared to previous reports of simulated spectral densities, the availability of ΔFLN data reporting the coupling strength of all vibrations in the frequency ranges of interest, i.e., from the environmental low-frequency region up to high-frequency intramolecular chlorophyll modes, allows us to critically assess the capability of the combined QM/MD strategy to describe the spectral density of this pigment-protein complex. Indeed, for Chl a the theoretical spectral density spans nearly the same range observed in the experiments, where the last peak appears at ∼1650 cm -1 . However, for Chl b the simulated spectral density predicts two peaks at ∼1730 cm -1 and ∼1870 cm -1 at frequencies somewhat larger than 1681 cm -1 , the highest-frequency peak observed in experiments. As discovered by other groups, simulations dramatically overestimate the coupling to high frequency intramolecular motions at ∼1600 cm -1 .
This problem, previously found in other pigment-protein complexes, has been largely attributed to the classical nature of the MD trajectory used, which neglects quantum features of such high frequency motions. 17 As shown by Jing et al. 26 the displaced equilibrium geometry of the pigments described by the classical force field compared to the corresponding QM one consistent with the excited state calculations leads to an overestimation of the coupling strength. The classical force field, however, can also misrepresent the distribution of vibrational modes, i.e. the position of the peaks, due to the simplicity of the force field energy function. Indeed, the rich structure of peaks observed in the 1000-1600 cm -1 range from ΔFLN spectra is lost in the simulated spectral densities, or at most strongly underestimated, with the exception of the exaggerated peak at ∼1600 cm -1 and the accurate peaks near ∼1550-1600 cm -1 . The underestimation of the spectral density in the 1000-1600 cm -1 range thus seems to arise from an incorrect description of the position of the high-frequency vibrations, which incorrectly accumulate near ∼1600 cm -1 in contrast to experimental observation, rather than due to an incorrect description of the coupling strength due the geometry mismatch problem. This effect is clearly observed in the total accumulated reorganization energies shown in Figure 4, which abruptly increase near ∼1600 cm -1 .
In sharp contrast, the simulated spectral densities at frequencies below 1000 cm -1 reproduce the experimental curve surprisingly well, with the notable exception of the exaggerated peak at ∼280 cm -1 predicted in the Chl a WSCP complex, and the peak at Thus, the middle and low-frequency motions of the chlorophylls seem to be quite accurately described by the classical simulations coupled to the QM/MM semi-empirical estimation of the excited states, though the modes beyond 1000 cm -1 are poorly described.
The adoption of QM/MM MD simulations immensely improves the spectral density in the highfrequency region (750-1500 cm -1 ), as shown in Figure 3. For the Chl a WSCP complex, all major experimental peaks in this range are observed with a corresponding peak in our simulations, although the intensity of some peaks is strongly overestimated, most notably at ∼800 cm -1 , ∼840 cm -1 , ∼920 cm -1 , and ∼990 cm -1 . In contrast to the simple force field energy function, thus, the PM6 QM/MM method seems to recover the proper distribution of pigment intramolecular vibrations in this frequency range suggested by experiments. There is a high degree of matching in terms of both frequency and intensity for ∼1150 -1500 cm -1 . Nevertheless, some highfrequency peaks are still inaccurate. In particular, the feature at ∼1550 cm -1 overestimates the electron-vibrational coupling, and the strong peak at ∼1750 cm -1 is not observed by experiment.
In sharp contrast, the description of the low and mid-frequency regions (50- The quality of the overall predicted spectral densities can be judged by inspection of Figure 4, in which we provide an estimate of the total accumulated reorganization energies as a function of the vibrational frequency range used for the integration. In the theoretical estimates based on classical MD trajectories, the agreement with experiment is excellent up to ∼1000 cm -1 , illustrating the fact that the pigment-protein coupling is described better than in QM/MM simulations, probably due to the inaccuracies in the QM/MM boundary parameters previously discussed. After ∼1000 cm -1 , however, the slope of the curve is better when the QM/MM MD trajectories are adopted, whereas the classical MD curves show a flat trend until ∼1600 cm -1 , when the strongly exaggerated peak observed in the spectral densities gives a step-like increase in the reorganization energy not observed in the experimental counterpart. As we argue previously, this suggests that the classical force field is unable to describe the smooth

Excitonic parameters
The relative simplicity of WSCP compared to other pigment-protein complexes makes it a  Table 1    Lastly, another potential problem is the assumption of equal site energies in the experimental estimate derived from excitonic splittings. 39 If the best coupling estimates are taken as 38 cm -1 and 35 cm -1 , which include dielectric screening effects and are consistent with Knox analysis, an asymmetry of about ∼150 cm -1 would be needed to reproduce the experimental excitonic splittings. Such asymmetry could arise due to different packing of the phytyl chains in the hydrophobic pocket. 41 Indeed, our simulations predict energy differences among Chl-1/Chl-2 or Chl-3/Chl-4 of ∼100-150 cm -1 , as shown in Tables 2 and 3    There is much room for improvement in such theoretical studies; for instance, adopting DFT methodologies could lead to even more accurate results.
Finally, we provide an estimate of the site energies and the electronic coupling strength corresponding to the chlorophyll pigments in WSCP. The results suggest that the polarizable environment of the protein is quite heterogeneous, leading to a substantial range of dielectric screening effects as experienced by the different chlorophylls in the complex, with effective dielectric constants in the 1.5-3 range. We also critically compare our results with the experimental estimate of the electronic coupling in the strongly coupled dimer provided by Pieper and co-workers from spectral hole-burning experiments performed in class IIa WSCP by assuming symmetric site energies. 39 The results indicate that effective dipoles equal to 6.5 Debye This information is available free of charge via the Internet at http://pubs.acs.org.