Four Resonance Structures Elucidate Double-Bond Isomerisation of a Biological Chromophore

Photoinduced double-bond isomerisation of the chromophore of photoactive yellow protein (PYP) is highly sensitive to chromophore-protein interactions. On the basis of high-level ab initio calculations, using the XMCQDPT2 method, we scrutinise the effect of the chromophore-protein hydrogen bonds on the photophysical and photochemical properties of the chromophore. We identify four resonance structures – two closed-shell and two biradicaloid – that elucidate the electronic structure of the ground and first excited states involved in the isomerisation process. Changing the relative energies of the resonance structures by hydrogen-bonding interactions tunes all photochemical properties of the chromophore in an interdependent manner. Our study sheds new light on the role of the chromophore electronic structure in tuning in photosensors and fluorescent proteins.

is highly sensitive to chromophore-protein interactions. On the basis of high-level ab initio calculations, we scrutinise the effect of the chromophore-protein hydrogen bonds on the photophysical and photochemical properties of the chromophore. We identify four resonance structures -two closed-shell and two biradicaloid -that elucidate the electronic structure of the ground and first excited states involved in the isomerisation process. Changing the relative energies of the resonance structures by hydrogen-bonding interactions tunes all photochemical properties of the chromophore in an interdependent manner. Our study sheds new light on the role of the chromophore electronic structure in tuning in photosensors and fluorescent proteins.

INTRODUCTION
Photoactive yellow protein (PYP) 1 is a remarkable model system for studying double-bond isomerisation in photoreceptor proteins, as this photoreceptor has amply been characterised using a wide variety of advanced methods such as time-resolved X-ray crystallography, 2-5 neutron crystallography, 6 NMR, 7 and ultrafast spectroscopy. [8][9][10][11][12] These studies provided important insights into photoinduced double-bond isomerisation of the PYP chromophore that triggers the PYP photoresponse. The chromophore is derived from the anionic phenolate form of the p-coumaric thioester (pCTM − ), which is in the E(trans)-configuration and stabilised by hydrogen bonds (Hbonds) with the protein (Fig. 1a). 13 The excited-state dynamics of the pCTM − anion in the protein is different from that in solution and the gas phase. [14][15] Despite numerous studies, the exact mechanism of the protein control of the chromophore's isomerisation remains elusive. In particular, it is not clear how to describe the impact of the protein on the chromophore at the electronic structure level.
The moderate molecular size and high sensitivity to intermolecular interactions render the pCTM − chromophore (model 1 in Fig. 1b) particularly suitable for accurate quantum-chemistry calculations addressing the chromophore's tuning by intermolecular interactions. [16][17][18][19][20][21][22] Similarly to other ionic chromophores of photoresponsive proteins, pCTM − undergoes a pronounced chargetransfer, from the phenolic moiety to the carbonyl fragment, upon photoexcitation. 23 Twisting the chromophore's central double bond (C7=C8) or single bond (C4−C7) results in localisation of the negative charge either on the phenolic or carbonyl moieties, respectively, in the excited (S1) state. 17,20,[24][25] The localisation of the negative charge in the ground state (S0) is reversed to that in the S1 state at the twisted geometries. Polar interactions of pCTM − with the protein or solvent profoundly influence the ratio of the S1-state trajectories returning to the S0 state via the single-and double-bond twists. 16,18,24,26 Experimentally it has been established that PYP photoactivation is triggered by the double bond isomerisation yielding a Z(cis) isomer with a quantum yield of 0.35, 27 which proceeds via a one-bond flip (OBF) mechanism destabilizing the carbonyl H-bond. [2][3]9 Activation of the single-bond rotation has been reported for a mutated variant with weakened phenolic H-bonds. 3 The mutation leads to a slight red shift of the absorption maximum, 28 reduced charge-transfer character of the S0−S1 transition, 23 and a so-called hula-twist (HT) isomerisation combining rotations around the single and double bonds. 3 In contrast to the OBF isomerisation, the HT isomerisation preserves the carbonyl H-bond in the resulting cis photoproduct, which impairs the signal transduction. 3 Here, by performing quantum-chemistry calculations for models 1-3 ( Fig. 1b) we study the impact of H-bonds on the S0 and S1 energy at the geometries of pCTM − constituting the single and double bond isomerisation pathways, including the geometries of the S0/S1 conical intersections (CoIns). The accurate description of the models is attained by using a high-level ab initio method, the extended multi-state multi-configurational quasi-degenerate perturbation theory of second order (XMCQDPT2), 29 which accounts for both static and dynamical electron correlation. In this context, the present work critically extends the previous studies performed using the completeactive-space self-consistent field (CASSCF) method 21,24,30 and the single-reference second-order approximate coupled-cluster (CC2) method. 17,20,22,25 To explain the trends revealed from contrasting the three models, we deduce four resonance structures, two closed-shell (CS) and two biradicaloid (BR), that describe charge delocalisation and charge transfer upon the S0−S1 excitation and the bond twists. Stabilisation and destabilisation of these structures by H-bonds account for interdependent changes of all photochemical properties of the pCTM − chromophore. We compare our findings to previously published results on PYP, fluorescent proteins and rhodopsin photoreceptors noting multiple parallels among the properties of ionic chromophore in these proteins. Figure 1. (a) Chromophore binding pocket (active site) of PYP featuring the anionic E(trans)-pcoumaroyl chromophore bound to PYP via a thioester linkage of the Cys69 residue and hydrogenbonded to the backbone of Cys69 and side chains of Tyr42 and Glu46. (b) Three models considered here: the anionic p-coumaric thiomethyl (pCTM -) chromophore (1), pCTMhydrogen-bonded with one water molecule at the carbonyl oxygen (2), mimicking the H-bond with the backbone of Cys69, and pCTMhydrogen-bonded with two water molecules at the phenolic oxygen (3), mimicking the H-bonds with the Tyr42/Glu46 residues.

