Spectral variability in phycocyanin cryptophyte antenna complexes is controlled by changes in the α polypeptide chains

Quantitative models of light harvesting in photosynthetic antenna complexes depend sensitively on the challenging determination of the relative site energies of the pigments. Here we analyze the basis of the light harvesting properties of four antennae from cryptophyte algae, phycocyanines PC577, PC612, PC630 and PC645, by comparing two alternative theoretical strategies to derive the excitonic Hamiltonian. The first is based on molecular dynamics simulations and subsequent polarizable quantum/molecular mechanics (QM/MMPol) calculations, whereas the second is based on three-layer QM/MMPol/ddCOSMO calculations performed on optimized geometries of the pigments, where the water solvent is described using the ddCOSMO continuum model. We find the latter approach to be remarkably accurate, suggesting that these four phycobiliproteins share a common energetic ordering PCB82 < PCB158 < DBV51/61 for pigments located in the highly-conserved β chains, whereas bilins in the more divergent α chains originate their spectral differences. In addition, we predict a strong screening of the coupling among central DBVs in “open” form complexes PC577 and PC612 compared to “closed” form ones, which together with the increased interpigment separation explains the attenuation of coherence beatings observed for these complexes.


Introduction
Pigment-protein antenna complexes play an important role in the photosynthetic machinery by collecting sunlight and transporting it to reaction centers (RC), where the absorbed energy is used to drive charge separation events. [1] Whereas reaction center complexes are quite conserved, antenna complexes encompass a rich variety of structures and pigment compositions adapted to the needs of the photosynthetic organisms depending on their particular habitats. The remarkable quantum efficiency of the electronic energy transfer (EET) process that drives the energy to the RCs has inspired a continued effort aimed at understanding the subtle details that relate the EET process with the structure of the underlying biomolecule. [2][3][4] However, since the first high-resolution crystal structure of an antenna complex, the Fenna-Matthews-Olson (FMO) complex, appeared 40 years ago, [5] a particularly challenging aspect has been an accurate determination of the site energies -the uncoupled excitation energies of the pigments. [6] A rather successful strategy relies on the simultaneous fitting of a variety of optical spectra. [7] However, those fittings are not free from ambiguities, even in antenna complexes with a relatively low number of sites. [8] In addition, they do not allow drawing a relation between the light harvesting properties of the complex and its underlying structure, a desirable insight in order to establish structural blueprints for the design of artificial light harvesting systems. An attractive alternative to empirical models consists in the theoretical calculation of site energies from the structure of the complex, but this has proven to be a remarkable challenge, because the errors associated to present computational methods are similar to the relative site energy differences. In the last decade a variety of groups have attempted the calculation of site energies in a number of photosynthetic systems, with variable degrees of success. [2][3][4] The differences among these computational approaches are mainly related to the quality of the description used to model the pigments and the surrounding protein and solvent environment, and thus their interactions. In addition, these strategies also differ on whether the solvatochromic shift arises purely from i) pigment-protein interactions, ii) deformation of the pigment internal geometries in the protein scaffold, or iii) both effects. Calculations of the site energies in chlorophyll-containing antennae have been quite successful by limiting the mechanism of site-energy tuning to the direct electrostatic pigment-protein interactions, [4] even if contradicting reports indicate an important role played by the chlorophyll macrocycle deformation. [9] In contrast, in phycobiliproteins (PBP), which contain quite flexible linear tetrapyrrole pigments called bilins, the site energy tuning exerted by the protein has been previously ascribed mainly to the constrained conformation of the bilins in the protein pocket. [10] Internal deformation effects are however hard to estimate in a reliable way, which partially explains why neglecting this effect has led to good results in chlorophyll-based systems. Indeed, estimating this effect from the pigment geometries in a crystal structure is in general not a reliable strategy, due to the limited resolution of the crystal. [11] In order to overcome this difficulty, accurate enough geometries to carry out excited-state calculations can be derived, in principle, by using quantumchemical (QM) methods, often coupled to a classical molecular mechanics (MM) description of the environment in so-called hybrid QM/MM approaches. [12] However, proteins are characterized by complex energy landscapes, so thermal effects can have an important impact on the results. In other words, performing calculations on a single structure can result in a limited description of the properties of a thermally disordered biological system. Here, we investigate the performance of alternative theoretical strategies to determine the light harvesting properties of photosynthetic complexes by studying four PBPs from cryptophyte algae, the phycocyanins PC577, PC612, PC630 and PC645. [13] These systems are well-suited for this task, because both the PC577/PC612 and PC630/PC645 pairs share the same pigment composition, therefore their changes in spectral properties arise from the changes in the underlying protein sequences. In particular, PC577 and PC612 contain six phycocyanobilins (PCB) and two 15,16-dihydrobiliverdins (DBV), whereas PC630 and PC645 contain four PCBs, two DBVs, and two mesobiliverdins (MBV). Moreover, a single-residue insertion switches the "closed" quaternary structure from PC630/PC645 to an "open" structure for PC577/PC612, [14] as shown in Figure 1, leading to a reduction of the exciton coupling of the central DBV-DBV pair that explains the attenuation of coherence beatings observed for these PBPs. We compare two different theoretical strategies in order to derive the excitonic Hamiltonian of these systems, which mainly differ in the way the effect due to the internal geometries of the pigments is handled. Our first strategy is based on classical molecular dynamics (MD) simulations of the complexes, [15] and a post-processing of the sampled structures using polarizable QM/MMPol calculations [16] of the relevant excited states and corresponding electronic coupling values. This strategy was shown to provide a good description of the properties of the similar phycoerythrin PE545 complex, although it required a different shift in order to correct for systematic errors in the site energy predictions for the phycoerythrobilin (PEB) and 15,16dihydrobiliverdin (DBV) pigments in the complex, probably owing to the limitations of the QM method used and the underlying force field adopted in the MD. [10,16] A similar MD-QMMM strategy has been used by the Kleinekathöfer group to study the properties of the PE545 and PE555 antennae [17,18] We then explore a cost-effective strategy based on QM/MM geometry optimization of the pigments in the protein scaffold, followed by multiscale three-layer QM/MMPol/ddCOSMO [19] calculations of the excited states and electronic coupling values, in which the surrounding solvent, absent in the crystal, is modeled through the novel linear-scaling domain decomposition solution of the COSMO equations (ddCOSMO). [20] Recently, the Coker group has shown that combining an intermediate strategy based on multiple geometry optimizations of the pigments along structures sampled from MD can provide a reliable description of the spectra in PE545 and PC645. [21] Our results show that the underlying force field used in the classical MDs introduce significant errors in the site energies, and significant systematic errors between DBVs and the other PCB and MBV pigments are found from these simulations. On the other hand, QM geometry optimizations provide more reliable relative energies between PCB, MBV and DBV pigment types, but the neglect of thermal effects leads to an underestimation of the energies in the MBV pigments in PC630 and PC645. We also compute the spectral densities of excitonphonon coupling for the bilins in the complexes using the Vertical Gradient (VG) approach, [21] based on a normal mode analysis of the pigment vibrations, which lead to excellent emission lineshapes for the majority of the pigments. Finally, we compare electronic couplings and solvent screening factors computed from the MD-QM/MMPol and the QM/MMPol/ddCOSMO methods. Our results indicate a drastic screening of the coupling in the central pair of DBVs in "open" form PBPs, which together with the increased DBV-DBV separation explains the attenuation of coherence beatings observed for these PBPs. Moreover, our results suggest that these PBPs share a common ordering PCB82 < PCB158 < DBV51/61 for the site energies of pigments located in the highlyconserved β chain. The spectral variability among PBPs thus mostly arises from the energetic location of the bilins in the more divergent α chain, which in PC577 and PC612 contribute to the middle-energy region of the spectra (PCB20A and PCB20C), whereas in PC630 and PC645 they contribute to the low-energy one (MBV19A and MBV19C).

