Phytochromes are widely found in plants and microorganisms, which are responsible for fundamental photo-response functions [1-4]. They play crucial roles in a variety of light-regulated biological processes, ranging from seed germination, growth and flowering in plants to pigmentation and phototaxis in bacteria [1, 5, 6]. The phytochromes can be engineered into optical probes for the deep-tissue imaging of living cells, thus their photoinduced properties are considerable significance in photobiological researches [7-9].

The chromophores in phytochromes, as the photoreaction centers, are mainly classified as three categories, namely biliverdin (BV), phytochromobilin (PΦB) and phycocyanobilin (PCB). Their structures are very similar and all share a methane-bridged open-chain conjugated tetrapyrrole structures with differences on the two terminal rings (Scheme 1). Their basic photochemical reaction mechanism is the reversible photo- and thermal-conversions between two isomers, i.e., the red-absorbing form Pr and the far-red-absorbing form Pfr. The primary step of the Pr-to-Pfr photoisomerization mechanism is governed by the nonadiabatic transitions. Once red-light absorption at ~660 nm, the chromophore in Pr state undergoes an ultrafast Z-to-E photoisomerization accompanied by the internal rotation of pyrrole rings, then follows conformational relaxations and finally leads to the formation of the Pfr isomer via the thermal processes. The reverse Pfr-to-Pr conversion experiences the E-to-Z isomerization upon the photoexcitation with far red light at ~730 nm. The forward and backward isomerization reactions induce the large-sized conformational evolution of the whole photoreceptor protein, giving regulating signals to adjust relevant biological functions. Therefore, their photoinduced reactions have been extensively investigated experimentally and theoretically [10-27].

Although numerous researches were conducted for many years, there are still several fundamental open questions in the photoinduced Pr-to-Pfr processes. The first critical argument is which molecular motion dominates the Z-to-E isomerization in the primary excited-state reaction. Several studies suggested that in protein environments the photoisomization takes place around the C₁₅C₁₆ methine bridge, resulting in the torsion of nearly-free side Ring D along this bond [11, 13, 28, 29]. However, different isomerization mechanisms were also proposed. For instance, Ulijasz et al. discovered that the internal rotation of the Ring A via the twisting motion of the C₄C₅ double bond triggers the excited-state dynamics of the PCB chromophore in photoreceptor protein [21]. It is also important to clarify the microscopic photoisomerization mechanism of the isolated chromophore at the all-atomic level. Some studies pointed out that the initial excited-state evolution of the isolated chromophore may be governed by the torsional motion at the C₁₀C₁₁ bond [20, 30, 31]. Zhuang et al. confirmed the pivotal role of the C₁₀C₁₁ double bond, and proposed that C₉C₁₀ single bond is also an key rotational degree of freedom in the photoisomerization of the isolated PΦB chromophore [20]. The existence of different photoisomerization mechanisms for the isolated chromophore and chromophore in phytochrome protein pocket is possibly due to the geometrical constrains added by the interactions between the chromophore and the surrounding protein. In this sense, the deep understanding of the role of the protein pocket in the modification of the photoinduced dynamics of phytochrome chromophore still requires additional efforts. Furthermore, hydrogen bonding plays a crucial role in many biological processes and accounts for many important physical properties of compounds [32-36]. Several works pointed out that the photochemistry of phytochrome chromophore is also strongly dependent on the hydrogen-bond network, as well as the protonation status of the pyrrole rings and surrounding residues [16, 22, 37-41]. A vibrational spectroscopic study on bathy phytochromes demonstrated that the proton transfer between the tetrapyrrole and a nearby histidine residue directly modules the photobiological functions by the generations of the photoactivated and deactivated states [22]. Osoegawa et al. revealed that the deprotonation of the Ring B/C distorts the methane bridge between the Ring B and the Ring C of the PCB chromophore, and such structural difference further adjusts the green or red light absorption [40]. Modi et al. clarified that the protonation/deprotonation status may regulate the light aborption of the BV chromophore in the protein pocket [24]. Kennis and co-workers proposed that the excited-state hydrogen transfer from the Ring A, B or C directly leads to the nonadiabatic decay to the ground state in the BV chromophore, resulting in the decreasing of the fluorescence signals [37]. Recently, Liu et al. demonstrated that the different orientations of the aspartate residue give the distinctive hydrogen-bond network connectivity that regulates the excited-state reaction channels [16]. As the summary, the protonation/deprotonation of the hydrogens attached to the pyrrole rings plays a critical role in the nonadiabatic dynamics of the phytochrome chromophore. It is interesting to clarify whether the removing of the hydrogen atom connected with the pyrrole nitrogen significantly modifies the nonadiabatic dynamics of the phytochrome or not. Particularly, it is also necessary to identify which key coordinates are responsible for the photoinduced dynamics of different phytochrome derivatives obtained by simple chemical modifications. These questions can be clearly addressed by the simulation of the nonadiabatic dynamics at all-atomic level, in which the active degrees of freedom in the excited-state dynamics are clarified by tracking the real-time molecular motion. In recent years, the on-the-fly ab initio nonadiabatic dynamics simulations were employed extensively to explore the ultrafast nonradiative reactions [42-47]. These simulations directly give the rich information, such as the time-dependent population dynamics, the nonadiabatic decay pathways and the essential molecular structure evolution.