Computational details
The energies and geometries for models 1-3 in the S0 and S1 states were obtained from the XMCQDPT2 calculations. As a zero-order wave function employed in the XMCQDPT2 calculations, we used a state-averaged CASSCF wave function, SA2-CASSCF (12,11). The active space comprised 12 electrons over 11 π molecular orbitals (MOs), as previously used by Boggio-Pasqua and Groenhof. 21 The state-averaging was done for the two lowest electronic states, S0 and S1, with equal weights. The subsequent XMCQDPT2 calculations using the SA2-CASSCF (12,11) zero-order wave function further accounted for dynamical electron correlation, up to the second order of perturbation theory, and interaction of the S0 and S1 states, via their mixing in a 2×2 Hamiltonian. The resolution of the identity (RI) approximation [31][32][33] implemented in Firefly version 8.2 afforded a substantial speed-up of the XMCQDPT2 computations and a decrease of memory requirements, significantly facilitating numerical evaluation of the energy gradients utilizing supercomputer resources. An intruder state avoidance (ISA) energy shift 34 of 0.02 Hartree was applied throughout. In all our calculations, the correlation consistent double-zeta basis set of Dunning (cc-pVDZ) was used, [35][36] resulting in 231 to 277 MOs depending on the model. All MOs were included in the MCSCF procedure, whereas the chemical-core MOs (1s of carbon and oxygen and 1s, 2s and 2p of sulphur) were not correlated in the XMCQDPT2 calculations.
The search of stationary points on the S1 and S0 PESs and points of S1/S0 minimum energy conical intersections (CoIn) was performed using numerical XMCQDPT2 energy gradients computed with finite differencing of second order. In the search, internal coordinates were used without imposing constraints, i.e, the full geometry optimisation was performed. The root mean square (RMS) gradient for optimisation convergence was 10 -3 Hartree/Bohr. In the CoIn geometry optimisation, a penalty function method was used, with a penalty function as proposed by Martinez and co-workers. 37 The threshold of the energy gap in the CoIn optimisation was 10 -3 Hartree.
Because of exceedingly high computational demands for numerical evaluation of the Hessians at the XMCQDPT2 level, vibrational analysis at the stationary points was not performed. The type of the stationary points was inferred based on previously published computations. 17,20 For the saddle points, the existence of a negative curvature was confirmed by calculating relaxed energy scans along the coordinate for which the negative curvature was expected. The planar structures of models 1 and 2 (E-S1-Sad) and the single-bond twisted structure of model 3 (α-S1-Sad) were demonstrated to correspond to saddle points by computing relaxed-energy scans along the singlebond twisting coordinate (α-torsion). The double-bond twisted barrier (β-S1) in all the models was demonstrated to correspond to a saddle point by computing relaxed-energy scans along the doublebond twisting coordinate (β-torsion).
For all the stationary points, the basis set superposition error (BSSE) corrections were computed for the S0 and S1 energies of models 2 and 3 using the protocol described in Ref 38 . Overall, the BSSE corrections for the relative/excitation energies are small, not exceeding 0.1 eV, with the largest correction (-0.085 eV) found for the α-S1-Sad saddle point of model 3. Hereby, the corrections for the twisted structures are larger than that for the planar structures, which obviously relates to the different charge distribution for the planar and twisted configurations. Importantly, accounting for the BSSE corrections does not lead to changes in the relative order of the stationary points, either within a particular model or the different models. In the following, we therefore neglect the BSSE corrections, presenting and discussing the plain XMCQDPT2 results.
Oscillator strengths, Mulliken atomic charges and free valences on atoms (number of unpaired electrons) 39 were computed using the S0 and S1 zero-order QDPT electron densities. For analysis of the S0 and S1 PESs in the vicinity of CoIn, an approach suggested by Olivucci and co-workers was used. [40][41] A series of points circumscribing a CoIn point (loops) lying in the branching plane were computed. The branching plane is defined by the gradient difference (x1) and derivative coupling (x2) vectors. 42 We computed the x1 and x2 vectors with the SA2-CASSCF(12,11) method using Gaussian09, 43 with solving the coupled-perturbed equations in the MCSCF cycle being skipped ("nocpmcscf" approximation). 44 The error introduced by the latter decreases with the decrease of the S0−S1 energy gap; therefore it might be important at the Frank-Condon region and negligible at CoIn. Except computations of the branching vectors, all calculations were performed using version 8.2 of the Firefly quantum chemistry package, 45 which is partially based on the Gamess-US source code. 46 The computations utilised resources of the high performance computing cluster JURECA. 47  minimum (E-S0-Min). For each model, the cross-section is composed of several stationary points on the S0 and S1 PESs and points of S1/S0 minimum-energy conical intersections (CoIns).