MD simulations and geometry optimizations
The structural models and the details of the MD simulations are reported in Ref. [22]. These are based on crystal structures solved for PC645 (PDB 4LMS, resolution 1.35 Å), [14] PC612 (PDB 4LM6, resolution 1.70 Å), [14] PC577 (resolution 1.00 Å) and PC630 (resolution 1.60 Å). PC577 and PC612 are organized as (αβ)2 homodimers and contain 6 PCBs and 2 DBVs, whereas PC630 and PC645 are organized as α1βα2β heterodimers and contain 2 MBVs, 4 PCBs and 2 DBVs. Note that DBVs, in contrast to the other singly-linked pigments, are linked through two Cys residues to the protein. Based on the analysis in Ref. [22], all amino acids were considered in their standard protonation state except His22A and Glu65C in PC630 and His21A and Glu65C in PC645, which were considered protonated, whereas bilin pigments were modeled with a fully protonated tetrapyrrole backbone and anionic propionic side chains. 10-ns MD simulations reported in Ref. [22] were extended to a total of 200 ns, and we extracted a total of 100 snapshots at regular time intervals from the last 50 ns to be used as input structures in QM/MMpol calculations. All MD runs were performed with the Amber 14 suite of programs. [23] Geometry optimizations of the bilins in their native protein environment were performed using the ONIOM method with electrostatic embedding, [24] using B3LYP/6-31G(d) and a classical Amber description for the QM and MM regions, respectively, as implemented in the Gaussian code. [25] The MM region was described using the same force field adopted in the MD. We first performed optimizations, starting from the crystal structures, including the full bilin pigments and selected amino acid interacting with them in the QM region, which was fully relaxed keeping the MM region frozen. In particular, for PCBs and DBVs we included the Asp or Glu side chains coordinating the central rings, whereas for MBVs we included selected amino acids (Lys22A and Lys23A for MBV19A and His22C and Glu26C for MBV19C in PC630; Asn22A for MBV19A and His21C and Glu25C for MBV19C in PC645). QM/MM boundaries were defined at the residue-residue and Cys-bilin bonds using the link atom scheme. [26] We then further optimized the bilin structures limiting the QM region only to the the bilin tetrapyrrole backbone, and thus keeping all amino acids and the propionic side chains frozen in the MM region (capping the pyrroles with methyl groups) in order to limit the NMA vibrational analysis to the tetrapyrrole backbone and in order to prevent contamination of excited-state calculations from charge-transfer effects.

