Xuchun Yang

Madushanka Manathunga

Samer Gozem

Jérémie Léonard

Tadeusz Andruniów

Massimo Olivucci

The activation of rhodopsin, the light-sensitive G-protein coupled receptor responsible for dim-light vision in vertebrates, is driven by an ultrafast excited state double-bond isomerization with a quantum efficiency of almost 70%. The origin of such a high light sensitivity is not understood. A key unanswered question is whether and how in-phase nuclear motion controls the quantum efficiency value. Here, we employ hundreds of quantum-classical trajectories to show that, 15 femtoseconds after light absorption, a degeneracy between the reactive excited state and a neighboring state, causes the splitting of the rhodopsin population into subpopulations propagating with different velocities and leading to distinct contributions to the quantum efficiency. We also show that such splitting is modulated by the protein electrostatics, thus linking amino acid sequence variations to quantum efficiency modulation. Finally, we discuss how such a linkage that, in principle, could be exploited to achieve higher 2 quantum efficiencies, would simultaneously increase the receptor thermal noise leading to a trade-off that may have played a role in rhodopsin evolution.

1 Quantum-classical simulations of rhodopsin reveal excited-state-population splitting and its effects on quantum efficiency 3 Rhodopsin (Rh) features an 11-cis retinal protonated Schiff base (rPSB11) chromophore covalently bound to a protein (opsin) cavity formed by seven α-helices. [START_REF] Ernst | Microbial and Animal Rhodopsins: Structures, Functions, and Molecular Mechanisms[END_REF] As illustrated in Fig. 1a, absorption of a photon induces a sub-picosecond isomerization of rPSB11 to its all-trans stereoisomer (rPSBAT) which initiates the receptor photocycle and, ultimately, visual transduction. [START_REF] Ernst | Microbial and Animal Rhodopsins: Structures, Functions, and Molecular Mechanisms[END_REF][START_REF] Khorana | Rhodopsin, photoreceptor of the rod cell. An emerging pattern for structure and function[END_REF][START_REF] Peteanu | The first step in vision occurs in femtoseconds: complete blue and red spectral studies[END_REF][START_REF] Polli | Conical intersection dynamics of the primary photoisomerization event in vision[END_REF] The Rh photoisomerization has a low thermal noise [START_REF] Gozem | The Molecular Mechanism of Thermal Noise in Rod Photoreceptors[END_REF] that, when combined with its 0.67 quantum efficiency (Φcis-trans) value [START_REF] Dartnall | The photosensitivities of visual pigments in the presence of hydroxylamine[END_REF] results in an extremely high light sensitivity. The photoisomerization has also been reported to display a coherent (i.e. phased) nuclear motion of the molecular population along the reaction coordinate. [START_REF] Wang | Vibrationally coherent photochemistry in the femtosecond primary event of vision[END_REF][START_REF] Kukura | Structural Observation of the Primary Isomerization in Vision with Femtosecond-Stimulated Raman[END_REF][START_REF] Johnson | Local vibrational coherences drive the primary photochemistry of vision[END_REF][START_REF] Mathies | A coherent picture of vision[END_REF] Such phased motion is attributed to the dynamics induced by the potential energy surface (PES) of the reactive excited state. In other words, the nuclear forces created by the PES slope impose specific phase relationships between the modes describing the rPSB11 transfer from the first singlet excited (S1) to the ground (S0) state. Three rPSB11 modes are implicated in the S1 to S0 transfer (see Fig. 1b). [START_REF] Schapiro | The Ultrafast Photoisomerizations of Rhodopsin and Bathorhodopsin Are Modulated by Bond Length Alternation and HOOP Driven Electronic Effects[END_REF][START_REF] Gozem | Theory and Simulation of the Ultrafast Double-Bond Isomerization of Biological Chromophores[END_REF] These are the C11=C12 counterclockwise twist (α), the hydrogen out-of-plane deformation of the HC11=C12H hydrogens (β) and the skeletal bond length alternation (BLA) stretching of the conjugated backbone. In Fig. 1c we also conveniently define the dihedral angles δop, representing the hydrogen out-of-plane deformation relative to the carbon skeleton, and τ, representing the orbital overlap across the isomerizing bond. As shown in Fig. 1c, the S1 to S0 transfer is described in terms of progression along a multi-mode coordinate connecting the S1 vertical excitation region of the PES to a decay region in the vicinity of the S1/S0 intersection space (ISS1/S0). [START_REF] Gozem | Theory and Simulation of the Ultrafast Double-Bond Isomerization of Biological Chromophores[END_REF] Past results [START_REF] Schnedermann | Evidence for a vibrational phase-dependent isotope effect on the photochemistry of vision[END_REF] have shown that the relationship between the phases of specific modes determines if rPSB11 would proceed towards rPSBAT and form the bathorhodopsin (bathoRh) photoproduct or relax back to its starting configuration. [START_REF] Ernst | Microbial and Animal Rhodopsins: Structures, Functions, and Molecular Mechanisms[END_REF] More specifically, at the singlemolecule level, the relationship between the phases of the velocities of α and δop at decay determines if the isomerization is successful or unsuccessful. [START_REF] Schapiro | The Ultrafast Photoisomerizations of Rhodopsin and Bathorhodopsin Are Modulated by Bond Length Alternation and HOOP Driven Electronic Effects[END_REF][START_REF] Schnedermann | Evidence for a vibrational phase-dependent isotope effect on the photochemistry of vision[END_REF][START_REF] Weingart | Product formation in rhodopsin by fast hydrogen motions[END_REF] The δop change appears to be critical. Its physical meaning implies that, in order to generate bathoRh, the H11/H12 hydrogen rotation described by the angle β (see Fig. 1b) shall decrease faster than the decrease in skeletal isomerization described by α. Only when this condition, yielding a positive δop change, is satisfied, the orbital overlap (proportional to τ) changes in direction thus allowing for the formation of a trans C11=C12 double bond (see ref. 14 for a more detailed explanation). Thus, as illustrated in Fig. 1d, when the C10/C13 carbon and H11/H12 hydrogen rotations are inphase, the τ velocity is negative, and the decay to S0 results in a product-like overlap leading to rPSBAT formation. Since, τ=α-δop/2 and α mainly change in the negative counterclockwise direction (dα/dt<0), a successful isomerization is equivalent to having a positive δop velocity (dδop/dt>0) at decay. The existence of a link between the phases of the α and δop velocities and light sensitivity in Rh has been experimentally demonstrated by showing that deuterium substitutions at H11 and/or H12 modulate Φcis-trans. [START_REF] Schnedermann | Evidence for a vibrational phase-dependent isotope effect on the photochemistry of vision[END_REF] In principle, a population displaying a phased nuclear motion along BLA, α and δop would be able to exploit the mechanism revised above to maximize Φcis-trans. As illustrated in Fig. 1e, this limit would be reached when the largest possible population fraction decays with dδop/dt>0. In order to find out how closely Rh follows such an ideal case, we have employed a hybrid multiconfigurational [START_REF] Sen | A comprehensive benchmark of the XMS-CASPT2 method for the photochemistry of a retinal chromophore model[END_REF][START_REF] Gozem | Mapping the Excited State Potential Energy Surface of a Retinal Chromophore Model with Multireference and Equation-of-Motion Coupled-Cluster Methods[END_REF][START_REF] Andersson | Second-order perturbation theory with a CASSCF reference function[END_REF] quantum mechanical / molecular mechanics (QM/MM) model of Rh [START_REF] Gozem | The Molecular Mechanism of Thermal Noise in Rod Photoreceptors[END_REF][START_REF] Luk | Molecular bases for the selection of the chromophore of animal rhodopsins[END_REF] to study the S1 dynamics of an initially room-temperature S0 population using 200 quantum-classical trajectories. The results show that the nuclear forces generated by the steep S1 PES synchronize the trajectories that initially describe a compact "in phase" progression. However, such a progression starts to break down just 15 fs after photoexcitation through a fractional but still organized process. The analysis of such process in terms of coupled nuclear and electronic changes discloses a theoretical framework that appears instrumental for the comprehension and control of the observed Φcis-trans value. This report is organized as follows. First we describe the result of Rh dynamics simulations showing that: i) population desynchronization/splitting occurs along the α mode due to a mixing of the S1 and S2 electronic states, ii) the Φcis-trans is the result of contributions associated with each generated subpopulation and iii) the splitting is modulated by the opsin electrostatics.

