Photochemistry and transient intermediates in a bacteriophytochrome photocycle revealed by multiscale simulations
Phytochromes are photoreceptors responsible for sensing light in plants, fungi and bacteria. Their photoactivation is initiated by the photoisomerization of an embedded chromophore, which triggers a large conformational change in the structure of the entire protein. Although phytochromes have been subject of numerous studies, the photoisomerization mechanism and the following reaction path leading to the final active state remain elusive. Here, we use an integrated computational approach that combines nonadiabatic surface hopping and adiabatic ground-state molecular dynamics simulations to gain atomistic details on the photoactivation mechanism of Deinococcus radiodurans bacteriophytochrome. Our simulations show that the ps-scale photoisomerization of the chromophore proceeds through a hula-twist mechanism that forces a counterclockwise rotation of the D-ring. The initial photoproduct rapidly evolves in an early intermediate which we characterize through IR spectroscopy simulation. The early intermediate then evolves on the ns-to-µs scale to a late intermediate, characterized by a more disordered binding pocket and a clear weakening of the aspartate-to-arginine salt bridge interaction, whose cleavage is essential to interconvert to the final active state.
Introduction<!>Methods<!>Surface Hopping<!>Ground-state MD simulations<!>Simulation of IR spectra<!>Results and discussion<!>The photochemical process<!>76<!>From the photoproduct to the early Lumi-R intermediate<!>45,71<!>From the early to the late Lumi-R intermediate
<p>Photoreceptors are the special intermediaries responsible for sensing light and activating the photoinduced biological function in many different organisms [1][2][3] . Phytochromes are an example of such complex biological machines, ubiquitous in plants, fungi and bacteria [4][5][6][7][8][9][10][11][12][13][14] . They are red-absorbing homodimeric proteins, carrying a photoswitchable bilin as a chromophore.</p><p>Most phytochromes exist in two distinct photoreversible forms: the red light-absorbing form, also known as Pr state, and the far-red light absorbing form, the Pfr state 12,[15][16][17][18] , which differ by both chromophore stereochemistry and protein structure. In canonical phytochromes the most thermodinamically stable Pr state is the resting (inactive) form. Red-light absorption triggers the conformational switch towards the Pfr state. The "reverse" reaction from Pfr to Pr proceeds via a thermal pathway or by absorption of a second photon of far-red light.</p><p>Phytochromes of different species share a generally conserved photosensory module (PSM, Fig. 1A), composed by three domains PAS (Per/Arndt/Sim), GAF (cGMP phosphodiesterase/adenyl cyclase/FhlA) and PHY (Phytochrome specific). The PSM is then linked to an output module, which exerts a specific function depending on the species [19][20][21] . Due to such a modular domain architecture, and because they absorb and fluoresce within the so-called "transparent window", phytochromes are exploited in emerging fields as bioimaging and optogenetics 13,[22][23][24][25][26] .</p><p>In bacterial phytochromes, such as the phytochrome from Deinococcus radiodurans 27-30   (DrBph), a biliverdin IXα (BV) chromophore, embedded in the GAF domain, is covalently bound to the protein through a cysteine residue (Cys24), belonging to the PAS domain.</p><p>Noteworthy, the "tongue" structural motif of the PHY domain connects to the chromophorebinding pocket through a salt-bridge involving a conserved aspartate belonging to the GAF domain (Fig. 1B). In response to red light, the chromophore is electronically excited and photoisomerizes at the C 15 =C 16 double bond 31 (Fig. 1A,C). Following the photoactivation, the salt bridge breaks, allowing the tongue to refold completely from a β-sheet to a α-helix, and finally leading to downstream changes in the entire phytochrome. The sequence of events that leads to such structural changes however is still not understood, due to the multiple time and space scales involved in the process. 16,20,21,28,29,32-40   Thanks to spectroscopic investigations, we know that the photocycle of bacteriophytochromes involves at least two intermediate states, namely Lumi-R and Meta-R 36,[41][42][43] . The former arises from the photoisomerization; then, on the microsecond time scale, thermal relaxation processes lead to the formation of the Meta-R intermediate. However, the structures of these intermediates have not been resolved, and it is still not clear how the chromophore structural change propagates first to the binding pocket, and then to the entire protein.</p><p>Already at the level of the photochemistry, there is no consensus on the direction of rotation around the C 15 =C 16 double bond. Comparison between circular dichroism (CD) spectra of the Pr and Pfr forms 44 of DrBph suggested a clockwise rotation of D-ring (Fig. 1C).</p><p>This hypothesis is supported by steric reasons, as the clash between the two methyl groups in C 13 and C 17 (Fig. 1C) would hinder a counterclockwise rotation. The crystal structure of the Pfr state 32 seems to support this hypothesis. On the other hand, a more recent study using time-resolved serial X-ray experiments 45  Through an integrated multiscale approach that combines non-adiabatic QM/MM excitedstate dynamics with ground-state molecular dynamics simulations, we investigate the photochemical pathway of DrBph, revealing the structure of the Lumi-R intermediate(s).</p><p>Our results show the presence of a conical intersection (CI) between S 0 and the first electronic singlet excited state of BV as a non-radiative decay channel. Such a CI is reached through a counterclockwise rotation of the D-ring with a hula-twist mechanism, resulting in a structure that is here interpreted as the early Lumi-R intermediate. Notably, this intermediate is characterized by the interaction of the chromophore D-ring with a tyrosine residue (Tyr263) in line with mutational studies 38,46 that demonstrated the tight connection between Tyr263 and the chromophore state. The IR spectrum calculated on the proposed structure shows a good agreement with the experimental data. Finally, the thermal relaxation into the late Lumi-R state is simulated, confirming that the latter is characterized by a more dynamic environment around the chromophore. In particular, we observe a weakening of the GAF-PHY Asp277• • • Arg466 salt bridge interaction, whose cleavage is essential to interconvert from Pr to Pfr.</p><!><p>We have used a strategy which relies on the integration of semiclassical surface hopping (SH)</p><p>and Born-Oppenheimer (BO) molecular dynamics (MD). In particular, we have used MD simulations based on both Molecular Mechanics (MM) and quantum mechanics/molecular mechanics (QM/MM) potentials.</p><!><p>The excited state dynamics was simulated by means of a semiclassical trajectory surface hopping 47 approach in a hybrid QM/MM scheme with electrostatic embedding. Due to large size of the protein (Fig. 1A) and the associated high computational cost, we simulated a model system, consisting of one monomer and a sphere of water molecules centered on the chromophore (Fig. 2B). The coordinates for this model were extracted from a 4 µs-long MM MD of the Pr state 48 . The resulting system was neutralized with sodium ions. In our QM/MM scheme, the BV chromophore was treated as QM subsystem, while the rest of the model system was treated with the AMBER ff14SB force field 49 (see Fig. 2B). To describe the electronic structure of BV, we used the multi-reference FOMO-CASCI 50,51 method in a semiempirical AM1 framework, with an active space of (6,6) and a Gaussian width of 0.05 a.u. for FOMOs. Details on the selection of the computational method are provided in the SI (Section S1.1).</p><p>The initial nuclear ensemble for the SH trajectories was obtained by propagating three ground-state QM/MM simulations (from here on, R1, R2 and R3) of the model system (Fig. 2B) for 20 ps, with a time step ∆t = 0.1 fs and the Bussi-Parrinello stochastic thermostat at 300K. The first 10 ps of each replica were discarded to avoid nonequilibrium effects due to the heating step. The last 10 ps were used for sampling the initial conditions (nuclear coordinates and momenta) for the swarm of non-adiabatic trajectories (from here on, Set1, Set2, Set3) 51,52 .</p><p>Both ground-state and surface hopping MDs used the same integration time step of 0.1 fs and were performed with a modified version of the semiempirical program MOPAC2002 53 , interfaced with the molecular dynamics TINKER 54 software package. In surface hopping calculations the first three adiabatic singlet states were considered, in order to account for the possible role of higher-lying states. The local diabatization algorithm 50,51 was used for the integration of the time-dependent Schrödinger equation for the electrons, and quantum decoherence was taken into account applying an energy-based correction 55 . More details are provided in the SI (Section S1.1.2).</p><!><p>Final configurations were extracted from the surface hopping trajectories to propagate adiabatic MD simulations in the ground state. Both QM/MM and MM descriptions were used For the QM/MM simulations, the QM subsystem (the same as for surface hopping simulations) was treated using the AM1 semiempirical Hamiltonian, while the rest of the system (monomer and water box) was treated at the MM level (with the AMBER ff14SB forcefield for the protein, while water was described with the TIP3P model). Four trajectories were propagated for 5 ns each and average-quantity analysis were made on 5000 frames per trajectory.</p><p>For MM MDs, the GAFF force-field was used for BV and for the bonding interactions between BV and the covalently bonded cysteine residue Cys24. Ten independent replicas were run for 1 µs each. All MD simulations were run with AMBER18 56,57 . More details are provided in the SI (Section S1.2).</p><!><p>In order to perform the IR spectroscopy simulation of the Lumi-R intermediate, we have extracted 57 frames from the BO QM/MM MDs. From these frames, we extracted the protein and a shell of 390 water molecules centered on the chromophore. On the selected subsystem, we have performed a ONIOM(QM:MM) 58  The complete pipeline for the calculation of the frequencies is reported in the SI (Section S3). All geometry optimizations were performed using the Gaussian 16 code 59 .</p><!><p>In the following, we first present and discuss results of the non-adiabatic MD simulation of the photoisomerization process. Then, we present the results of the ground state MD trajectories, which started from the structures reached through the non-adiabatic dynamics.</p><p>The structure of the intermediate identified in this analysis is then characterized in terms of its IR spectrum, which is compared with that of the Pr state. The resulting difference spectrum is compared with the available experimental data. Finally, the further evolution of the intermediate in the µs time scale is explored through µs MM-MD trajectories.</p><!><p>According to the literature and to our previous MD simulations 36,48,[60][61][62][63][64][65][66][67] , the Pr state is conformationally heterogeneous. In order to capture such heterogeneity, the initial nuclear ensemble of the system was prepared via three ground-state QM/MM MD simulations (R1-R3 in the following) performed on the model system (Fig. 2B), as described in Section 2. The heterogeneity of the Pr conformation is reflected in our QM/MM MD replicas. Although many hydrogen-bond interactions are conserved in all the replicas, the residues close to the Dring of BV show a different behaviour (Fig. 2A). In particular, we observe a closer interaction between the D-ring carbonyl (O D ) and His290 in trajectory R2. In the other trajectories this interaction is either absent (R1) or unlikely (R3). In addition, two methionine residues (Met174 and Met264) are closer to the D-ring in R2 and R3 (Fig. 2C,D), which suggests a more compact structure of the chromophore binding domain in these replicas.</p><p>The different local environments are also reflected in the distribution of the dihedrals D5 and D6 of BV (Fig. S6, Fig. 1C), which are expected to be involved in the isomerization. A shift in the equilibrium values of D6 and D5 are observed in R2 with respect to R1 and R3, which can be related to the stronger hydrogen bond with the His290 residue.</p><p>Configurations and momenta extracted from the ground-state QM/MM MDs were used as initial conditions for surface hopping trajectories, which were all initialized in the lowest S 1 excited state. We analyzed 3648 non adiabatic trajectories in total and we classified them as Set1, Set2 and Set3 on the basis of the starting R1, R2, R3 trajectories, respectively.</p><p>In all simulations, the excited-state dynamics proceeds from the Franck-Condon point through a counterclockwise rotation around the double bond C 15 =C 16 , i.e. by increasing D6 (Fig. 3A). The simultaneous clockwise rotation of D5 reveals a concerted isomerization around the double and single bonds (for the movies of the whole dynamics, see Supplementary  Population (%) Movies 1 and 2). Such an opposite rotation of two adjacent dihedrals corresponds to a hula-twist mechanism, a space-saving isomerization inside a tight binding pocket 68,69 . This mechanism has been previously proposed for other conjugated systems, including photoactive protein chromophores 68,70 .</p><p>In all trajectories, the nonradiative decay from S 1 to S 0 proceeds through a conical intersection (CI), reached when D6 is around 90 degrees (Fig. 3A). At the same time, D5 is almost planar. Representative structures of the Franck-Condon point and CI are compared in Fig. 3C, D. After reaching the ground state through the CI, most trajectories return to the initial conformation, i.e. the Pr state (from here on, non-reactive trajectories), while a small population continues the rotation towards larger values of D6 (from here on, reactive trajectories). The final geometry is represented in Fig. 3E.</p><p>As shown in Fig. 3B, the opposite motion of dihedrals D6 and D5 causes the D ring to rotate counterclockwise. Namely, the angle between the planes of rings C and D increases by ∼55 degrees on average. None of the 3648 trajectories here analyzed (see Fig. 3A and Fig. S8) showed clockwise rotation as a possible reaction path. On the contrary, both reactive and non-reactive trajectories showed a similar initial evolution of dihedrals D5 and D6 (Fig. 3A, Fig. S8, Supplementary Movie 2). These findings are in disagreement with a previous study based on CD experiments 44 , which proposed clockwise rotation as the most probable mechanism. However, the clockwise mechanism was suggested on the comparison of the initial and final states (Pr and Pfr) and on steric and chromophore-confined arguments, without taking into account the conformational dynamics of the system or the role of the protein environment. Conversely, a counterclockwise rotation was directly observed in timeresolved X-ray experiments 45,71 , which are consistent with the time scales simulated here.</p><p>In addition to the rotation of D5 and D6 dihedrals, also the propionyl group linked to the B ring changes its conformation, reaching the one observed in the final Pfr state (Fig. S9, Fig. S10). However, it may be noted that in our simulations the D ring does not reach the representative conformation of the Pfr state. In fact, the clash between the two methyl groups of rings C and D (Fig. 3E, Fig. 1C) prevents from reaching this conformation.</p><p>In many trajectories, the hydrogen bond between the D-ring carbonyl and His290 is lost upon reaching the CI region, while the interaction with a water molecule remains unchanged (Fig. 3D and Fig. S7). In particular, for the trajectories generated from the R2 replica (Set2) the probability of the O D • • • His290 hydrogen bond changes from 31% to 10%, while in Set3 it changes from 7% to 4%. On the contrary, none of the trajectories generated from R1 (Set1) show this interaction. Interestingly, when BV is back in the ground-state, we observe the formation of a hydrogen bond between the D-ring amide hydrogen and Tyr263 (Fig. 3E, Fig. S7). The time evolution of the electronic excited state populations is shown in Fig. 3F. All three sets of trajectories show exponential decays with a time constant τ in the order of 0.4-2.3 ps, with Set2 being characterized by a significantly slower decay than Set1 and Set3 (see Tab. S3). According to the literature, the protein environment has a critical role in defining the excited state lifetime of phytochromes 67,[72][73][74][75] . Our results are in line with this picture.</p><p>In fact, comparing the results obtained starting from R2 with those of the other replicas, it is clear that a more compact structure of the protein pocket hinders the isomerization, and induces a slower dynamics.</p><p>The three sets of surface hopping trajectories also agree with the low experimental photochemical yield of the Pr to Lumi-R transition, which is around 10-20% 60,67 (Table S3).</p><p>The experimental low quantum yields of fluorescence and photoproduct formation (Lumi-R) suggest that the largest decay channels for phytochromes are the non-productive ones</p><!><p>(neither fluorescent nor Lumi-R-forming channel). Our results are in agreement with this picture: the excited-state of BV decays via a conical intersection, but only 10-20% of the trajectories are reactive. The trend observed for the different quantum yields of the three sets of SH trajectories (see Tab. S3) is reflected in their decay rates. In fact, Set1 has a shorter decay time than Set3, which in turn decays faster than Set2 (Fig. 3F). The same trend can be seen in both reactive and non-reactive trajectories (Tab. S3). Therefore, the Pr heterogeneity is reflected on both quantum yield and decay rates.</p><p>There are two factors that concur in determining the photoisomerization quantum yield:</p><p>the "early decay" and the angular inertia. Early decay means that the S 1 →S 0 hop occurs when the dihedral D6 has not yet reached ∼90 • , so that the molecule is still on the reactant side: once on the ground state, it returns to the starting configuration. In reactive trajectories, the decay occurs closer to the CI than in the nonreactive ones (Tab. S3). The angular inertia refers to the angular velocity of the reaction coordinate at the time of the hop. We see that reactive trajectories have a larger average velocity along D6 than the nonreactive counterparts. In other words, when reaching the conical intersection, only the fastest com-ponent proceeds towards the product. The trend observed for the decay rates of the three sets of SH trajectories correlates with both the average D6 value and its time derivative at S 1 →S 0 hops (Tab. S3). This in turn determines a lower quantum yield for the set that decays more slowly.</p><p>Our non-adiabatic dynamics simulations show that the dihedral D6 degree of freedom is not sufficient to describe the photoisomerization of the chromophore. In fact, the latter occurs through a hula-twist mechanism, which involves simultaneous D5 rotation forcing the D ring to undergo a counterclockwise motion (Fig. 3B). This process is accompanied by the loss of hydrogen bond between the N D H group and a water molecule, and between the O D carbonyl and His290. A new interaction is formed between the N D H group and Tyr263. This primary photoproduct is distinct from the final Pfr state, and features a new conformation around the C 14 -C 15 single bond. A full rotation of D ring is not yet possible, due to the steric restraints imposed by the protein pocket.</p><p>Finally, we note that our estimated timescales for the photochemical process are significantly shorter than what seen in transient absorption measurements, which for the PSM of DrBph phytochrome report an excited-state lifetime of about 170 ps 63 , attributed to a concomitant reorganization of the hydrogen-bonding network in the chromophore binding pocket. On the other hand, in support of our results, the more recent femtosecond X-ray crystallography study 45 has detected a twist of the D-ring already after a delay time of 1 ps after photoexcitation. To conclude, there is not yet consensus about the excited state lifetime of phytochromes, as it has been shown to strongly depend on the type of construct (i.e., chromophore binding domain, photosensory module or full-length system) 63,77 . Moreover, a more complex excited state decay with multi-exponential kinetics has been seen in some bacterial phytochromes 60,63,64,73,78 .</p><!><p>Using as starting point the final ground state configurations of the SH reactive trajectories, we have performed adiabatic simulations to characterize the photoproduct and its time evolution. To perform this analysis we considered both QM/MM ground-state MDs and the first 10 ns of the MM MD simulations. The selected time windows was chosen taking into account the experimental data suggesting the ns scale for the formation of the early intermediate.</p><!><p>Both sets of simulations indicate a stable structure of the chromophore, with dihedrals D5 and D6 remaining close to the values obtained in the photoisomerization (Fig. 4B,C). Only for D6 the MM simulation describes a more planar structure. All the other dihedrals remain nearly the same as in the Pr state (Fig. S10), similarly to what observed in the analysis of SH simulations. The only exception is represented by the dihedral of the propionyl group (D1B), which is stabilized at Pfr-like values in the QM/MM MDs (Fig. S10).</p><p>To better characterize the structures reached by the QM/MM and MM MDs, we have calculated the root-mean-square displacement (RMSD) of the chromophore binding pocket using as a reference the Pr crystal structure (Fig. S11A). This analysis shows that both QM/MM and MM MDs are stable around an equilibrium value of about 1.6 Å. In order to compare the configurational space explored by the QM/MM and MM MDs, we have focused on the interaction between the chromophore and the closest residues (Fig. 4 and Fig. S12).</p><p>The two sets of simulations seem to explore roughly the same space, but with different distributions. For example, a water molecule interacts more tightly with the O D atom in the MM-based simulations. On the other hand, both QM/MM and MM MDs show a stable hydrogen bond between the D ring nitrogen (N D ) and Tyr263 (Fig. 4F). This interaction, already present in the photoproduct, is likely to stabilize this early intermediate.</p><p>To have a more quantitative comparison of the configurational space spanned by QM/MM and MM MDs, we have used a principal component analysis (PCA) based on intermolecular distances involving key-residues in the chromophore binding pocket (Fig. S13). This analysis reveals that the QM/MM MDs explore a subset of the MM MD configurations (Fig. 4G). We also analyzed the configurational space through a hierarchical clustering algorithm. Details on the clustering procedure and on the selection of the clusters are provided in the SI (Section S2). Three clusters (named 0, 1 and 2) were identified and projected in the PCA space (Fig. 4H). Cluster 2 differentiates from the others due to a large distance between Tyr263 and Asp207 (Fig. S14). Cluster 1, on the other hand, is mainly characterized by a large distance of the D-ring carbonyl oxygen to the nearest water molecule (Fig. S14). Both QM/MM and MM MDs are characterized by an evolution of the cluster population within the first 5 ns (Fig. S15). Initially the system is mainly described by cluster 1. For QM/MM MD the description remains unchanged up to 4 ns, when the system evolves from cluster 1 to cluster 0. The MM MD, on the other hand, evolves much earlier, and after 10 ns populates only clusters 0 and 2.</p><p>The structural ensemble reached after photoisomerization is characterized by a stable BV and dynamic, but still recognizable, interactions within the chromophore-binding pocket. We assign this conformation to the early Lumi-R state. Mutagenesis studies 38,46 have suggested that the Y263F mutation hinders the formation of Lumi-R, thereby reducing the photoconversion (Pr→Pfr) yield: our simulations show that Tyr263 is indeed fundamental for stabilizing the Lumi-R state through a hydrogen bond with the N D atom of BV.</p><p>With the exception of the interaction with Tyr263, the hydrogen bonding network characterizing the Pr state 61 is conserved in our simulations of the intermediate. This finding suggests that a more complex rearrangement of the chromophore-binding pocket has to occur for the system to reach the Pfr state.</p><p>In order to further validate our proposed early Lumi-R structure, we have computed the IR spectrum as described in the Methods section.</p><p>Experimentally 37 , the difference spectrum between the Lumi-R and Pr states was measured by step-scan IR spectroscopy: also here we report the difference with respect to the Pr state. Moreover, here we compare our results to the spectrum recorded in D 2 O, as the measurements in H 2 O also feature signals arising from the protein modes, which we do not compute here.</p><p>Since we are interested in the region of carbonyl stretches, the spectra reported in Fig. 5 use only the CO D and CO A normal modes (for more details see Section 2.2.1). The Lumi-R spectra have been obtained for the three different clusters identified in the previous section.</p><p>In particular, we have used 18, 26, and 13 configurations belonging to clusters 0, 1, and 2, respectively.</p><p>In Fig. 5A we can see how cluster 0 (Cl0) and cluster 1 (Cl1) give a rather similar representation, except for a shoulder in Cl0 at frequencies around 1700 cm −1 . Cluster 0 is characterized both by configurations in which the CO D is at hydrogen bonding distance with a water molecule and configurations in which such interaction is absent (Fig. S16). Cluster 2 (Cl2), in which the Tyr263• • • Asp207 hydrogen bond is absent, shows a major CO D peak at frequencies very close to the shoulder of Cl0. In fact, this signal is due to configurations in which a water molecule is hydrogen bonded to the carbonyl. On the contrary, configurations in which such interaction is absent show a peak at frequencies comparable to Cl1.  (B) In the upper panel, the averaged spectrum of the Lumi-R was represented together with that of the Pr state. In the lower panel, the theoretical (calc.) and experimental 37 (exp.) Lumi-R−Pr difference spectra were represented. Since, experimentally, a shift of 8 cm −1 was observed as a result of deuteration, our reference spectrum has been shifted by that amount.</p><p>An average spectrum was obtained through a weighted average of the three clusters (Fig. 5B). We have used the weights which reflect the population of the three clusters in the QM/MM MDs (0.1825, 0.8085 and 0.0090 for Cl0, Cl1 and Cl2, respectively). Finally, we computed the difference spectrum (Lumi-R − Pr) and compared it with experiments (Fig. 5B, bottom panel). The calculated difference spectrum reproduces quite well both and shape of the main positive and negative CO D signals. In particular, the correct reproduction of the frequency shift experimentally observed when moving from Pr to the Lumi-R, further increases our confidence on the validity of the obtained structure of this intermediate. It has in fact to be stressed that the difference spectrum is determined by a delicate interplay between the structural change in the chromophore and the modified chromophore-residues interactions. As shown in Fig. S17, the internal degrees of freedom of BV do not show noteworthy differences compared to what observed for the early Lumi-R structure. On the other hand, significant differences are found in the binding pocket, as shown by the larger RMSD fluctuations. (Fig. S11B). Going into detail, some trajectories evolve towards a configuration in which Arg466 has an enhanced mobility and, at times, it is very close to the carbonyl D-ring (Fig. 6B). In fact, analyzing the distribution of distances between BV and the closest residues (Fig. S18), it is clear that some of these distances present multimodal distributions. For example, the distance between the D-ring carbonyl and Arg466 has a bimodal distribution, indicating a dynamic behaviour of Arg466 (Fig. 6B), while in the early Lumi-R intermediate the narrow distribution indicates a rigid configuration. The dynamic behaviour of this residue has been also confirmed by a recent study 37 . This is of considerable interest, as in the Pr state Asp207 forms a strong salt bridge with Arg466, whereas this contact is replaced in Pfr by a hydrogen bond with Ser468.</p><!><p>This increased mobility for Arg466 has a knock-on effect on the residues around the Dring: the salt-bridge interaction between Arg466 belonging to the PHY domain and Asp207 between Asp207 and Arg466. In particular, we define a contact if the oxygen (for Asp207), and nitrogen (for Arg466) are at a distance lower than 3 Å. All the analysis were performed on 10000 frames.</p><p>belonging to the GAF domain, a hallmark feature of the Pr state, becomes weaker (Fig. 6G), as well as the interaction between Asp207 and Tyr263 (Fig. 6F). A water molecule shields the interaction between Asp207 and Tyr263, causing its cleavage and the opening of the "cage", formed by a hydrogen bonding network among Arg466, Asp207 and Tyr263, with the final exit of the pyrrole water (Fig. 6F). As a result of the weakening of the interaction between Asp207 and Tyr263, the hydrogen bonding interaction between the D-ring and Tyr263 becomes weaker as well (Fig. 6E), and the aforementioned water molecule takes the place of Tyr263 in the interaction with the D-ring (Fig. 6F). Therefore, according to our simulations, the pyrrole water is not ejected during the photochemical event, but when the Asp207• • • Tyr263 interaction is broken. Moreover, the hydrogen bond between B-ring propionate and Ser257 becomes weaker (Fig. 6D), while the one with Tyr216 is strengthened.</p><p>In conclusion, our simulations show two sequential reaction intermediates: early and late Lumi-R. Both structures are stabilized through an interaction with Tyr263. The steric clash between methyl groups in C-and D-rings has not been overcome yet. However, it is worth observing the bond angle between the C-and D-rings (here defined by the atoms C 1 4-C 1 5-C 1 6, see Fig. 1C): the distribution of this angle is centered at larger values than in the Pr or Pfr states (Fig. S10). In fact, such coordinate is expected to play an important role in overcoming the steric clash between the methyl groups.</p><p>The largest structural changes between early and late Lumi-R are represented by intermolecular features, mainly involving the Arg466 residue, which couples the chromophore to the tongue, in accordance with the literature. The salt bridge Arg466• • • Asp207 has a central role in unleashing the activation pathway of the phytochrome since this breaking would free the tongue for refolding. Thus, at a longer time scale, we can expect the detachment of Arg466, allowing the tongue refolding and the formation of the Pfr state.</p>