In present work, we wish to provide the detailed description for the photoisomerization of the isolated BV chromophore by using the on-the-fly nonadiabatic dynamics. The simulation was run with the trajectory surface hopping (TSH) approach [48] and the electronic structure calculations were performed at the semi-empirical OM2/MRCI level (the orthogonalization-corrected OM2 Hamiltonian combined with multireference configuration interaction) [49, 50]. Many studied have verified that the nonadiabatic dynamics at the OM2/MRCI level is a powerful simulation tool in the description of excited-state processes for media-sized and large polyatomic systems with the balance of computational accuracy and efficiency [51-56]. By following the excited-state molecular motions, we found that the protonation/deprotonation states of the pyrrole rings show the strong impact on the excited-state lifetime and reaction pathways. The underline reason is that the hydrogen bond network deeply modules the intramolecular motion, leading to different profiles of excited-state potential energy surfaces. As the consequence, the distinct population dynamics appears and different conical intersections are involved. Such finding confirms the important role of the hydrogen bond in the photoinduced nonadiabatic dynamics of the phytochrome chromophore. In other words, the understanding of the photoisomerization mechanism of phytochrome chromophere must take the hydrogen bond network into account, beside the steric effects caused by its surrounding protein environments. This work provides new insights into the photobiological functions of phytochrome chromophere and suggests the possible ideas to control its photoconversion processes for the further rational engineering of near-infrared fluorescent proteins.

We focus on the nonadiabatic dynamics of the Pr conformation of the BV chromophore in DrBphP (Deinococcus radiodurans) bacterial phytochrome. The ZsZsZa (C₅-Z, syn C₁₀-Z, syn C₁₅-Z, anti) stereoisomer (Fig. S1 in Supporting information) is chosen as our target configuration, which was identified as the dominant isomer by crystallography and NMR spectroscopy [57-59]. We built the simplified BV model to perform the nonadiabatic dynamics simulation, in which some side groups (-CH₂COO⁻ and -Cys-S) were replaced by hydrogen atoms to reduce the computational cost. As it is rather difficult to assign the protonation status by experimental studies, we consider all possible single deprotonation derivatives. In order to systematically rationalize the influence of phytochrome chromophore protonation state on the photoinduced processes, the four most probable protonation models (Fig. 1 and Fig. S2 in Supporting information) were investigated, namely:

(1) BV_N: all nitrogen atoms protonated;

(2) BV_C: nitrogen in the Ring C deprotonated and the Rings A, B and D protonated;

(3) BV_B: nitrogen in the Ring B deprotonated and the Rings A, C and D protonated;

(4) BV_A: nitrogen in the Ring A deprotonated and the Rings B, C and D protonated.

We first optimized the structures of the S₀ and S₁ minima (S₀_min and S₁_min) for the four protonated/deprotonated BV chromophore models (BV_N, BV_C, BV_B and BV_A), and the key internal-coordinate parameters of these geometries are listed in Tables S1-S4 (Supporting information). The most striking feature here is that the different protonation patterns clearly build distinguishing hydrogen-bond networks, as shown in Fig. 1, Figs. S3 and S4 (Supporting information). Here, the interaction region indicator (IRI) method was employed to investigate the intramolecular hydrogen bond (Fig. S3) [60]. For instance, at the S₀_min, only a weak hydrogen bond between O_A of the Ring A and H_CD atom (bridged the Ring C and the Ring D) is found in BV_N. However, the more complex hydrogen-bonding interactions are observed due to the presence of the multi-hydrogen-bond connectivity in the other three models (BV_C, BV_B and BV_A), as the results of the deprotonation of different pyrrole ring.