Based on findings i-iii, we then proposed a mechanistic framework for the modulation of Φcistrans in opsin-based receptors. Finally, the same results are employed to hypothesize a possible trade-off between photochemical and thermal isomerizations which may have shaped the Φcistrans value observed in Rh.

Results

S1 population splitting. The validation of the constructed Rh model and the stability of the population dynamics simulation are documented in Supplementary sections 1, 2 and 3. In Fig. 2a we report the evolution of α for the entire trajectory set. It is apparent from inspection of the figure that the S1 to S0 decay occurs when α has a value of -90º±30º (i.e. when the π-bond of C11=C12 is broken) and lasts for ~150 fs starting with an early ~35 fs decay event. The trajectories can be used to compute a number of observables including the S1 lifetime, photoproduct appearance time, and Φcis-trans. Supplementary Table 1 shows that the computed values are consistent with the experimental observations including the relevant S1 vibrational frequencies [START_REF] Johnson | Local vibrational coherences drive the primary photochemistry of vision[END_REF] and the isotopic modulation of Φcis-trans. [START_REF] Schnedermann | Evidence for a vibrational phase-dependent isotope effect on the photochemistry of vision[END_REF] Furthermore, the simulation complies with the relationship linking the dδop/dt sign with successful isomerizations (see Supplementary Fig. 4a) as 73% of the successful trajectories have, at decay, dδop/dt>0.

