Conversion of light-energy into molecular strain in the photocycle of the photoactive yellow protein †

The Photoactive Yellow Protein (PYP) is a light-driven photoreceptor, responsible for the phototaxis of halophilic bacteria. Recently, a new short-lived intermediate (pR0) was characterized in the PYP photocycle using combined time-resolved X-ray crystallography and density functional theory calculations. The pR0 species was identified as a highly contorted cis-intermediate, which is stabilized by hydrogen bonds with protein residues. Here we show by hybrid quantum mechanics/classical mechanics (QM/MM) molecular dynamics simulations, and first-principles calculations of optical properties, that the optical shifts in the early steps of the PYP photocycle originate from the conversion of light energy into molecular strain, stored in the pR0 state, and its relaxation in subsequent reaction steps. Our calculations quantitatively reproduce experimental data, which enables us to identify molecular origins of the optical shifts. Our combined approach suggests that the short-lived pR0 intermediate stores B1/3 of the photon energy as molecular strain, thus providing the thermodynamic driving force for later conformational changes in the protein.


Introduction
The photoactive yellow protein (PYP) is a small (14-kDa) watersoluble protein found in purple photosynthetic bacteria that live in halophilic environments. 1PYP functions as a blue-light receptor, in which the light-absorption triggers a signaling cascade that switches the rotational direction of the bacterial flagellum, and directs the bacteria to swim away from potential harmful high energy light. 2,3PYP has become a key photobiological model system for studying light-absorption induced changes in proteins, due to its small size but relatively complex photocycle. 4The photoactive center of PYP comprises a p-coumaric acid (pCA) chromophore that is covalently linked to a cysteine residue (Cys-69) via a thioester bond (Fig. 1).The pCA is stabilized by hydrogen bonds from Tyr-42, Glu-46 and Thr-50, whereas the backbone of Cys-69 stabilizes the chromophore in early intermediates of the photocycle.Arg-52 shields the chromophore from the surroundings, and may electrostatically stabilize the anionic chromophore.[7][8][9] Recently, Schotte et al. 10 characterized a new photocycle intermediate using picosecond time-resolved X-ray crystallography and density functional theory (DFT) calculations. 10,11This new short-lived state, pR 0 , appeared on a 150 ps timescale (Fig. 2B), and was characterized as a highly strained cis-species with the double bond ca.301 out of plane.3][14][15] Recently Jung et al. 16 also characterized a new short-lived species, I T , forming on time-scales similar to the pR 0 state.In contrast to Schotte et al., 10 however, the refined species was suggested to be a transition state structure between pG and later intermediates, with the C2QC3 double bond refined to ca. 901 out of plane.Jung et al. also used DFT calculations to study the likelihood of their species, but in contrast to their refined X-ray structure, the obtained DFT-structure was in close structural agreement with the pR 0 species described by Schotte et al. 10 Differences between these two structures were recently discussed by Kaila et al. 11 and Jung et al. 17 However, the focus of this study concerns optical tuning effects in the PYP photocycle studied by QM/MM calculations.
9][20][21][22][23] By stabilizing the electronic ground state or excited state, the absorption can be shifted towards the blue (higher energy) or red (lower energy) regime of the spectrum.The strain energy stored in the recently discovered pR 0 species is expected to lead to a destabilization of the electronic ground state, which in turn would further redshift the optical properties.Moreover, relaxation of the strain in the subsequent photocycle steps would be expected to blueshift the absorption.
To study the structure and dynamics of the early photocycle intermediates in PYP, we use here a combination of QM/MM molecular dynamics simulations, first-principles computation of optical properties at second-order approximate coupled cluster (CC2) and time-dependent density functional (TDDFT) theory levels to elucidate how the photon energy is converted into molecular strain in the PYP photocycle, and to characterize the dynamics of the subsequent intermediates.