Fig. 2 presents the time-dependent fractional occupations of the S₀ and S₁ electronic states for the above four models, assuming that all trajectories start from the S₁ state. Their lifetimes are far from each other and their corresponding time scales are even differed by the order of magnitude. The BV_N model with the saturated protonation displays the ultrafast excited-state population decay. The S₁ population decay to ~50% around 400 fs and the S₁ population becomes very low at about 1000 fs. As the contrast, the excited-state population decay of the BV_C model takes place at a much longer time scale and less than 50% of S₁ population decays to the S₁ state within 2.0 ps. The S₁ lifetimes of the latter two models (BV_B and BV_A) are extremely long and only the very minor population decay is observed within the simulation time duration. The remarkably distinctive population decay behaviours among these four BV models imply that their excited-state reaction mechanism may be drastically different. In this scenario, two open questions naturally arise: (1) How to explain the different population dynamics in these BV models and do they display different decay channels? (2) How to correlate their decay mechanism and the protonation status? Next, we try to clarify both questions.

We begin with the BV_N model to explore the mechanism of the nonadiabatic dynamics. By analysing all S₀/S₁ hopping geometries and the CI optimization results, we found that two S₀/S₁ minimum-energy conical intersections (CI₁ and CI₂ in Figs. 3a and b) govern the nonadiabatic transitions. The former CI₁ is the main decay channel, since 93.6% trajectories decay via it, as shown in Fig. 3c. In the beginning of the photoinduced process, both HOMO and LUMO are delocalized on two central rings (Rings B and C) at the Frank-Condon region. When the torsional motion drives the system to access the CI₁ and CI₂, the HOMO and LUMO tend to become the localized one. As shown in Fig. S5 (Supporting information), when such torsional motion take place at the C₁₀C₁₁ or C₉C₁₀ bond, HOMO and LUMO start to locate in two moieties separated by these two bonds. This type of biradical or charge-transfer character widely exists in photoisomerization of double-bond conjugated systems [61].

The first conical intersection CI₁ is characterized by the strong twisting motion of the C₁₀C₁₁ bond, along with the minor torsion of the C₁₅C₁₆ bond. From S₀_min to CI₁, the dihedral angle of C₉C₁₀C₁₁N_C between the Ring B and the Ring C changes from −12.4° to −85.4°, and the dihedral angle of C₁₄C₁₅C₁₆N_D between the Ring C and the Ring D also changes slightly from 6.4° to 19.0°. For the trajectories via this channel, the rotation of whole C-D double-rings moiety takes place in the nonadiabatic dynamics, while the minor torsion of the terminal Ring D with respect to Ring C also exists.

The second conical intersection CI₂ is characterized by the significant twisting of the C₉C₁₀ bond and also the minor torsion of the C₅C₆ bond. From S₀_min to CI₂, the dihedral angle of N_BC₉C₁₀C₁₁ between the Ring B and the Ring C increases from −16.5° to −93.4°, while the dihedral angle of C₄C₅C₆N_B between the Ring A and the Ring B shows the minor decreasing from 28.2° to 17.6°. In this decay channel, the obvious rotation of the whole A-B double-ring moiety and the minor torsion of the terminal Ring A exist.

The reaction pathways from the S₀_min to two conical intersections were built by employing the linear interpolated method. Fig. 3d predicts that the CI₁ pathway should be more preferable, although both pathways are barrierless. Moreover, the S₁_min optimization directly gives the CI₁ geometry. These findings clearly explain why the CI₁ channel is dominant in the nonadiabatic dynamics of the BV_N model (Fig. 3c).