As discussed above (see Fig. 1e), to achieve top quantum efficiency, Rh would generate a S1 population decaying when dδop/dt>0. However, the data in Fig. 2a show that, along α, the population remains compact for less than 30 fs and then splits. The splitting produces subpopulations decaying at different times forming peaks I-IV in Fig. 2b which feature, as expected, increasing δop values for groups of successful decay events (see Fig. 2c). Peak I corresponds to the fastest subpopulation which decays between 35 and 80 fs (curve a in Fig. 2a). In contrast, the slowest subpopulation inverts the initial counterclockwise to a clockwise rotation pushing α back (curve b in Fig. 2a) and then re-inverts to decay at times ranging from 110 to 170 fs forming peaks III-IV. The occurrence of splitting does not appear to depend on the specific QM/MM model employed or type of visual pigment, although peak height and distribution may differ substantially. Such behavior is supported in Supplementary Fig. 5 that displays data from different QM/MM models.

The results above indicate that Φcis-trans is the sum of fractional contributions coming from different subpopulations. Since each subpopulation hits ISS1/S0 at a different time and δop phase, it would display a different dδop/dt distribution yielding a different Φcis-trans contribution. It is thus apparent that the splitting mechanism, which determines the number of subpopulations and their dynamics, is critical. More specifically, it is apparent that the nuclear forces causing the splitting, as well as their relation to the S1 electronic structure, have to be investigated.

The ultrashort timescale of the splitting event suggests that Franck Condon forces are its primary cause. Such forces are experimentally fingerprinted by the Rh Resonance Raman (RR) spectrum displayed in the top panel of Fig. 2d. The simulation of such spectra may thus help to formulate a correct mechanistic hypothesis. We found that, when only S1 forces computed at the equilibrium Rh structure were employed to simulate the spectrum (see bottom panel of Fig. 2d), certain observed spectral regions were poorly reproduced suggesting the involvement of other electronic states in the population dynamics. Indeed, as shown in the same figure, the RR spectrum simulated using an effective gradient with an increased weight of the S2 gradient with respect to S1 gradient displayed a substantial improvement (see also Supplementary section 10).