Computational methods
Time-resolved X-ray structures of the pG, pR 0 , pR 1 , and pR 2 states (PDB ID: 2ZOH, 4B90, 4BBT and 4BBU) were obtained from the Brookhaven Protein Databank and used as starting points for the QM/MM simulations.We also performed QM/MM simulations on the refined I T structure by Jung et al. 16 (PDB ID: 4I38) to study the X-ray model dependence of our data.The QM region (143 atoms) included the pCA chromophore, Tyr-42, Glu-46, Thr-50, Arg-52, Phe-62, Val-66, Ala-67, Pro-68, Cys-69, Thr-70 and Phe-96, and was described at the D-BP86/ def2-SVP level of theory.The MM region, comprising the surrounding protein residues and water molecules, was described using the CHARMM27 force field. 37Each state was structure optimized, followed by 5 ps QM/MM MD simulations using a time step of 1 fs at 310 K. From the QM/MM trajectories, we extracted 100-500 snapshots for each state along the MD trajectories for calculation of ground state and optical properties.The absorption spectrum of the early PYP photocycle states was computed at the time-dependent density functional theory Fig. 2 Structure of photocycle intermediates.The relative energy (E rel in kcal mol À1 ) stored within the pCA chromophore is also shown.The C 2 QC 3 double bond in pCA is trans-in the pG state and cis in all pR states.This double bond is strained in pR 0 , which upon breaking the hydrogen bond to Cys-69 converts to pR 1 , which is much more planar and less strained.The transition from pR 1 to pR 2 involves a rotation of the C b -S g -C 1 -C 2 dihedral angle, which transforms from syn-syn to anti-anti, and further reduces strain in the pCA.
This journal is © the Owner Societies 2016 (TDDFT) level using CAM-B3LYP/def2-SVP, and at the secondorder approximate coupled cluster level (CC2/def2-TZVP).The first six vertical excitation energies (VEEs) at TDDFT level of theory were computed for each snapshot leading to evaluation of 5000 VEEs per method (10 000 VEEs in total), whereas at the CC2 level, we computed the first two VEEs.To speed up the evaluation of the VEEs at the CC2 level, we employed the reduced virtual space (RVS) approach 38 where virtual orbitals above 70 eV from the highest-occupied molecular orbital (HOMO) were disregarded from the correlation calculations, leading to a RVS 38 cutoff error of 0.02 eV relative to frozen core CC2 calculations (Fig. S1, ESI †).For the isolated chromophore, the CC2 and CAM-B3LYP produced similar VEE profiles, but CAM-B3LYP was found to systematically blueshift the absolute excitation energies relative to CC2 (Fig. S2, ESI †).All CAM-B3LYP profiles were therefore redshifted by a 0.37 eV offset-factor obtained from the CC2 fit.For calculation of optical properties, we reduced the size of the QM region to explicitly include the pCA chromophore, Cys-69, Tyr-42, Glu-46 and Arg-52, for the TDDFT calculations, and the pCA chromophore for CC2 calculations, while the remaining protein and water atoms were described using point charges obtained from the QM/MM trajectories.For the estimation of electrostatic-induced shifts, we calculated the optical spectra of the chromophore along the QM/MM trajectories, by removing all protein residues around the chromophore but by keeping the chromophore structure unchanged.Strain-induced optical shifts were calculated as the difference between the absorption maximum of the isolated (trans or cis) chromophore structures obtained from gas-phase and the protein simulations, without considering the electrostatic polarization of the protein environment.The absorption spectra were calculated as an average of 500 VEEs obtained from the MD trajectories and weighted with their respective oscillatory strengths for each VEE calculation.
The distribution of weighted VEE along the trajectories was fitted with a Gaussian width of 100 hc/E h eV, where E h is the excitation energy in eV, and combined to give a spectra by scaling the maximum intensity to 1 or by matching to the maximum experimental intensities.The two first excited states were included in the calculations due to larger oscillator strengths at TDDFT level.The spectra of each state were also calculated by computing histograms of the observed VEEs along the MD trajectories, and by normalizing the probability distribution.We obtain very similar spectra using both the Gaussian fitting and the histogram approaches, showing that the calculated line shapes arise from the statistical sampling of the ground state dynamics.Both approaches were found to give similar absorption spectra.QM/MM simulations were performed using the Q-Chem/ CHARMM. 37,39TDDFT and CC2 (RVS) calculations were performed with Q-Chem v4.1 40 and TURBOMOLE v6.3, 41 respectively.