Characterisation of single-and double-bond twists
Specifically, there are the E(trans) and Z(cis) planar minima in the ground S0 state, two minima and two saddle points of planar and twisted geometries in the S1 state, and two twisted CoIn points.
All the points were obtained from a full geometry optimisation using the XMCQDPT2 method. In the following, we refer to the rotation around the C4−C7 bond as α-twist and around the C7−C8 bond as β-twist ( Fig. 1b) according to the previously adopted notation. 17,20,25 The geometries of the models at the stationary points are depicted in Fig. S1 in the Supporting Information (SI). The energies given in Fig. 2 are not corrected for the basis set superposition error (BSSE), which according to our estimation introduces only marginal changes in the cross-section energies (Table S3 in SI).
At first, we consider the results for model 1, which relate to the intrinsic properties of the pCTM − chromophore (blue cross-section in Fig. 2). At the trans configuration E-S0-Min minimum, the S0−S1 vertical excitation energy (VEE) is 2.60 eV and the corresponding oscillator strength is close to unity From the Franck-Condon (FC) point, the chromophore first undergoes a bond-length alteration (BLA) geometry relaxation toward the planar E-S1-Sad saddle point followed by the αtwist relaxation toward the 90°-twisted α-S1-Min minimum. Structurally close to the α-S1 minimum but higher in energy, there is a S1/S0 α-CoIn, mediating decay of the excited-state population back to the S0 state. The excited-state α-twist was previously implicated in internal conversion. 20,22,24 Activation of the β-twist involves passing a small barrier, β-S1-Sad, which is 0.12 eV above VEE, and is associated with the change of the BLA coordinate and β torsion.
Beyond the barrier, the chromophore relaxes by undergoing further β-twist toward the 90°-twisted β-S1-Min minimum. Structurally close to the β-S1 minimum but higher in energy, there is a S1/S0 β-CoIn, through which the chromophore decays back to the S0 state either forward to the Z-S0-Min minimum or backward to the E-S0 minimum. The ratio of the Z and E products depends on many factors, one of which is the topology of the S0 PES in the vicinity of β-CoIn. Both the α-and β-CoIns demonstrate a sloped topology, as each of them was found in the vicinity of a 90°-twisted S1 energy minimum.
The results of models 2 and 3 (green and red cross-sections, respectively, in Fig. 2) demonstrate that the H-bonds lead to considerable energy changes and, moreover, to a change of the S1 PES topology in the case of model 3. Specifically, the carbonyl H-bond in 2 reduces VEE, stabilises the α-pathway and destabilises the β-pathway. In contrast, the phenolic H-bonds in 3 increase VEE, destabilise the α-twist and stabilises the β-twist. Notably, the α-twist energies are more affected by the H-bonds than that of the β-twist. In 3, destabilisation of the α-twist results in appearance of the planar E-S1-Min minimum, with the α-S1 α-twisted geometry becoming a saddle point (red cross-section in Fig. 2).