Accordingly, we hypothesized that a mixing of the S1 and S2 electronic characters during the initial S1 dynamics is implicated in the population splitting. S2 and S1 mixing and splitting mechanism. The mixing of the S1 and S2 electronic characters was investigated by focusing on two sets of trajectories. As shown in Fig. 3a, the first (49 trajectories) and second (48 trajectories) set represent the fastest and slowest subpopulations respectively. At time 0 the two sets are indistinguishable in terms of α, δop and BLA distributions, indicating a weak dependence of the population splitting on the initial conditions (see Supplementary Fig. 6). However, the progression along BLA and δop reported in Fig. 3b shows that the two sets start to split along BLA after just 15 fs. We could also demonstrate that the splitting is modulated by the electrostatics of the opsin. Indeed, Fig. 3c and3d show that the splitting is reduced when the simulation is repeated after setting all opsin atomic charges to half of their original values (see Supplementary Fig. 8) and, therefore, effectively decreasing the electrostatic field acting on the chromophore. After calculating the full set of 200 trajectories with the half-charge model, we confirmed that the S1 population remains more compact and decays to S0 on a faster timescale. The fact that such a behavior is closer to a single population decay (see Fig. 1e), led us to recompute Φcis-trans that was found to be 0.72 and therefore higher than the 0.68 value obtained for the unaltered Rh model. This suggests that: (i) the splitting decreases the Φcis-trans, value and (ii) such value is not maximized in Rh.

The hypothesis that the mixing of the S1 and S2 electronic characters during the initial S1 relaxation leads to a slower α progression, was initially investigated by computing the evolution of the S2-S1 gap and S1 electronic character along the fast and the slow trajectory sets. In Fig. 3e we show that both the fast (top panel) and slow (bottom panel) sets enter a S2/S1 degeneracy region (i.e. the ISS2/S1 region) in ~10 fs. However, while such region is abandoned in the first case, the slow population dwells in that region for over 50 fs (see also the Supplementary Fig. 9). These differences are reflected in the behavior of the corresponding variations in the chromophore charge distribution and, therefore, electronic character; while the fast trajectories (Fig. 3f, top panel) gradually change their charge distributions (i.e., for instance, the charge of the N-containing backbone moiety of Fig. 4a top-left) in smooth fashion, the trajectories of the slow set (bottom panel) show multiple sudden charge variations pointing at multiple changes of the S1 electronic character consistently with transits close to ISS2/S1. Motion in the region of a S2/S1 conical intersection. The results of Fig. 3e and 3f can be interpreted by recalling that at time 0, the S1 and S2 states of Rh are dominated by 1Bu and 2Ag electronic characters respectively. [START_REF] Gozem | Theory and Simulation of the Ultrafast Double-Bond Isomerization of Biological Chromophores[END_REF] As illustrated in Fig. 4a, the 1Bu character is described by a charge transfer (ionic) configuration (central structure) where the C11=C12 bond acquires single-bond character prompting torsional deformation and reactivity. On the other hand, the 2Ag character is dominated by a double excitation and has a diradical character (bottom structure) featuring a C11=C12 bond with residual double-bond order. Accordingly, the slow set would be associated with 1Bu/2Ag character mixing and inversion restraining double bond twisting. Remarkably, such an inversion is not seen along the trajectory sets computed using the Rh model with halved opsin charges indicating that the change in opsin electrostatics drives the two trajectory sets away from the ISS2/S1 region.

The points belonging to the ISS2/S1 region are represented by the S2/S1 minimum energy conical intersection (MinCIS2/S1) of Fig. 4b which is located below the FC point. From inspection of the MinCIS2/S1 branching plane cross section in Fig. 4c, it is clear that the fast and slow trajectory sets travel on opposite sides of MinCIS2/S1. While the fast set propagates along the side dominated by the 1Bu character, the slow set enters the opposite side featuring a region with a larger 2Ag character and then re-enters a region with 1Bu character. The slow set should therefore transit in a more bounded region where its progression is restrained. Remarkably, both fast and slow trajectory sets propagate on the same side of MinCIS2/S1 when the opsin charges are halved (see Supplementary Fig. 10). In fact, the scaling of the opsin electrostatics leads to the displacement of MinCIS2/S1 to a higher energy location (see inset in Fig. 4c) that is avoided by the majority of the trajectories.