After the internal conversion of BV_N, the excessive kinetic energy leads to the high excitation of vibrational motions, and the system can easily overcome barriers on the electronic ground state to form various photoproducts. To explore possible reactions products, we ran the additional 5-ps Born-Oppenheimer molecular dynamics on the electronic ground state after the 2-ps nonadiabatic dynamics. By connecting the nonadiabatic and adiabatic dynamics, a few of representative trajectories are given in Fig. 4. Starting from the S₁ state, the trajectories on the excited state move towards two S₀/S₁ CIs. Since the CI₁ channel is dominant, we mainly show the geometry evolution of some representative trajectories decaying via it. At first, the trajectories propagate on the S₁ state. Once approaching the CI₁ where two electronic states show similar potential energies, the trajectories decay to the ground state. After the internal conversion, the initial BV_N isomer can be reformed though the reversed torsion of the whole C-D double-rings moiety, returning back to the starting BV_N configuration. Next, in some trajectories the vibrational energy flow may furtherly induce the internal rotation of the Ring D along the C₁₅C₁₆ bond, or the mixed Ring A and the Ring D rotations (Figs. S6-S8 in Supporting information). At the end of simulation, the ratio of these three channels is about 14.4%: 32.1%: 19.3%. Other isomers were also found while their probabilities are small (Fig. S9 in Supporting information). Since only 5.88% trajectories passing the CI₂, it is not necessary to discuss their final photoproducts. The current finding clearly confirms that the rotation of other rings (such as Ring D) may be highly possible even when such motion is not dominant in the excited-state dynamics. Nevertheless, such multiple-step reaction mechanism provides a possible explanation on why different isomerization patterns were discovered in previous studies [11, 13, 20, 21, 28-31].

Next, let us consider the deprotonated BV_C model, which shows the longer excited state lifetime than the BV_N model (Fig. 2). Different to BV_N, no significant torsional motion at the single and double carbon and carbon bonds between the Ring B and the Ring C was observed for BV_C in the trajectories. This difference can be understood by the presence of the strong hydrogen bond between deprotonated N_C and atom H_B in BV_C, which locks these torsional motions. As the consequence, the preferred isomerization pathways in BV_N are not observed in the excited-state dynamics of BV_C.

To clarify which motion is responsible for the nonadiabatic decay of the BV_C model, we collected the hopping geometries and located the key internal coordinates. As shown in Figs. 5a and b, two dihedral angles N_CC₁₄C₁₅C₁₆ and H_BN_BC₆C₉ are found to play the significant roles. Starting from these hopping geometries, the optimization of the minimum-energy S₀/S₁ CI geometry also confirms the importance of these two internal coordinates. Therefore, the synergy of the rotation of the Ring D and the out-of-plane motion of the H_B atom associated with the Ring B (H_BOOP) opens a new pathway in the nonadiabatic dynamics in the deprotonated BV_C model, different to the BV_N model.

The excited state reaction pathway of the BV_C model was given in Fig. 5c. Starting from Frank-Condon region, there is a barrierless pathway to the S₁_min. According to the linear-interpolated pathway from the S₁_min to the minimum-energy geometry in the S₀/S₁ CI (Fig. 5d), we noticed that the CI lies much higher than the S₁ minimum. Also considering the large size of the current chromophore, the vibrational relaxation on the excited state may lead to that the system tends to stay at the S₁ minimum region for very long time before it enters the CI region. This may explain the much slower nonadiabatic decay of the BV_C model with respect to the BV_N model. At the same time, this CI is accessible from the S₁_min only after the rather large geometrical adjustment on the excited state. Several motions must be involved, which include the significant H_BOOP motion (∠H_BN_BC₆C₉ is 125.7° at CI, 159.9° at S₁_min and 176.9° at S₀_min), the breakdown of the strong hydrogen bond N_CH_B and the torsion of the N_CC₁₄C₁₅C₁₆ dihedral angle by ~10°. As the consequence, it takes a longer time for trajectories to access this CI region, resulting in the slow population decay of the BV_C model.

Now we pay attention to the deprotonated BV_B and BV_A models, which show significantly longer excited state lifetime (Fig. 2). As the excited-state reaction pathways were clarified for BV_N and BV_C systems, it is necessary to examine whether these channels are possible in the BV_B and BV_A models.

In BV_B model, the hydrogen-bond network is observed clearly (Fig. S3). Similar to the BV_C case, the N_B-H_C hydrogen bond limits the torsional motion between the Ring B and the Ring C, preventing the ultrafast decay channel observed in the BV_N model. Meanwhile, the rotational motion of the Ring D is also blocked by the hydrogen bond between atom O_A and atom H_C-D. Moreover, as the nitrogen atom in the Ring B is deprotonated, there is no H_BOOP movement in this system. Therefore, BV_B cannot experience the similar decay channel as BV_C.