Four resonance structures rationalizing effect of H-bonds
In the previous studies, the effect of the H-bonds on the S0 and S1 energies has been explained by stabilisation and destabilisation of the negative charge either on the phenolic or carbonyl groups of pCTM − . 17,20,[24][25] At the planar E-S0 minimum, the charge de-localisation is explained by the resonance of two closed-shell (CS) anionic forms. 17,25 At the same time, the origin of the reversed negative charge localisation at the α-and β-twisted S1 minima has not been linked to any resonance forms so far. 17,20,25 To elucidate this issue, we analyzed distribution of the negative charge and the number of unpaired electrons (NUE) 39 at the different geometries (stationary points) in the S0 and S1 state (Fig. 3, Tables S5-S7 in SI). The analysis indicates that the S0−S1 electronic excitation results in a significant charge transfer and increase of NUE, which is consistent with a biradicaloid (BR) structure of the S1 state (Table S8 in SI). The key result of our analysis is realizing that the changes of the phenolic charge and NUE at all the stationary points ( Fig. 3) can be accounted for by changing contributions of four resonance structures − two CS structures, characterizing the S0 state as mentioned above, and two BR structures, characterizing the S1 state (Scheme 1). The resonance structures are denoted according to the presence of the phenolic (Ph) or quinonic (Qn) π-system. Noteworthy, QnBR relates to PhCS, and PhBR relates to QnCS via an electron excitation. Precisely, QnBR describes a single-electron charge-transfer excitation of PhCS, whereas PhBR describes a single-electron charge-transfer excitation of QnCS. To the best of our knowledge, the QnBR and PhBR structures have never been explicitly considered for pCTM − before.
It is worth noting that the four resonance structures are consistent with those introduced in the biradicaloid theory of Bonačić-Koutecký et al. [48][49] and theory of twisted intermolecular chargetransfer (TICT) of Rettig and co-workers. [50][51] The TICT theory treats the "hole-pair" (hp) configuration of a polar chromophore as electron-acceptor and electron-donor moieties connected via a central single bond. The electron transfer from the donor to the acceptor upon the HOMO-LUMO excitation (typically corresponding to the S0−S1 transition) yields a "dot-dot" (dd) structure with a translocated molecular charge. By twisting the linking single bond, the HOMO and LUMO involved in the charge-transfer transition become localised on the mutually orthogonal molecular fragments, and thus do not interact. Therefore, the twisting decreases the energy of the dd configuration and increases the energy of the hp configuration, which altogether significantly reduces the S0−S1 energy difference and even leads to S0/S1 state crossing. 48, 51 Figure 3. Distribution of (a) the net negative charge on the phenolic moiety in the S0 and S1 states and (b) the number of unpaired electrons (NUE) at the C4, C7 and C8 atoms in the S1 state in models 1-3 at the different stationary points along the αand β-twisting pathways. The colour code is the same as in Fig. 2. The atomic and net charges on the phenolic and carbonyl moieties, as well as NUE are given in Tables S5-S7 in SI.
It is worth noting that the four resonance structures are consistent with those introduced in the biradicaloid theory of Bonačić-Koutecký et al. [48][49] and theory of twisted intermolecular chargetransfer (TICT) of Rettig and co-workers. [50][51] The TICT theory treats the "hole-pair" (hp) configuration of a polar chromophore as electron-acceptor and electron-donor moieties connected via a central single bond. The electron transfer from the donor to the acceptor upon the HOMO-LUMO excitation (typically corresponding to the S0−S1 transition) yields a "dot-dot" (dd) structure with a translocated molecular charge. By twisting the linking single bond, the HOMO and LUMO involved in the charge-transfer transition become localised on the mutually orthogonal molecular fragments, and thus do not interact. Therefore, the twisting decreases the energy of the dd configuration and increases the energy of the hp configuration, which altogether significantly reduces the S0−S1 energy difference and even leads to S0/S1 state crossing. 48,51 In Scheme 1, PhCS and QnCS correspond to conjugated hp configurations featuring the phenolic and carbonyl negatively charged electron-donor moieties, respectively. Electron excitation of the PhCS and QnCS hp configurations populates the QnBR and PhBR dd configurations, respectively.
At the same time, the α-(β-) twist leads to the QnBR (PhBR) structure in the S1 state, and the PhCS (QnCS) structure in the S0 state. This picture is fully confirmed by Fig 4, which displays the difference of the S1 and S0 total electron densities at the α-S1 and β-S1 twisted minima. As can be clearly seen in the figure, the negative charge at the α-S1 minimum resides on the alkene-carbonyl fragment in the S1 state and on the phenolic fragment in the S0 state (Fig. 4a), whereas it is vice versa at the β-S1 minimum (Fig. 4b). The same holds for models 2 and 3. . Difference of the S1 and S0 electron densities at the α-twisted (a) and β-twisted (b) minima of pCTM − (model 1). Orange and purple correspond to S1 and S0 electron densities, respectively.
Resonance of the PhCS and QnCS structures determines the single-and double-bond lengths at the planar E-S0-Min minimum and the charge distribution in the S0 state. Figure 5 shows the BLA of the C4−C7−C8 methine bridge in the anionic models in comparison to that in the protonated phenol and enol isomers of pCTM. The shorter C7−C8 distance in comparison to C4−C7 (Fig. 5) and the amount of the negative charge on the phenolate (Fig. 3a) indicate the dominance of PhCS over QnCS in the three models. The prevailing phenolic character of the anionic models suggests that PhCS has a lower energy than QnCS in all the models under consideration. The decrease of the phenolic charge and increase of the C7−C8 bond length in the order 3-1-2 indicate the increase of the QnCS weight, and hence the increase of the PhCS and QnCS mixing in the same order. At the FC point in the S1 state, QnBR dominates over PhBR as it follows from the low negative charge on the phenolate (Fig. 3a) and higher NUE at C4 and C7 than at C8 (Fig. 3b). The weight of QnBR increases in the order 2-1-3, as indicated by a simultaneous decrease of the C8 NUE and phenolic negative charge in this order. In the same order, the S0−S1 charge-transfer character increases, as demonstrated by the increase of the S1−S0 phenolate charge difference (Fig. 3a) and the permanent dipole moment difference (Table S1 in SI). The increase of VEE in the same order of models (Fig. 2) is explained by the higher VEE of PhCS in comparison to that of QnCS.
Thus, mixing of PhCS and QnCS determines de-localisation of the molecular charge in the S0 state, the extent of BLA at the S0 equilibrium geometry, the S0−S1 VEE and the extent of the charge transfer upon S0−S1 excitation. Obviously, the more PhCS prevails in the S0 state (small PhCS/QnCS mixing), the larger the BLA and phenolic charge in the S0 state and the larger the S0−S1 VEE and the extent of the charge-transfer. In contrast, the less the PhCS prevails in the S0 wave function (large PhCS/QnCS mixing), the smaller the BLA and phenolic charge in the S0 state, and the smaller the S0−S1 VEE and the charge-transfer upon excitation.