The above conclusions are supported by Mayer's analysis of the S2 and S1 wavefunctions (see Supplementary Fig. 15) showing distinct free valences for the fast set (i.e. at the putative radical centers C9 and C14 the free valence is ca. 0.5 and 0.2 in S2 (2Ag) and in S1 (1Bu) respectively) but not for the slow set (e.g. around 0.3 at both C9 and C14 for both S2 and S1) indicating character mixing.

The presence of 1Bu/2Ag mixing in the slow trajectory set indicates not only increased bonding along the S1 PES, but also the possibility of non-adiabatic transitions. This hypothesis has been investigated by recomputing the initial 50 fs dynamics of the entire population at a level of theory explicitly incorporating the S2 PES (see the Method section). As reported in Fig. 4d-f, when such PES is included in the calculation, hops from S1 to S2 are detected along the trajectories suggesting that non-adiabatic effects may contribute to restrain the motion of the slow set with respect to the fast set (see the Supplementary sections 14-16). Indeed, Fig. 4d show that the trajectories associated with a slower α progression display a larger percentage of sequential S1→S2 and S2→S1 transitions. This result is in line with the hypothesis that larger non-adiabatic coupling (NADC) characterizes the slow α motion. This is confirmed by mapping the NADC modulus along S1 PES cross-sections spanned by the critical modes α, δop and BLA.

The results displayed in Extended Data Fig. 1 in the Supplementary Information shows that, while at the initial α values (i.e. at the beginning of the S1 relaxation) there is no remarkable difference between slow and fast populations, at α = -20 degrees the slow set moves closer to the region with high NADC modulus.

While the results above would point to a non-adiabatic or mixed adiabatic/non-adiabatic origin of the S1 population splitting, a comparative analysis of the initial population dynamics computed without and with S2 PES, points to a very limited difference in geometrical progression and time scale and, therefore, population splitting. This points to a dominant role played by the topography of the adiabatic S1 PES (i.e. its slopes, barriers and paths) which appears not to be qualitatively altered by the employed level of theory. In other words, while it is possible that the splitting involves a mixed adiabatic/non-adiabatic mechanism (i.e. comprising the two types of trajectories depicted in Fig. 4f) the hop events occurring along several trajectories do not qualitatively change the dynamics driven by the S1 PES.

Effect of the protein electrostatics. In order to explain how the scaled-opsin charges modify the S1 PES and, therefore, the trajectory progression and splitting, we have mapped the protein electrostatic potential (ESPopsin). It is apparent that, while ESPopsin has a complex structure (see Fig. 5a), a suitable cross-section (see Fig. 5b) shows that ESPopsin is more negative in the Schiff base rather than β-ionone region. It is thus possible to hypothesize that ESPopsin better stabilizes the 2Ag character relative to 1Bu character leading to a decrease of the S2-S1 energy gap (see Supplementary Fig. 15). In fact, the chromophore S2 and S1 charges in the Schiff base region (see Fig. 5b) show a larger charge in the S2 state indicating 2Ag character. The ESPopsin decrease would thus lead to a destabilization of S2 and displacement of MinCIS2/S1 to higher energies (from -1 to +14 kcal/mol with respect to FC) far from the S1 relaxation path (compare the blue and green circles in Fig. 4c).

In the Supplementary Fig. 15-17 we demonstrate that the described effects of ESPopsin are independent from the employed levels of theory including multistate multiconfiguration levels (XMS-CASPT2 and XMCQDPT2 [START_REF] Granovsky | version 8.0.0[END_REF] ) and expanded basis set (ANO-L-VDZP). To further support our interpretation, XMS-CASPT2/ANO-L-VDZP numerical gradients have been employed to locate the S2 energy minimum (MinS2) in a simplified Rh model. The calculation indicates that MinS2 is located 2 kcal/mol below the FC point (see Supplementary Fig. 18) and it is not far from degeneracy (i.e. 10 kcal/mol S2-S1 gap).

Discussion and Conclusions

The results above connect an opsin-induced S1 population splitting to the Φcis-trans value in a dim-light visual pigment. Such a connection is based on two findings. The first is that Φcis-trans is the sum of contributions associated to different subpopulations and the second is that the splitting generating the subpopulations is governed by the opsin electrostatics. We now discuss the Rh photoisomerization mechanism on the basis of a theoretical framework incorporating these results.