QM/MMPol and QM/MMPol/ddCOSMO calculations
Excited-state energies and electronic couplings were computed using the multiscale two-layer QM/MMPol [27] and the three-layer QM/MMPol/ddCOSMO [20] models. In both models, a QM description of the chromophores is complemented with a polarizable MM description of the environment based on point charges and induced dipoles. In the QM/MMPol/ddCOSMO, however, the solvent is described as a continuum medium characterized by its macroscopic dielectric constant. The polarization equations that stem from the continuum solvation model were solved using the recently developed linear-scaling domain decomposition algorithm for COSMO. [28][29][30] Both models account for full mutual polarization effects among the two or three layers. The QM/MMPol model is based on the following effective Hamiltonian: ,-+ ,-+ + " )) #+ + " )) ,-+ + " ()/../01)0 + " ))/../01)0 (2) We use the extension of both approaches to linear response time-dependent density functional theory (TD-DFT) to perform excited-state calculations. The corresponding transition densities are then used to compute electronic couplings as a sum of Coulomb and environment-mediated terms: [2,31] = /-3+ + #45 (3) The #45 , which in the present formulation includes contributions mediated by both the MMPol and the ddCOSMO environment, typically leads to an overall attenuation of the coupling, and an effective screening factor can be defined as: [32,33] where ε eff is the effective optical dielectric constant of the environment, defined in analogy to the Forster expression of the coupling. QM/MMPol calculations along MD trajectories and QM/MMPol/ddCOSMO calculations based on the optimized geometries of the pigments described in the previous section were computed at the TD-DFT CAM-B3LYP/6-31G(d) level of theory. [34] For the sake of comparison, QM/MMPol/ddCOSMO site energies were also computed using the M06, M06-2X and wB97XD functionals. QM/MM boundaries in the bilin-S bonds were treated using the link atom scheme. [26] The MMPol region was described using the Amber pol12 AL parameters. [35,36] Atomic charges for water were computed from RESP fits at the MP2/aug-cc-pVTZ level of theory on a TIP3P geometry, whereas polarization-consistent ESP charges for bilins were derived at the B3LYP/aug-cc-pVTZ level, based on the crystal geometries, using the Polchat tool. [37] Explicit polarization in the MM region was only included for atoms within a cutoff radius of 18Å from the heavy atoms of the QM region. All calculations were performed using a locally modified development version of the Gaussian package. [25] In this work, the magnitude of the transition dipole moments predicted for the different bilins was in the range ~12-15 D (see Tables S3-S10 in the Supporting Information). This is only ~10-30% larger than an experimental estimate of 11.24 D, [38] verifying that our approach gives correct magnitudes for transition densities and coupling values.