Energies of the resonance structures − impact on the S1 state and S1/S0 conical intersections
The energies of the resonance structures can be estimated for the 90°-twisted geometries. 52 Indeed, the four resonance structures are uncoupled at the α-S1 and β-S1 twisted minima, with the negative charge being completely localised either on the phenolic or carbonyl fragment, as evident in Fig   3a. At the α-S1 minimum, the S0 state corresponds to the PhBR structure, whereas the S1 state is represented by the QnBR structure, with the negative charge being localised on the phenolic moiety and carbonyl fragment in the S0 and S1 states, respectively ( Fig. 3a and 4a), and the C8 NUE has the smallest value (Fig. 3b). At the β-S1 minimum, the S0 state is represented by the QnCS structure, while the S1 state is associated with the QnBR structure, with the negative charge being localised on the carbonyl fragment and phenolic moiety in the S0 and S1 states, respectively ( Fig. 3a and 4b), and the C8 NUE has the largest value (Fig. 3b).
By comparing the energies at the 90°-twisted geometries (Fig. 2), we infer that PhCS has the lowest energy among the four resonance structures in all the models, which is in agreement with the analysis of the resonance mixing at the planar geometries described above. The carbonyl H-bond stabilises QnCS and QnBR and destabilises PhCS and PhBR, whereas the phenolic H-bonds stabilise PhCS and PhBR and destabilise QnCS and QnBR. Among the four resonance forms, the energy of QnBR demonstrates the largest variation due to the H-bonding interactions.
As the geometries of the stationary points demonstrate, the BLA, α-twist and β-twist are the main coordinates driving the relaxation of the pCTM − chromophore in the S1 excited state. If the QnBR energy is lower than VEE, the initial relaxation from the FC point involves the BLA coordinate followed by the α-twist. If that condition is not fulfilled, which is the case for model 3 (see Fig.2), the planar E-S1 minimum appears on the S1 PES.
At the planar E-S1-Min minimum (model 3), the S0−S1 oscillator strength is 0.9, and the computed radiative rate is ~0.26 ns -1 , suggesting a sizable emission with a Stokes shift defined by the extent of the BLA relaxation. At the twisted α-S1 minimum in models 1 and 2, the computed radiative rate is several orders of magnitude smaller, ~0.51 s -1 , because of the substantial decrease of the S0−S1 energy and oscillator strength. Based on these results, we predict that the dynamics of the fluorescence signal in PYP correlates with the dynamics of the chromophore H-bonds stabilizing/destabilizing the planar E-S1 minimum. The smaller BLA variation at the planar S1 minimum as compared to that at the S0 minimum is a consequence of the chromophore instability with respect to the α-twist. As the α-twist virtually quenches the fluorescence, the variation of the emission maximum in colour-tuned PYP mutants should be significantly smaller than that of the absorption maximum, which indeed was experimentally observed by Hoff and co-workers. 28 Finally, the energies of the resonance structures at the α-S1 and β-S1 minima determine the energies of α-CoIn and β-CoIn, which in turn define whether the excited chromophore returns to the S0 state via a single-or double-bond twist. As Fig 2 demonstrates, the CoIn energies correlate with the S0−S1 energy gap at the 90°-twisted geometries, i.e., the difference in energy of the corresponding CS and BR resonance structures. The smaller the energy gap at the 90°-twisted S1 minimum (or saddle point), the lower the CoIn energy and the smaller the distortion of the CoIn geometry in comparison to the twisted minimum (or saddle point). As an example, we consider model 3, for which the PhCS and QnBR energy difference is rather high, and consistently the highenergy α-CoIn in this model demonstrates a substantial β-twist in addition to the α-twist (Fig. S1 in SI). In contrast, stabilisation of PhBR and destabilisation of QnCS by the phenolic H-bonds lowers the energy of the β-CoIn and reduces its distortion with respect to the β-S1 minimum. In models 1 and 2, the PhCS and QnBR energy difference is smaller than in model 3, the α-CoIn energy therefore decreases, whereas the PhBR and QnCS energy difference is larger, and hence the β-CoIn energy increases. The variations of the α-CoIn energy among the models suggests that accessibility of α-CoIn and therefore activation of the single bond rotation and internal conversion via the single-bond twist can be efficiently controlled by chromophore-protein interactions.
Stabilisation of PhCS and PhBR by these interactions reduces the resonance mixing and favours the double-bond OBF isomerisation.