As outlined in Fig. 6a, immediately after vertical S0→S1 excitation the Rh population starts a phased motion along the S1 reaction coordinate. However, the population enters a region with a reduced S2-S1 gap where a set of subpopulations form in 15 fs. Each subpopulation then propagates with different timescale and δop phase. Thus, the observed Φcis-trans value must be the consequence of the synchronization between population splitting at ISS2/S1 and generation of favorable dδop/dt vs. dα/dt phases at ISS1/S0. This suggests that the Φcis-trans should be studied, as mentioned in early works by Warshel et al., [START_REF] Warshel | The dynamics of the primary event in rhodopsins revisited[END_REF] by using simulations accounting for multiple modes and states. In fact, the described mechanism appears to be closely related to the one reported in Lim et al. [START_REF] Lim | Experimental probing of conical intersection dynamics in the photodissociation of thioanisole[END_REF] who have proposed that the behavior of molecules reaching the reactive intersection space (in our case ISS1/S0) is determined by early events (i.e., a sort of memory effect).

According to our calculations the involvement of S2, with its 2Ag character, in the Rh dynamics is limited to the initial relaxation. However, such temporary 2Ag character not only justifies the population splitting, but also its sensitivity to the opsin electrostatic field and, therefore, amino acid sequence. In other words, a variation in the sequence, such as one caused by natural selection, may transform a single population (Fig. 1e) into multiple subpopulations (Fig. 6a) and modify the resulting Φcis-trans value. Notice that, due to the complexity of this process, the details of the QM/MM model alter the number and distribution of the simulated subpopulation, but not the fact that they are generated nor its sensitivity to the opsin electrostatics.

It remains to be seen if the mechanism presented above has a significance for the adaptation of a visual pigment to dim-light conditions. Such a pigment would need maximum lightsensitivity and, therefore, the conclusion that population splitting makes the optimization of Φcis-trans complex or, as suggested by our results, decreases its value, is counterintuitive.

However, light sensitivity in Rh depends on both Φcis-trans and the noise associated with rPSB11 thermal isomerization. As we will now explain, it is possible to argue that population splitting is the byproduct of a trade-off between photochemical and thermally isomerizations leading to a maximum light sensitivity.

A trade-off hypothesis can be formulated starting from the reported [START_REF] Gozem | The Molecular Mechanism of Thermal Noise in Rod Photoreceptors[END_REF] inverse proportionality between λmax and E T (compare the left and right parts of Fig. 6b) where E T is the energy barrier controlling the thermal (S0) isomerization of rPSB11. This inverse proportionality implies that Rh may have evolved to absorb blue light to maximize E T and, in turn, decrease the rate of thermal Rh activation. However, shorter λmax are achieved by selecting an opsin sequence whose electrostatics destabilizes S1 with respect to S0. Since S0 and S2 have similar charge distributions (see Fig. 6c), such destabilization must lead to smaller S2-S1 gaps and, as a consequence, population splitting. This is in line with the documented effects of the ESPopsin scaling causing, simultaneously a decrease in ΔES1-S0 (i.e. a red-shift in λmax) and an increase in ΔES2-S1. The small S2-S1 gap found in Rh would thus be the result of sequence variations designed to decrease its thermal noise but also inducing a counterproductive Φcis-trans decrease due to population splitting. Such a hypothesis remains to be experimentally assessed.