In BV_A model, the hydrogen bond between the Ring A and the Ring B freezes the rotation of C₄C₅C₆N_B. However, the twisting motion of the C₅C₆ bond is clearly involved in the early stage of the nonadiabatic dynamics for the BV_N model (Fig. S10 in Supporting information), although the torsional motions at the bonds between the Ring B and the Ring C are dominant. As such torsional motion is fully locked in the BV_A model (Fig. S11 in Supporting information), the corresponding channel is closed. In addition, the same hydrogen bond also restricts the H_BOOP motion of the Ring B, and the presence of the other strong hydrogen bond (S₁_min, O_AH_C-D = 1.629 Å in Fig. S3) prevents the torsional motion of the Ring D. Therefore, the decay pathway in the BV_C model cannot exist for the BV_A model.

Finally, we give some comments on the role of the hydrogen transfer in the nonadiabatic dynamics of the deprotonated BV_C, BV_B and BV_A models. Although the proton transfer may take place due to the existence of the hydrogen-bond network in the above three model (Table S5 in Supporting information), this motion is not important in the photoisomerization dynamics. Let us take the BV_C model with the highest nonadiabatic transition probability as an example, most hopping structures still remain as the BV_C configurations, and only 6.1% of them become the BV_B configurations. Even in this case, the existence of the hydrogen bond still locks the torsional motion at the single and double CC bonds between the Ring B and the Ring C. For the BV_A and BV_B models, the hydrogen transfer is also not important for the nonadiabatic decay dynamics and thus we do not repeat the discussions here. Overall, the proton transfer has very minor effects on the non-adiabatic transitions processes for the deprotonated models, possibly due to the fact that the overall hydrogen bond network is not destroyed after the proton transfer motion.

In this work, we employed the on-the-fly surface hopping dynamics simulation at the semiempirical OM2/MRCI level to understand the photoinduced processes for the BV phytochrome chromophore. Several critical geometries, including electronic state minima and conical intersections, were obtained. We found that the deprotonation of different pyrrole nitrogen leads to distinct photophysical and photochemical behaviors, such as the substantial changes in the excited state lifetime and nonadiabatic decay pathways.

The BV_N model shows the fastest excited-state decay because there are barrierless pathway from Franck-Condon region to the two CIs, and its excited-state process is dominated by the rotation of whole A-B double-rings moiety or whole C-D double-rings moiety. As the contrast, the excited-state population decay of the BV_C model takes place at a much longer time scale than that of BV_N model. The Ring D rotation and HOOP motion in the Ring B play efficient roles in the nonadiabatic process in the BV_C model. Interestingly, the other two models (BV_A and BV_B) show only extremely minor excited-state population decay within the simulation time duration.

We attribute the striking difference to the formation of different intramolecular hydrogen bond connective patterns, which affect the preferable reaction pathways and accessible conical intersections, and give diverse decay channels from the excited state to the ground state. At the end, we emphasize that in addition to the protein steric effect, the variation in deprotonation state also greatly impacts on the excited-state dynamics of the phytochrome chromophore. This finding not only deepens our understanding of the primary photoinduced reaction of the phytochrome chromophore, but also provides new design ideas for the efficient engineering near-infrared fluorescent proteins based on the phytochrome. Certainly, it is highly interesting to examine whether the protein environments can modify the protonation states of the BV chromophores and the protein pocket can regulate the photoisomerization dynamics. Clearly, the nonadiabatic dynamics mechanism of different phytochrome chromophore (BV, PΦB and PCB) is strongly modified by hydrogen bond networks, as well as their interactions with the surrounding protein and solution environments. Such theoretical simulations may involve the further large scale QM/MM nonadiabatic dynamics simulations. This represents the challenging and on-going topic in the future study.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

This work is supported by National Natural Science Foundation of China (NSFC, Nos. 21903030, 21933011 and 21873112).

Theoretical methods, BV chromophore, S1_min of four models, typical trajectories, reaction products of the BV_N model, distributions of selected internal coordinates for all hopping geometries in the BV_N model, time-dependent evolution of the distributions of a few key dihedral angles in the nonadiabatic dynamics, configurations at hop events, benchmark of the vertical energies, reaction pathways in the OM2/MRCI(24, 22) level, cartesian coordinates and additional implementation details are available. Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.cclet.2022.107850.