Double-bond isomerisation of the pCTM − chromophore
The β-CoIn provides the excited-state decay channel for the OBF double-bond isomerisation of pCTM − . We remind that a minimum energy CoIn geometry corresponds to a local energy minimum on the 3N-8 dimensional hypersurface (N is the number of atoms) where the S1 and S0 states are degenerate. 42 The degeneracy is lifted in the first order along the gradient difference vector (x1) and the derivative coupling vector (x2), forming a so-called branching plane. [53][54][55] Inspection of the energies in the branching plane unveils relaxation on the S0 PES after a nonadiabatic S1-S0 transition at the β-CoIn. We computed x1 and x2 at the β-CoIn geometries and constructed branching-plane loops following the procedure suggested by Olivucci and coworkers. 56 In all our models, x1 is dominated by BLA, whereas x2 features the β-torsion and, in addition, hydrogen out-of-plane (HOOP) 20,22,57 excursions, which mix with the BLA changes in models 1 and 2 (see Fig. S2 in SI). At the β-CoIn geometries, the HOOP value tends to be "ahead" of β-torsion, i.e., smaller than β-torsion, when β > 90° (in models 1 and 2), whereas it is "concerted" with β-torsion (close to β-torsion) when β < 90° (model 3) (Fig. S3 in SI). Fig. 6a presents the x1 and x2 vectors for model 3 (for the other models see Fig. S2 in SI). Fig. 6b demonstrates that β-CoIn has the β-torsion value slightly smaller than 90° in model 3, whereas it shifts toward a larger β-torsion value (i.e., closer to the initial E(trans) geometry) in models 1 and