In conclusion, the documented connection between population splitting and Φcis-trans advances the construction of the theoretical framework necessary for the understanding of Rh lightsensitivity and reveals a novel effect of the opsin electrostatics on the primary event in vision. 27- 30 Such connection is not likely to be a common feature across the vast rhodopsin family, even if evidence exists for its involvement in microbial rhodopsins. [START_REF] Gozem | Excited-State Vibronic Dynamics of Bacteriorhodopsin from Two-Dimensional Electronic Photon Echo Spectroscopy and Multiconfigurational Quantum Chemistry[END_REF] On the other hand, the presented isomerization efficiency theory not only goes beyond the seminal one-dimensional Landau-Zener model 32 holding that a high speed along α determines a high Φcis-trans, but provides a distinct and multi-mode scenario for the involvement of the S2 state with respect to past hypotheses. [START_REF] Gai | Chemical dynamics in proteins: The photoisomerization of retinal in bacteriorhodopsin[END_REF] At the same time, it raises important questions about our ability to quantitatively predict the effect of mutations on Φcis-trans as this would require a comprehension of the coupled electron-nuclear dynamics [START_REF] Albert | Coupled electron-nuclear quantum dynamics through and around a conical intersection[END_REF][START_REF] Olivucci | Unlocking the Double Bond in Protonated Schiff Bases by Coherent Superposition of S1 and S2[END_REF] and possible geometric phase effects [START_REF] Xie | Up to a Sign. the Insidious Effects of Energetically Inaccessible Conical Intersections on Unimolecular Reactions[END_REF] [START_REF] Lin | Vibrational Assignment of Torsional Normal Modes of Rhodopsin: Probing Excited-State Isomerization Dynamics along the Reactive C11C12 Torsion Coordinate[END_REF] Bottom panel: Comparison between the resonance Raman spectra simulated using the S1 gradient and using an equal mixture of S1 and S2 gradients according to the weighted-gradients method. [START_REF] Jarzȩcki | Porphyrin Distortion from Resonance Raman Intensities of Out-of-Plane Modes: Computation and Modeling of N-Methylmesoporphyrin, a Ferrochelatase Transition State Analog[END_REF] The arrows point at regions requiring states for being correctly described. for 200 fs followed by a 2-root-state-average CASSCF/6-31G*/Amber level on the S0 PES for 50 fs. The final geometries and velocities from those ground-state QM/MM dynamics comprise the initial conditions for subsequent S1 population dynamics. We note that such initial conditions correspond to the effect of an ultrashort light pulse treated in the so-called "classical" limit of an instantaneous optical transition in the Condon approximation.

We propagated 200 quantum-classical trajectories starting from the 200 initial conditions defined above for 200 fs each at the 2-root-state-average CASSCF/6-31G*/Amber level using the Tully surface-hop method and equal S0 and S1 weights. [START_REF] Granucci | Critical appraisal of the fewest switches algorithm for surface hopping[END_REF] The trajectory yielded a simulation of the dynamics of the entire population (population dynamics). We also rerun all 200 trajectories from the original initial conditions but with the opsin atomic charge to be half of the original value. In order to study the impact of surface hopping on the initial 50 fs dynamics describing the population splitting the full trajectory set has been recomputed at both the 3-root-state-average CASSCF/6-31G*/Amber level with equal S0, S1 and S2 weights as well as with 0.0, 0.5 and 0.5 weights (equivalent to a 2-root-state-average CASSCF/6-31G*/Amber level with equal S1 and S2 weights).

Raman spectra simulations. The Franck Condon factors were obtained within the harmonic approximation assuming identical normal modes and vibrational frequencies in the S0 and electronic excited states. The normal modes and vibrational frequencies of S0 equilibrium geometry were calculated using the QM/MM model of Rh seen above and the CASSCF/6-31G*/Amber level of theory, S1 and S2 gradients were obtained using the same QM/MM model at the 2-root-and 3-root-state-average CASSCF wave function, respectively.

More details of model construction, population dynamics simulations and Raman spectra simulations are provided in the Supplementary Information.

Data availability. The authors declare that the data supporting the findings of this study are available within the main article and the Supplementary Information. Additional data are available from the corresponding author upon request.

  occurring very early during the dynamics of biological chromophores.

Fig. 1 .Fig. 2 .

 12 photoexcited S1 population

Fig. 3 .Fig. 4 .

 34 Fig. 3. Analysis of the subpopulation dynamics. (a) α progression along the fast and slow trajectory sets. The framed region represents the investigated α (-10 to -40 degrees) and time (0-50 fs) window. (b) Progression along the average BLA and δop values for the fast and slow sets. (c) and (d) Same as part A and B but for a Rh model with halved opsin atomic charges. (e) Time evolution of S2 and S1 potential energy profiles for the fast and slow sets. (f) Time evolution of the S2 and S1 charge transfer magnitude for the fast and slow sets. The 1Bu and 2Ag labels refer to the symmetry of the electronic wavefunction of a reference planar all-trans polyene.