Spectral densities and modeling of steady-state spectra
Spectral densities of the exciton-phonon coupling were obtained with the Vertical Gradient (VG) approach. [39] This method assumes that the ground and excited-state potential energy surfaces (PES) are described by the same harmonic function, but with displaced equilibrium position, and allows extracting the reorganization energy from the difference between the vertical and adiabatic excitation energies. Normal modes were computed based on the optimized geometries described in Section 2.1. Then, vertical gradients were calculated from the vertical excitation energies obtained at the same level of theory. Huang-Rhys factors H of each mode k were obtained by projecting the gradient of the excited state J on the ground-state normal mode k: thus obtaining the following expression for the spectral density: where H is the frequency of the mode k. This model allows estimating the intramolecular part of the spectral density. The continuous part due to slow environmental motions was added a posteriori using an overdamped Brownian oscillator, defined by λ and τ. While τ was fixed to a value of 50 fs, λ was fitted for each protein in order to match the experimental Stokes shift between absorption and emission. Simulations of absorption (OD), circular dichroism (CD) and fluorescence (FL) spectra were performed using the site energies, electronic couplings and spectral densities, as reported in Ref. [10]. Static disorder was modeled by averaging realizations of the spectra over a random distribution of the site energies characterized by a given standard deviation σ, which was adjusted in order to reproduce the broadening of the experimental emission lineshapes for each complex. Lifetime broadening was accounted for based on transfer rates computed using the modified Redfield tensor. [40] We used the following excitonic Hamiltonian describing the multichromophoric system: where N is the number of interacting chromophores, 4 is the excitation energy of the n th chromophore and 4` is the coupling between the n th and m th chromophores. Exciton states | ⟩ were then obtained from diagonalization of the excitonic Hamiltonian: where H is the energy of the k th exciton state and 4 H describes the participation of the n th chromophore to the k th exciton state.

Spectral densities
In Figures S1-S4 of the Supporting Information we show the spectral densities (SD) computed using the VG method for the four PBPs considered, whereas in Tables S1 and S2 of the Supporting Information we report the corresponding reorganization energies. On the other hand, in Figure 2, we show the total accumulated reorganization energies as a function of the vibrational frequency. In the development of quantitative models of light harvesting, it is common to assume equal SDs for all pigments in a given complex. This can be a strong approximation in systems with different kind of pigments, whose coupling to the environment can be different, as is the case of PBPs. [41] Our results clearly indicate a larger coupling of DBVs to the vibrations compared to PCBs. This could be related to the fact that DBVs are covalently linked to the protein backbone through two Cys residues instead of one. The comparison with MBVs is less straightforward, as our simulations provide some clear outliers on the general trends observed in the four complexes, characterized by total reorganization energies in the range λ~410-690 cm -1 : the two DBVs in PC612 and the MBV19A in PC630, which display values λ ~900-1300 cm -1 (see Tables S1 and S2). Also the value for MBV19C in PC645 is somewhat low, with a λ = 350 cm -1 . This suggests that in a few particular conformations of the pigments, as optimized in the protein scaffold, the vertical gradients computed suffer significant errors. Indeed, our calculation on the DBV50/61B pigment in PC630 failed to provide reasonable results. Notably, the main outliers correspond to the lower-resolution structures of PC612 and PC630, indicating that the VG approach is sensitive to inaccuracies in the starting structure.
In Figure 2, one can observe that such outliers arise mostly from inaccuracies in the contribution to λ of low frequency vibrations, probably contaminated by the constrained optimization in the protein scaffold. For DBV50/61B and MBV19A, however, unusual steps in λ are also observed at frequencies ~1500-1600 cm -1 , the region characterized by C=C stretchings. Nevertheless, as we show in the next section, in general the computed SDs lead to simulated emission lineshapes in excellent agreement with experiment, once we use the SD of MBV19C in place of the MBV19A outlier pointed out above. Note that in the spectral simulations of PC612 we used the SD values for DBVs computed for PC577, but this does not impact the emission spectra of this complex as the DBVs have the higher energies in the complex. Therefore, although the calculated SDs cannot account for the low-frequency part due to the environmental motions, they provide an efficient and reliable method to estimate vibronic couplings compared to much more costly MD-QMMM approaches, which usually lead to a strong overestimation of high-frequency peaks when the MD is based on an approximate classical force field. [15,17,[41][42][43][44][45][46][47][48][49][50] A more reliable description consists in performing the MD at the QM/MM level instead of a purely MM potential, [51][52][53] however, the computational cost associated to the calculation of the SD is drastically increased if such a strategy is used.

Site energies
In Figure 3 we report the relative site energies for the bilins in the PBP complexes computed using MD-averaged QM/MMPol values compared with QM/MMPol/ddCOSMO calculations performed on the optimized geometries of the pigments in the protein scaffold. If we focus on PC577 and PC612, which only contain PCB and DBV chromophores, the results show a clear overestimation of the energy difference between these two kinds of molecules when the calculations are based on structures extracted from the MD, whereas the use of QM optimized geometries leads to better relative site energies. Similar systematic deviations between DBV and PEB pigments were observed previously when applying the same approach to the PE545 PBP complex. [10,16] It seems clear that this problem arises from the use of an approximate classical force field in the MD, which is expected to lead to different systematic errors for different kind of pigments. For PC630 and PC645, this issue is further complicated by the presence of MBVs. Regarding DBVs, again the MD-based values are largely overestimated compared to the PCB energies as found for PC577 and PC612. On the other hand, MBVs have the largest degree of π conjugation along the tetrapyrrole backbone, and could thus be expected to be the lowest-lying pigments in these complexes. The use of QM optimized geometries indeed leads to such a result, whereas MD-averaged values lead to MBV energies in ranges similar to those of PCBs, which could be seemingly ascribed again to systematic errors arising from the use of MD-derived geometries. Overall, thus, the force field underlying the MD QM/MMPol protocol seems to be rather problematic in order to derive energy differences for pigments with diverse underlying chemical structure, whereas the use of QM-optimized geometries looks more reliable. This latter approach, however, suffers from the neglect of thermal effects. Indeed, whereas the relative energies between DBVs and PCBs computed in this way are rather accurate, those of the MBVs are clearly underestimated in light of the experimental absorption spectra. Previous calculations based on the crystal structure of PC645, however, suggested that MBVs lie in between the energies of DBVs and PCBs, [54] and QM/MMPol calculations performed on crystal geometries lead to similar results (data not shown). Thus, it seems that thermal effects can be particularly important to properly describe the MBV properties because of its largest degree of π conjugation among all bilin types. In other words, thermal fluctuations are expected to distort the π conjugation that is otherwise exaggerated based on the optimized geometry. Indeed, recent calculations presented by the Coker group showed similar energies for MBVs and PCBs when QMoptimizations were done along different conformations of the protein complex sampled along an MD. [21] Thus, whereas it seems clear that DBVs populate the high-energy part of the absorption spectra of these complexes, it is difficult from the present simulations to determine the precise relative positions of MBV and PCB bands.  . Experimental [14,[54][55][56][57] and simulated spectra of a) open (PC577, PC612), and b) closed (PC630, PC645) cryptophyte antenna complexes. Simulations were done using MD-averaged site energies and couplings computed from QM/MMPol calculations. OD, FL and CD spectra were shifted to match experimental spectra (PC577 1400 cm -1 , PC612 1200 cm -1 , PC630 1400 cm -1 and PC645 900 cm -1 ). Static disorder was modeled by applying a s = 100 cm -1 for PC612, PC630 and PC645, and s = 150 cm -1 for PC577, whereas l for the continuous part of the spectral density was adjusted to match the experimental Stokes shift (PC577 130 cm -1 , PC612 170 cm -1 , PC630 and PC645 75 cm -1 ).
In Figures 4 and 5 we show the absorption (OD), fluorescence (FL) and circular dichroism (CD) spectra of the PBP complexes simulated using both theoretical approaches compared to experiments. [14,55] As previously discussed, the OD spectra in Figure 4, based on MD-averaged site energies and couplings, clearly show that DBV site energies are overestimated. In addition, whereas the relative energies of the PCBs in PC577 are rather accurate, those for the other PBPs seem to span a range of energies that is too small, leading to a spectrum that is too narrow. On the other hand, the FL lineshapes are nicely reproduced by the SDs discussed in the previous section, as emission is not much affected by the relative site energies given that only populated lowest-lying states contribute to it. The main features of the CD spectra of PC612 and PC645 are also reproduced, although the high energy bands display very little intensities.
We remind here that PC577 and PC612, as well as PC630 and PC645, share the same pigment composition, thus spectral changes arise from changes in the amino acid sequence of the protein complex. In this light, the simulations in Figure 4 are in qualitative agreement with the displacement of intensity from the high-energy band to the low-energy one when passing from PC577 to PC612, which seems to arise from a lowering of the energies of PCB20A and PCB20C, located in the α polypeptide chains, whose structure is much more divergent compared to that of the highly-conserved β chains, identical in these complexes. A similar qualitative agreement is found if we compare the spectra of PC630 and PC645. In this case, the lowenergy band is shifted to lower energies in PC645, a trend also qualitatively reproduced by our simulations, which in this case seems to arise from a stabilization of the MBV19C as well as PCB82B and PCB82D. This suggests that the description of pigment-protein interactions based on the QM/MMPol method correctly reproduces these qualitative trends, so the main errors in the determination of the site energies seem to arise from the internal pigment geometries used for excited-state calculations.
For the spectral simulations based on QM-optimized geometries shown in Figure 5, the results for PC612, and especially for PC577, are in good agreement with experiments, although in this case the intensity shift to lower energies when passing from PC577 to PC612 is not qualitatively reproduced, because the energies of PCB20A and PCB20C are not lowered in PC612. Instead, the energies of PCB158B and PCB158D are slightly increased. Overall, our MD results indicate an approximate energy ladder PCB158 < PCB82 < PCB20 < DBV51/61 for PC577 and PC612, whereas results based on QM-optimizations exchange the lowest energy sites for the two PCB82s, and moreover predict a different ordering for the intermediate energy sites, leading to an ordering PCB82 < PCB158 < PCB20 < DBV51/61 for PC577 and PCB82 < PCB20 < PCB158 < DBV51/61 for PC612. Thus, it is unclear if the lowest-energy sites in these complexes are the PCB158 or the PCB82 bilins. The excellent results regarding the overall broadening of OD spectra and the slightly better shape of PC612 CD spectrum indicates nevertheless that results based on QM-optimized structures are more reliable.
On the other hand, in the spectra of PC630 and PC645 the lowest-lying pigments are the MBVs when QM-optimized geometries are used, whereas in MD-based results they show very similar energies compared to the low-energy PCBs, as found by Coker and co-workers for PC645. [21] Moreover, in this case the QM-based approach indicates a similar energy ordering MBV19 < PCB82 < PCB158 < DBV51/61 for both complexes, whereas MD-based data estimate rather similar energies for all low-energy sites: MBV19 ≈ PCB158 ≈ PCB82 < DBV51/61. Although the spectra based on QM-optimized structures are too broad, probably owing to the significant underestimation of MBV energies, the overall shape of the bands due to PCBs and DBVs is in much better agreement than that obtained from MD-based data. Thus, the energy ordering MBV19 ≈ PCB82 < PCB158 < DBV51/61 seems quite likely, once MBVs are shifted up in energy. Moreover, the CD spectrum of PC645 derived from QMoptimized structures is in remarkable agreement with experiment, reproducing the high-energy band that is otherwise missing in MD-based data.
The results on the four antenna complexes thus indicate that calculations based on QM-optimized geometries are remarkably more robust and accurate compared to those derived from classical MD simulations due to the limitations of the force field adopted in the latter. Taken globally, the more reliable QMbased results describe a common energetic ordering PCB82 < PCB158 < DBV51/61 for sites located in β chains, a finding supported by the fact that these chains are highly conserved. The spectral variability thus seems to arise from the energetic location of the bilins in the more divergent α chain, which in PC577 and PC612 seems to be in the middle-energy range (PCB20A and PCB20C), whereas in PC630 and PC645 seem to be in the low-energy region (MBV19A and MBV19C). Indeed, the β chains in PC577 and PC612 are identical. This overall picture further agrees with the recent simulations in the Coker group. [[21]] We further investigated the robustness of our QM-based findings against the choice of DFT functional by recomputing the site energies of PC577 and PC645 using different global and rangeseparated hybrid functionals with different degrees of exact exchange (M06, M06-2X and wB97XD). In Figure S5 of the Supporting Information we compare the OD spectra simulated with these functionals to that obtained using CAM-B3LYP, whereas in Tables S19 and S20 we report the corresponding site energy values. The results indicate that the overall range and ordering of site energies is rather unaffected by the functional choice, thus reinforcing our conclusions obtained with CAM-B3LYP.

Electronic couplings
The spectra of "open" (PC577, PC612) and "closed" form (PC630, PC645) PBPs discussed in the previous section have a fundamental difference related to the change in quaternary structure. In PC630 and PC645, the small center-to-center distance (~16Å) between central DBVs leads to a strong excitonic splitting of the high-energy bands, whereas the switch to the "open" form in PC577 and PC612 leads to a larger separation (~22Å) that drastically attenuates the coupling between DBVs, thus leading to mostly localized bands. [14] According to previous CIS/cc-pVTZ coupling estimates based on the transition density cube method, this leads to a change in DBV50/61B -DBV50/61D coupling from 647 cm -1 in PC645 to 29 cm -1 in PC612. [14] The change in quaternary structure also induces variations in other interdimer couplings, whereas intradimer ones remain similar. Note here that these values do not include the screening effect exerted by the environment, i.e., they only include the Coulomb term in Eq. 3. Our Coulomb terms based on the crystal structures lead to similar values 515/559 cm -1 for PC630/PC645 and 19/21 cm -1 for PC577/PC612, as reported in Tables S15-S18 in the Supporting Information. The MD-based values shown in Table S11-S14 are rather similar, except for PC577, where thermal effects increase the /-3+ value to 67 cm -1 . This is due to a ~2.5Å shorter average interpigment distance along the MD compared to the crystal, but also because structural fluctuations lead to skewed distribution of coupling values with a tail toward large coupling values over 100 cm -1 . The Coulomb coupling, however, is modulated by the additional term #45 of Eq. 3, which describes the interaction between the chromophores mediated by the environment. This term usually counteracts the /-3+ term, leading to a screening of the interaction. In Förster theory, for example, this effect is accounted for using a constant 1/n 2 factor, where n is the refractive index of the medium. Here, instead, we include this effect based on the MD-QM/MMPol and QM/MMPol/ddCOSMO methods, thus explicitly taking into account the heterogeneous nature of the protein environment. The latter can significantly modulate screening effects compared to Förster homogeneous assumption. [16,32] Moreover, in some situations the #45 term can even lead to an enhancement of the coupling depending on the relative arrangement between chromophores. [58] It is worth noting here that the results obtained with the QM/MMPol/ddCOSMO method differ from the MD-QM/MMPol values in two important aspects. First, the water solvent is described through the ddCOSMO continuum model instead of the atomistic MMPol description used for the protein. In addition, the coupling values are based on the crystal structure, so the protein is not relaxed like in the MD, a delicate issue that can lead to some exaggerated polarization interactions and thus induce some noise in the MMPol coupling term. For these two reasons, we thus consider the MD-based values to be more reliable and we only focus on those here. For PC630 and PC645, we predict screening factors equal to 0.60 and 0.50, respectively, for the central DBV pair, slightly smaller than the 0.64 value predicted previously for PC645 based on a continuum model. [59] In contrast, the environment leads to a striking attenuation of the coupling in PC577 and PC612, which screening factors 0.1 and -0.3, respectively, which correspond to /-3+ = 67 cm -1 and #45 = -54 cm -1 for PC577 and /-3+ = 23 cm -1 and #45 = -30 cm -1 for PC612. Thus, for PC612 the environment-mediated term is actually larger in magnitude than the Coulomb interaction between the DBVs. This strong screening effect thus further reduces the coupling in the central DBV pair beyond the increase in pigment separation, which explains the attenuation of coherence beatings observed for these PBPs. [14] Beyond this unexpected strong screening effect, most other pairs in the PBPs experience screening values closer to ~0.5, as expected in a protein environment if the optical dielectric constant is taken as 2, as shown in Figure 6. However, some pairs in PC577 and PC612 deviate from the general s~0.5-0.7 trend with s values close to 1 or even larger, showing a small screening effect or even an enhancement of the coupling mediated by the environment.

Conclusions
In this contribution we have provided an analysis of the spectral properties of four PBPs from cryptophyte algae, PC577, PC612, PC630 and PC645. We compared two different theoretical strategies to derive the excitonic Hamiltonian of the system. The first is based on QM/MMPol excited state calculations performed on geometries sampled from classical MD simulations. We then explored a cost-effective strategy based on multiscale QM/MMPol/ddCOSMO excited state calculations performed on geometries optimized for the bilins in the protein scaffolds. Our results indicate that the force field used in the classical MDs introduces significant systematic errors in site energies predicted for the different type of bilins DBV, PCB and MBV. The second strategy based on QM-optimized geometries gives more reliable energy differences, which are rather unaffected by the choice of functional in the TD-DFT calculations, although the neglect of thermal effects leads to a significant underestimation of the MBV energies. In addition, we compute the spectral densities of the pigments using the Vertical Gradient approach, which leads to emission lineshapes in excellent agreement with experiment. Regarding electronic interactions in the complexes, our results indicate a drastic screening of the coupling in the central pair of DBVs in "open" form PBPs, which together with the increased DBV-DBV separation explains the attenuation of coherence beatings observed for these PBPs. Overall, our results strongly suggest that these PBPs share a common ordering PCB82 < PCB158 < DBV51/61 of the site energies for pigments located in the highly-conserved β chain. The spectral variability among PBPs thus seems to arise from the energy tuning of the bilins in the more divergent α chain, which in PC577 and PC612 are expected to contribute to the middle-energy region of the spectra