The smaller β-torsion value at the β-CoIn, explained by stabilisation of PhBR and destabilisation of QnCS in model 3, is favourable for the E(trans)-Z(cis) isomerisation. The β-CoIn branching
plane of pCTM − is similar to that previously reported for the double-bond isomerisation of the retinal protonated Schiff base (RPSB) chromophore and analogues in rhodopsins. [58][59] Figure 7 shows the energy and property changes along the branching plane loops. In the x1 direction (θ=90°, 270°), the S0 and S1 wave functions change from one resonance structure to another. Specifically, the charge distribution and NUE (Fig. 7b) indicate that the S0 wave function changes from QnCS to PhBR whereas the S1 wave function changes from PhBR to QnCS. The x2 direction (θ=0°,180°) corresponds to mixing of PhBR and QnCS and, to a lesser extent, of PhCS, as indicated by the changes of the phenolic charge and NUE (Fig. 7b). The presence of PhCS follows from the steeper phenolic charge curve in comparison to the NUE curve. The PhCS weight significantly increases as the chromophore becomes more planar.
Two S0 energy minima emerge along x2 (Fig. 7a) as decay channels towards the Z(cis) and E(trans) minima. Notably, the Z(cis) decay channel lies below the β-CoIn energy only in 3 (at θ=0°, Fig. 7a), correlating with the smaller β-torsion value at β-CoIn in this model. In addition, there is a high-energy barrier, associated with the QnCS structure, separating the two channels in model 3, which is explained by the higher QnCS energy in this model. Overall, the topology of the S0 PES along the loops suggests that the double-bond isomerisation quantum yield is higher in 3 as compared to 1 and 2. Figure 6. Characterisation of the β-twisted CoIn: (a) The branching plane vectors x1 and x2 for model 3 (see also Fig. S2 presenting x1 and x2 for all three models) (b) Changes in the β-torsion and C7=C8 bond length along the 0.005-Å loops for models 1-3. The values of the angular coordinate θ at some points of the loops, and the orientation of the x1 and x2 vectors are indicated. The colour code is the same as in Fig. 2.
Finally, we discuss on some aspects related to the sloped character 60 of the β-CoIn. There are two stationary points associated with the β-CoIn, the β-S1 minimum on the S1 PES, and a saddle point on the S0 PES, which are of the PhBR and QnCS electronic structure, respectively. The S0 saddle point is associated with the aforementioned QnCS energy barrier separating Z(cis) and E(trans) channel (Fig. 7a). Our preliminary calculations for model 3 have identified such a point. This is the only saddle point connecting the Z(cis) and E(trans) minima on the S0 PES, and hence it mediates the thermal isomerisation. The sloped type of β-CoIn implies that the aforementioned two stationary points lie on the same side from the β-CoIn along x1. This is different from the peaked CoIn found for the RPSB chromophore and its analogues, where two saddle points of different electronic structure lying on the S0 PES flank the CoIn along the BLA coordinate. [61][62][63][64] Noteworthy, for the rhodopsin photoreceptor, the energy of the saddle point controlling the thermal isomerisation was found to correlate with VEE, such that the blue-shifted absorption maximum corresponds to a reduced thermal activation (higher energy of the saddle points). 61 Despite different β-CoIn topology, a similar correlation as for the RPSB chromophore should be expected for the pCTM − chromophore, as the increase of VEE correlates with the increase of the QnCS energy.

Comparison to other computational studies
Although our models 1-3 do not account for the real environment of the pCTM − chromophore in PYP, they correctly grasp the effects of the Cys69 and Tyr42/Glu46 H-bonds (Fig. 1a). The energy changes presented in Fig. 2 are qualitatively consistent with those predicted by the hybrid quantum-mechanics/molecular mechanics (QM/MM) models of PYP 16,[18][19]26 and solvated chromophore. 24 Specifically, the carbonyl H-bond lowers the α-CoIn energy, whereas the phenolic H-bonds lower the β-CoIn energy in the QM/MM models. 24 In the surface-hopping molecular dynamics simulations of PYP, 26 the presence of the protonated (positively charged) Arg52 stabilizing the phenolic negative charge, that is the PhCS and PhBR resonance structures with respect to the QnCS and QnBR ones, leads to stabilisation of the β-CoIn, and hence facilitates the OBF isomerisation.
We compare the energies of model 1 with the results of Boggio-Pasqua and Groenhof 21 obtained using the CASSCF (12,11) method, and the energies of model 3 with the results of Gromov et al. 17,20 obtained using the CC2 method. It is instructive to consider the energies of the resonance structures at the 90°-twisted geometries. Both CC2 and CASSCF favour QnBR over PhBR more than predicted by the XMCQDPT2 method. Consistent with this energy trend, the β-CoIn is higher in energy than α-CoIn at the CASSCF level of theory, whereas it is vice versa at the XMCQDPT2 level. Thus, employing the CASSCF method to describe the pCTM − chromophore in surfacehopping dynamics of PYP 16,26 underestimates the efficiency of the β-twist in the S1 state. For a similar anionic chromophore of the green fluorescent proteins (GFPs), it has indeed been reported that the CASSCF method underestimates the decay of the S1 state via the β-twist as compared to the more accurate XMSCASPT2 method. 65 In the QM/MM simulations of PYP, 26 this underestimation seems to be compensated by the effect of the positively charged residues, such as protonated Arg52 stabilizing, the PhCS and PhBR structures.