Fig. 5 .

 5 Fig. 5. Opsin electrostatic potential at MinCIS2/S1. (a) Protein electrostatic potential map on the solvent accessible surface (solvent radius 1.4 Å) of the Rh chromophore. (b) Protein electrostatic potential crosssection along the conjugated framework. The fraction of S1 and S2 charges on the chromophore Schiff base region (framed moiety in red) are given at the top of the panel. The corresponding FC charges are given in square brackets.

Fig. 6 .

 6 Fig. 6. The mechanism controlling Φcis-trans and its possible origin. (a) Left. Three electronic states and two intersection spaces, ISS2S1 and ISS1S0, are involved in the primary event in vision. ISS2/S1, which is accessed by relaxing along the lower intersecting state, mediates population splitting while ISS1/S0 is accessed from the higher intersecting state and mediates the rPSB11 branching occurring at decay. Right. The relationship between decay time and phase of the hydrogen-out-of-plane wagging determines the Φcistrans contribution of each subpopulation. (b) Trade-off between thermal noise (λmax values) and population splitting. (c) Inverse proportionality between the FC computed S0-S1 and S2-S1 changes upon decrease of the protein charges.

Acknowledgements

We are grateful to Richard A. Mathies for pointing our attention to the Rh resonance Raman spectrum, Belinda Chang for suggesting the idea of a trade-off and Stefan Haacke for helpful discussions. We are also grateful to Nicolas Ferrè for the frequency calculation code in Molcas/Tinker and Leonardo Barneschi with help in generating the plots of Figure 5. The research has been partially supported by the following grants: NSF CHE-CLP-1710191, NIH GM126627-01, USIAS 2015, the Ohio Supercomputer Center, the MIUR (Ministero

Code availability. The authors declare that the present research has been produced with distributed software available to the public as also detailed in the Supplementary Information.

Competing interests Statement

The authors declare no competing financial interests.

Figure Legends

Methods

Rh QM/MM model construction. The Rh QM/MM model was constructed from the 2.2 Å resolution crystallographic structure (PDB ID 1U19) [START_REF] Okada | The Retinal Conformation and its Environment in Rhodopsin in Light of a New 2.2Å Crystal Structure[END_REF] and following a reported protocol. [START_REF] Schnedermann | Evidence for a vibrational phase-dependent isotope effect on the photochemistry of vision[END_REF] The QM layer comprises atoms of the retinylidene chromophore, NH group and CεH3 substituent.

The rest of the protein forms the MM layer described by a modified AMBER94 force field featuring specific Lys296 side-chain parameters. [START_REF] Cornell | A second generation force field for the simulation of proteins, nucleic acids, and organic molecules[END_REF][START_REF] Ferré | Complete-active-space selfconsistent-field/Amber parameterization of the Lys296-retinal-Glu113 rhodopsin chromophore-counterion system[END_REF] The corresponding S0 equilibrium geometry is obtained by QM/MM optimization at CASSCF(12,12)/6-31G*/Amber level. During the optimization, the QM layer, all side-chains and waters with at least one atom within 4 Å from any atom of the retinal chromophore were kept flexible while all remaining atoms were kept frozen. All calculations were performed with the Molcas/Tinker package. [START_REF] Aquilante | Molcas 8: New capabilities for multiconfigurational quantum chemical calculations across the periodic table[END_REF][START_REF] Ponder | Tinker Molecular Modeling Package[END_REF] Population dynamics simulation. The QM/MM model constructed above was employed to simulate a room temperature Boltzmann-like (i.e. vibrationally incoherent) distribution described by 200 initial conditions. Briefly, a 20 ns dynamics simulation was performed using the software Gromacs 42 starting from the S0 equilibrium geometry at the MM level of theory at 298 K and 200 snapshots (geometries and velocities) were extracted along such dynamic simulation at 100 ps time interval. Each snapshot was propagated at HF/6-31G*/Amber level