Ground state dynamics of PYP photocycle intermediates
Our QM/MM molecular dynamics simulations suggest that all studied photo-cycle intermediates; pG, pR 0 , pR 1 , and pR 2 remain stable on picosecond time-scales.Fig. 2 shows how the energy stored in the dihedral bonds of the chromophore is released during the pR 0 -pR 2 transition, in which the dihedral angle around the C 2 QC 3 double bond relaxes from a B301 strained cisconformation in pR 0 to planar species in pR 1 and pR 2 .The pR 0 state stores ca.20 kcal mol À1 of strain energy relative to the pG state, which corresponds to about 1/3 of the photon energy (446 nm = 64 kcal mol À1 ).The pR 0 state is stabilized by a hydrogen bond between the carbonyl group of pCA and the backbone of Cys-69 (Fig. 2), as suggested by Schotte et al. 10 Table 1 Comparison of geometrical parameters between crystal structure (C) and average (AEstandard deviation) values obtained from 5 ps QM/MM (Q) trajectories of the four photocycle states Consistent with time-resolved X-ray crystallographic experiments, we also observe a transient formation of a planar pR 1 -like species during the trajectory, suggesting that the pR 0 /pR 1 species may reside in (quasi)-equilibrium during turnover (see Movie S1(C and D), ESI †).Although the strain energy of the pR 1 state is lower relative to pR 0 , the intrinsic energy of the chromophores in the two states are similar, due to dissociation of the backbone hydrogen bond between pCA and Cys-69 (Fig. 2).In the pR 2 state, we observe ca.5-10 kcal mol À1 release of strain energy relative pR 1 by further relaxing the C b -S g -C 1 -C 2 dihedral angle (Table 1).We also performed QM/MM MD simulations on the I T structure of Jung et al. 16 to study how the dynamics depends on the employed X-ray model.Consistent with our earlier observations 11 we find that the refined I T structure is unstable and the dihedral angle around the C 2 QC 3 rapidly relaxed from 801 to 341, and stays in this conformation during the 5 ps QM/MM simulation trajectory (Fig. S4, ESI †).The two X-ray models thus result in very similar conformational dynamics.
The hydrogen bonding distances between pCA, and Tyr-42 and Glu-46, are shown in Fig. 3 for the four studied photocycle intermediates.We observe short hydrogen-bonding distances (2.5 Å, Table 1) and partial proton transfer between Tyr-46 and the pCA, which becomes more pronounced in the pR 1 and pR 2 states as strain energy is released from the chromophore.A possible reason is that the proton affinity of the pCA increases with the decrease in strain due to flow of electron density from the double bond to the phenol ring (Fig. S3, ESI †).The hydrogen bond between Glu-46 and the pCA is also short (2.6 Å), but in contrast to Tyr-42, we do not observe proton transfer between the two former residues.
Yamaguchi et al. 42 refined the position of hydrogen (deuterium) atoms in the pG state using neutron diffraction experiments of deuterated samples, and observed an unusually short hydrogen (deuterium) bond of 1.38 Å between the Glu-42 and the pCA.Saito and Ishikita 43,44 studied the energetics of the proton transfer between Glu-42 and pCA by hybrid QM/MM optimizations, and observed short distances between the residues only upon dissociation of Tyr-42 from pCA.They also found a slightly uphill proton transfer barrier between Tyr-42 and the pCA, which seems consistent with our dynamics simulations, where the proton transiently fluctuates towards the pCA (Fig. 3).To study the distance distribution of hydrogen bonds observed in the different photocycle states, we computed the correlation between the proton donor (D -phenolic/carboxylic oxygen of Tyr or Glu) and proton acceptor (A -phenolic oxygen of pCA) distances, with the A-H (acceptor-proton) and D-H (donorproton) distances.This correlation plot, shown in Fig. 4, suggests that short hydrogen bonding distances, indicated by blue dotted lines, are observed for Tyr-42, but not for Glu-46.Our simulations further suggest the short distances are observed both for the D-H (closed circles) and A-H (open circles) bonds when the A-D distance is o2.5 Å.It is, however, possible that the experimentally observed short distances between the pCA and Glu-46 could arise from dynamically averaged fluctuations of multiple positions relative to the static heavy atom positions obtained from the X-ray structure in the neutron diffraction refinement.Such short bonds have, e.g., been observed in high resolution Xray structures of aquaporin. 45tical properties of PYP photocycle intermediates Fig. 5 shows the dynamically averaged absorption spectra computed at the TDDFT and CC2 levels of theory for the pG, pR 0 , pR 1 , and pR 2 states.The computed absorption maxima are compared to experimental spectra, 9 summarized in Table 2.We find an excellent agreement between the experimental and computational data.The expected absolute error of the employed electronic structure methods is around 0.2 eV (or 40 nm in the 450-500 nm regime) based on benchmarking calculations, which may lead to small discrepancies between the first-principles computed and measured spectra.The quantitative agreement allows us to computationally characterize the likely reason for the observed optical shifts.absorption of the trans-pCA at 460 nm. 46The effect is counterbalanced by a 0.08 eV (0.13 eV) blueshift due to electrostatic stabilization of the surrounding protein environment.Formation of the pR 0 intermediate leads to a large strain-induced redshifted of 0.41 eV (0.61 eV), consistently with the destabilization of the electronic ground state, discussed above.The large strain effect is somewhat counter balanced by the blue-shifting effect of 0.16 eV (0.18 eV) by the protein environment, resulting in absorption peak at around 485 nm (483 nm).Relaxation of this strain energy by 0.1-0.2eV leads to a gradual blueshift of the predicted spectra in pR 1 and pR 2 .Interestingly, we find that the partial proton transfer between Tyr-42 and pCA leads to a small population of blueshifted intermediates in the pR 0 and pR 1 states, an effect that becomes more pronounced when the pCA chromophore relaxes.This effect might arise from the increased oscillator strength with decreasing chromophore strain, as suggested by model calculation on the isolated pCA chromophore.
We also computed the absorption spectrum of the 5 ps QM/MM simulation of the I T model by Jung et al. 16 As expected from the similar conformational dynamics (Fig. S4, ESI †), we also obtain nearly identical spectra for the I T and pR 0 trajectories (Fig. S5, ESI †), further supporting that the early 150 ps intermediate leads to large redshift due to strain-induced effects in the pCA.
Our calculations suggest that the molecular strain redshifts all intermediates, while electrostatic polarization has a blueshifting effect on the absorption properties.However, we observe that in contrast to the electrostatic effect, which is 0.08-0.16eV in all intermediates, the main spectral shift of the pCA chromophore originates from the molecular strain effects, which has a maximum in the recently characterized pR 0 state.