Similarity with other biological chromophores
Similar to pCTM − , the anionic p-hydroxybenzylideneimidazolinole (pHBDI − ) chromophore of GFPs is described by resonance interactions of the PhCS and QnCS forms coupled along the BLA reaction coordinate. [66][67] Recently, Boxer and co-workers 67 suggested that the energy difference of the PhCS and QnCS forms can serve as a linear scale for quantitative prediction of the GFP photophysical properties using a small number of parameters intrinsic to the pHBDI − chromophore. This energy difference was implicated to explain all strong correlations of properties among a large number of systematically tuned GFP variants. 67 We note that the charge-transfer GFP model of Boxer and co-workers 67 is consistent with our computational results for the anionic pCTM − except that it does not invoke the QnBR and PhBR structures of pHBDI − . By analogy to pCTM − , we suggest that the energies and properties of the QnBR and PhBR resonance forms should be considered when addressing the S1 PES of pHBDI − in GFPs, in particular its photoactivation. Applying this prediction, stabilisation of PhBR in pHBDI − favours the OBF double-bond isomerisation, whereas stabilisation of QnBR gives a chance to the single-bond activation and HT isomerisation. 68 Mixing of the resonance forms describes charge delocalisation in protonated cationic chromophores such as the RPSB chromophore of rhodopsins and tetrapyrrole chromophores of phytochromes. The extent of mixing explains different sensitivities of these chromophores to electrostatic interactions with their proteins. The dominant resonance structure of RPSB features a positive charge localised on the Schiff-base fragment, giving rise to considerable charge-transfer character of the S0−S1 transition and high sensitivity to electrostatic interactions with the protein.
These interactions in turn, efficiently control the contributions of the resonance forms, as suggested by the correlations of the increased/decreased BLA and increased/decreased VEE. [69][70][71] In contrast, the tetrapyrrole positive charge is redistributed among the four pyrrole rings described by at least four resonance forms. 72 This charge delocalisation reduces the charge-transfer character of the S0−S1 transition and subsequently reduces the electrostatic effect of the protein on the chromophore.

CONCLUSIONS AND OUTLOOK
Given the prominent role of ionic chromophores in light activation of naturally evolved photoreceptor proteins, it is important to rationalise the effect of chromophore-protein interactions that are crucial for the light sensory function. Here, by analysing how the anionic pCTM − chromophore of the PYP bacterial photoreceptor is influenced by phenolic or carbonyl H-bonds we derived a qualitative model explaining interdependent changes of all photochemical properties related to the PYP photoactivation. To account for electronic correlation effects, the state-of-theart XMCQDPT2 method was employed to compute the S0 and S1 PES cross-sections. To rationalise the obtained trends, two pairs of resonance structures, each pair comprised of a CS structure and corresponding to it charge-transfer BR structure were invoked (Scheme 1).
Delocalisation of the negative charge and other properties of the S0 and S1 states are described by at the S0 minimum-energy geometry; there is also an increase of the excitation energy and stabilisation of the planar S1 minimum. Moreover, isomerisation around the central double bond becomes more favourable than the phenolic ring rotation around the single bond, and the energy barrier controlling the thermal double-bond isomerisation increases. In contrast, stabilisation of the quinonic structures by the carbonyl H-bond leads to delocalisation of the negative charge and a smaller BLA in the S0 state, smaller S0−S1 excitation energy and smaller charge transfer character, as well as activation of the single-bond rotation in the S1, increasing the probability of the HT photoisomerisation. Overall, our computational results predict that stabilisation of the pCTM − phenolic resonance structures, i.e. reduction of chemical resonance, is crucial for PYP photoactivation via the OBF double-bond isomerisation.
Comparison of our findings and conclusions to the previously published computational and theoretical results characterizing other biological chromophores identifies many common features.
In particular, our model derived from the high-level quantum-chemistry calculations in its essence is similar to a model proposed by Boxer and co-workers 67, 73 explaining spectral properties of systematically mutated GFP proteins. Yet, our model extends Boxer's model to the twisted geometries and invokes BR structures to describe the S1 state properties. The broad similarities among the ionic chromophores imply that the concept of chemical resonance is instructive in experimental and especially computational studies addressing photochemical mechanisms of photosensory proteins operating via the double-bond isomerisation.  Table S.1: Relative S 0 and S 1 energies (eV) and differences in the permanent dipole moments of the S 1 and S 0 states (∆µ = µ S 1 − µ S 0 , in Debye) at the different stationary points and points of S 1 /S 0 minimum energy conical intersections for models 1, 2 and 3. In parentheses, S 0 →S 1 oscillator strengths.  Table S.3: Basis set superposition error (BSSE) corrections to the S 0 and S 1 total and relative (in parentheses) energies (eV) at the different stationary points for models 2 and 3.   Table S.5: Net atomic charges (charges at the adjacent hydrogens are included) in the S 0 and S 1 states at the different stationary points for models 1 (first number), 2 (second number) and 3 (third number). For the atom numbers see Figure 1b of the paper.  Table S.6: S 0 and S 1 charge distributions at the different stationary points in terms of the net charges on the phenolic (first number) and carbonyl (second number) moieties of the chromophore for models 1, 2 and 3. The charge on the water molecule(s) is not included.  Table S.7: Number of unpaired electrons at the chromophore atoms in the S 0 and S 1 states at the different stationary points for models 1 (first number) 2 (second number) and 3 (third number). For the atom numbers see Figure 1b of the paper.