Conclusions
In conclusion, we have studied here the structure and dynamics of early photocycle intermediates of PYP, and identified how strain energy, originating from light-absorption is stored in the pCA chromophore using a combination of QM/MM molecular dynamics simulations, and first principles prediction of optical properties.Our data suggest that the recently identified pR 0 state is stable on picosecond timescales, and that the highly strained species is stabilized by hydrogen bonds to the protein environment, which allows the intermediate to store ca.1/3 of the photon energy as molecular strain.Optically, we find that this strained intermediate leads to formation of redshifted species, followed by gradual blueshifting of the spectra with relaxation of the strain.Our simulations also identified short hydrogen bonds and partial proton transfer between Tyr-42 and the pCA chromophore, leading to a small population of blueshifted species in both the computed and experimental optical spectra.While reproducing well the observed optical shifts, our work cannot be used to unambiguously assign the structures responsible for the early steps of the PYP photocycle (ref.10, 11,  16 and 17) due to the intrinsic computational errors involved in the employed computational methodology.Our work, nevertheless, clearly shows how the strain effects contribute to the optical tuning in the early PYP photocycle intermediates, providing a powerful methodology to study elementary energy conversion steps in the photocycle of PYP.

Fig. 1
Fig.1The structure of the pCA chromophore in the PYP dark state (pG).Left: pCA and its hydrogen bonding partners.Right: The structure of PYP showing the chromophore-binding region, modeled at the QM level of theory.

Fig. 3
Fig. 3 Dynamics of hydrogen bonding distances in the pG, pR 0 , pR 1 , and pR 2 states.The proton coordinate r = r 1 (TyrO-H/GluO-H) -r 2 (pCA:OÁ Á ÁH) (in Å), with the Tyr-42 and Glu-46 bonds shown in red and blue, respectively.r o 0 indicates that the proton is closer to pCA than to Glu-42/Tyr-46.Right inset: Probability distributions of r, showing that the Tyr-42 is shorter relative to the Glu-42 bond.

Fig. 6
Fig.6The computed ensemble averaged spectra of the four photocycle intermediates pG, pR 0 , pR 1 , and pR 2 obtained from hybrid QM/MM MD simulations.The absorption maximum of the isolated and protein embedded chromophore are shown in red and blue, respectively, indicating strain (arrow and red values in eV) and electrostatic (arrow and blue values in eV) induced optical shifts.The vertical dotted line shows the averaged absorption maximum for pCA in the trans (pG), and cis (pR 0 , pR 1 , and pR 2 ) conformations.The CAM-B3LYP/def2-SVP spectra are computed as a sum of individual Gaussians, weighted with their respective oscillator strength (solid lines) and as histograms of the VEEs obtained from the MD simulations (vertical bars).The spectra obtained at the RVS/CC2/def2-TZVP level is shown in transparent purple in the background.