On the Possible Role of an Intermolecular Charge Transfer State in the Excitation of the Biologically Active Bond of the Retinal Chromophore- counterion Pair

In non-polar solvents the protonation of the all-trans Schiff base of retinal (SBR+) using strong acids leads to formation of retinal chromophore-counterion pairs. Previously we have shown that the main non-reactive deactivation channel of these ion pairs involves internal conversion from the initially excited ππ* state to an inter-molecular charge transfer state (inter-CT) with subsequent dissociation of the chromophore-counterion pair. In solution this leads to the abortion of isomerization. Motivated by the recent X-ray diffraction experiments showing that the disruption of the hydrogen-bonded network of counterions precedes the isomerization of all-trans SBR+ in bacteriorhodopsin we decided to take a closer look at the dynamics of the chromophore-counterion pair in the inter-CT state. Using constrained non-adiabatic dynamics simulations in which the dissociation of the chromophore-counterion pair was impeded, we show that the charge distribution in the inter-CT state leads to site-specific elongation of the biologically active C13=C14 bond. On this basis we hypothesize that an inter-molecular charge transfer state involving the retinal chromophore and the H-bonded counterions (S2) may play an active role in the photoisomerization reaction in a constraint environment.


INTRODUCTION
HE photoisomerization of the all-trans Schiff base of retinal (SBR + ) in bacteriorhodopsin (bR) is one of the most studied photochemical reactions.  The reaction is ultrafast, very efficient (~65 %) and happens specifically around the C13=C14 double bond. [20,[28][29][30] There is general agreement that after photoexcitation of all-trans SBR + to the first excited electronic state (S1) charge migration takes place within the chromophore with positive charge migrating toward the β-ionone ring and negative one toward the Schiff base (SB) giving rise to an intra-molecular charge transfer (intra-CT) state. [23,31,32] However, how the reaction proceeds on the intra-CT state is an open question. Namely, in a very recent experimental and computational work [33] the ubiquitous two-state model has been replaced by the long conjectured three-state model [13,34,35] of photoisomerization in bR. Valuable new insight has also been obtained from the X-ray diffraction study of Nogly et al. in which the sub-picosecond response of the hydrogenbonded network of counterions surrounding the photoexcited all-trans SBR + was reported. [17] The time evolution of difference Fourier electron density light dark obs obs showed that the disruption of the H-bonded network surrounding all-trans SBR + precedes the formation of the 13-cis SBR + isomer.
Apparently there is little similarity between the photoisomerization of all-trans SBR + in bR and the reaction in solution. [36][37][38] The reaction in solution is an order of magnitude slower, inefficient and yields a mixture of 9-cis SBR + , 11-cis SBR + and 13-cis SBR + stereoisomers. [37,[39][40][41] The mechanisms leading to the non-reactive deactivation in solution have been investigated and it was found that two channels, a fast one operating on the subpicosecond to picosecond time scale and a slow one active on a time scale of several picoseconds, contribute to the non-reactive decay of the S1 state. [41][42][43][44] In non-polar solvents in which the electrostatic interaction between all-trans SBR + and the surrounding counterions is favored, the fast non-reactive channel is the dominant one. [41] In two joint experimental and theoretical studies we have uncovered the molecular mechanism behind the fast non-reactive channel. [32] First, by combining a 2D correlation analysis of the experimental infrared spectra with multidimensional anharmonic frequency calculations, [45,46] we have shown that chromophore-counterion pairs are indeed formed in the non-polar dichloromethane (DCM). Depending on the ratios of all-trans n-butylamine SBR + (nSBR + ) and the protonating trifluoroacetic (TFA) acid, two types of ionic aggregates were formed: [nSBR + TFA -] and [nSBR + TFA -]TFA. [47] Subsequent pump-probe measurements demonstrated that the excited-state dynamics slows down by 4 % in going from [nSBR + TFA -] to [nSBR + TFA -]TFA. From the computational side we have shown that the fast non-reactive decay involves a formation of an inter-molecular charge transfer (inter-CT) state via electron transfer from the counterion to the retinal backbone. [32] Altogether, previous studies showed that the electronic coupling between all-trans nSBR + and the counterions affects the reaction.
As the dissociation of the chromophore-counterion pairs occurs both in the protein binding pocket [17] and in solution, [32] we thought worthwhile to study this process in more detail. While the [nSBR + TFA -]TFA ion-pair is a very simplified model of the environment that all-trans SBR + experiences in bR, we believe that new insight into the steps preceding photoisomerization can be obtained from a thorough investigation of the dynamics in the inter-CT sate. Thus, from the relatively large body of non-adiabatic dynamics trajectories of [nSBR + TFA -]TFA (73) reported in Ref. [32] we selected representable trajectories that ended up the inter-CT state (27/73) and performed constrained non-adiabatic dynamics simulations in which the distance between nSBR + and the TFAcounterion was fixed. By prompting the system to reside for a longer time in the inter-CT state we could compute the average length of all double and single bonds in the chromophore and unravel the characteristic structural changes. We will see shortly that this includes a pronounced elongation of the C13=C14 bond.

COMPUTATIONAL METHODS
Non-adiabatic dynamics simulations of the photoexcited chromophore in solution were performed using the ONIOM(QM/MM) approach. [48][49][50] For a given molecular geometry the excitation energies were computed as n 0 ,ONIOM ,ONIOM ONIOM high low low model real model where S0 and Sn are to the electronic singlet ground and n th excited state, respectively. The chromophore-counterion pair [nSBR + TFA -]TFA constituted the model system. The energy of the model system, which is a subsystem of the real (entire) system is computed at a high (QM) level of theory and denoted high model

ΔE
. The solvent (real -model system) was described using a low level (MM) method. As the electronic excitations were localized on the chromophore-counterion pair, that is on the model system we used the constrained low-level state (CLS) approximation [49] and computed the ONIOM excitation energy as In other words, excited state calculations were performed only for the model system.
As a necessary compromise between accuracy and feasibility we performed non-adiabatic dynamics simulations in the framework of TDDFT using the hybrid PBE0 functional [51,52] and the SVP basis set. [53] Tully's fewestswitching surface hopping (FSSH) algorithm was used to switch between electronic states. [54] The non-adiabatic couplings, which are required in FSSH calculations, were computed by constructing a formal CIS wave function using Kohn-Sham spin-orbitals and evaluating the overlap between singly excited Slater determinants. [55][56][57][58] The CIS and molecular orbital coefficients were obtained from Turbomole 7.0. [59] The MM part was described using the general Amber force field (GAFF), [60] both in the ground and in the excited states. In all simulations periodic boundary conditions were used.
For the constrained QM/MM non-adiabatic dynamics simulations of [nSBR + TFA -]TFA in DCM we used the RATTLE algorithm. [61] The constraint was imposed once the system reached the inter-CT state, which happened in 38 % of unconstrained trajectories. The N-O distance was then fixed at the value attained in intra-CT/inter-CT conical intersection (CI) geometry. Thus, the non-adiabatic dynamics trajectories analyzed in this work differed from those in Ref. [32] for the time period that the system spent in the inter-CT state.
The initial conditions for 73 [nSBR + TFA -]TFA nonadiabatic trajectories were sampled from two room temperature, 10 ps long RI-PBE/TZVP//GAFF NVT trajectories in the ground electronic state. The PBE0/SVP//GAFF NA trajectories were propagated in a subspace spanned by the ground and three singlet excited electronic states. All trajectories were initiated from the first excited state and the trajectories were propagated for 800 fs or until the gap between the S1 and S0 states dropped below 0.15 eV. The classical equations of motion were integrated in steps of 0.5 fs using the velocity Verlet algorithm, while the Runge-Kutta method was used to integrate the time-dependent Schrödinger equation in time steps of 5 × 10 -5 fs.
Special attention was paid to the validation of the results of the dynamics. The ADC(2) method was used to assess the applicability of the PBE0/SVP//GAFF method. The validation of computational method is extensively discussed in section 3.4 of Ref. [32].

RESULTS AND DISCUSSION
We start by reviewing the salient features of non-adiabatic dynamics simulations performed in Ref. [32]. As most experiments were performed using an excess of acid, we focused on a rather large system encompassing nSBR + and two TFA counterions, [nSBR + TFA -]TFA. It was found that 28 out of 73 [nSBR + TFA -]TFA trajectories in DCM deactivated non-reactively to the ground state. All but one relaxed via a counterion-sensitive mechanism that we identified as the fast non-reactive channel. The remaining trajectories stayed in the S1(ππ*) state. Figure 1 displays a characteristic [nSBR + TFA -]TFA non-adiabatic trajectory relaxing via the fast non-reactive channel. The time variation of the potential energy of the ground (black line) and first two excited electronic states (red and blue) is shown. Red circles indicate the currently populated electronic state.
In the initial geometry the β-ionone ring is almost coplanar with respect to the nSBR + backbone. From the insets showing the differences in the electron density between S1 and S0 states one can infer the change of the character of the S1 state. In the initially excited ππ* state the excitation is delocalized along the backbone (leftmost). This sudden change in the electron density leads to the elongation of double bonds and contraction of single bonds. That is why in the first 50-60 fs electronic energies of all states evolve in a very similar way (see also Figure 2). With time, the dynamics in the S1 state leads to the strengthening of the intra-CT character of S1 (left). There is no dissociation of the ion-pair as the distance between the and the first two excited singlet states (S1, red; S2, blue) for the selected non-adiabatic trajectory. The circles indicate the populated state at a given time. The insets show the electron density difference between the S1 and S0 states with areas of depletion (increase) of electron density in S1 shown in red (blue). The character of the S1 state changes from ππ* (left) to intra-CT (middle) and inter-CT (right). The switch from intra-CT to inter-CT state at 130 fs (vertical dashed line) is followed by the dissociation of [nSBR + TFA -]TFA. (b) Change of the chromophore-counterion distance, dNO, during the dynamics.
nitrogen and the closest oxygen of the TFA counterion is virtually unchanged (Figure 1b). At 130 fs a CI between two excited states is reached and internal conversion takes place. From the change in the electron density maps it is evident that a hole (red areas) has been translocated from the β-ionone ring to the carboxylic group of the TFAcounterion. Alternatively, this can be viewed as a transfer of negative charge from the counterion to the retinal backbone. The net result is an increase of electron density (blue areas) at the Schiff base side of the chromophore (right and rightmost). Clearly the new electronic state corresponds to an inter-molecular charge transfer (inter-CT) state. The dissociation of the hydrogen bond between the chromophore and the counterion (Figure 1b) takes place in the inter-CT surface. This strongly destabilizes the formally ionic ground electronic state of [nSBR + TFA -]TFA and leads in the next ~100 fs to the CI with the ground state. The deactivation is non-reactive.

Constrained Non-adiabatic Dynamics
Having illustrated how the dynamics in the inter-CT state leads to non-reactive deactivation of all-trans nSBR + in DCM, we explore the effects of the nuclear motion on that potential energy surface for systems with frozen chromophore-counterion distance. Two orientations of the β-ionone ring will be considered. In the first case the β-ionone ring is almost coplanar with the retinal backbone. As the situation corresponds to the trajectory displayed in Figure 1. we start by imposing a constraint to that trajectory.
Thus, the non-adiabatic trajectory shown in Figure 2 differs from Figure 1 only after 135 fs, when the constraint was imposed (vertical dashed line). One immediately sees that from 135 onward the system experiences a strong and prolonged vibrational excitation. Compared to the unconstrained case the lifetime of the inter-CT state is 150 fs longer. This allows us to investigate in detail the change of the bond structure. The dynamics in the inter-CT state leads to vibrational excitation of the two (formally) double bonds at the ends of the π-conjugated system, C 5=C6 and C13=C14 (Figure 2b). The effect is directly related to the hole translocation from the β-ionone ring to the TFAcounterion. On the β-ionone side, one observes a sudden contraction of the C5=C6 bond by almost 0.08 Å as well as an elongation of the C13=C14 bond to a maximum value of 1.49 Å and on average by 1.43 Å. Subsequently, the vibrational excitation is converted into torsion of the β-ionone ring, varying between 7° and 70°, and redistributed along the retinal backbone. The flexing of the retinal backbone leads to nonreactive deactivation to the ground state. Figure 2. Localization of vibrational excitation in the nSBR + chromophore induced by hole translocation. (a) Potential energy profiles of the electronic ground state (S0, black) and the first two excited singlet states (S1, red; S2, blue) for the non-adiabatic trajectory analysed in Figure 1. The constraint is imposed at 135 fs (vertical dashed line) and the nSBR + to TFAdistance is fixed to dNO = 2.90 Å. The insets show representative [nSBR + TFA -]TFA geometries along the trajectory. (b) Time variation of the C5=C6 (grey) and C13=C14 (sky blue) bond lengths.
As a second example we selected a trajectory in which the β-ionone ring is almost perpendicular to the retinal backbone at the intra-CT/inter-CT CI. Note that initially the chromophore was in the 6s-cis form, which is more stable than the 6s-trans conformer in DCM. For consistency we reanalyze the trajectory shown in Figure 3 of Ref. [32], but of course under constraint conditions. The constrained trajectory is displayed in Figure 3. The constraint is imposed at 216 fs (vertical dashed line) when the system crossed to the inter-CT state. Here again the contraction of the C5=C6 bond (0.06 Å) and the elongation of the C13=C14 bond (0.04 Å) in the inter-CT state is clearly visible. However, as the β-ionone ring is effectively decoupled from the conjugated double bonds structure, the vibrational excitation of the C5=C6 bond persists in time.
In the two constraint non-adiabatic trajectories shown in Figures 2 and 3 we highlighted the weakening of the biologically active C13=C14 double bond. This creates the conditions for regiospecific isomerization and one would expect that the isomerization about C13=C14 may take place. However, in solution the vibrational excitation of the C13=C14 bond is transferred into the unimpeded largeamplitude motion of the all-trans nSBR + frame (Figures 2  and 3, rightmost inset). Flexing of the nSBR + skeleton destabilizes the ground state to such an extent that the S1-S0 energy gap closes down and the system returns nonreactively to the ground state. This happens, however, at longer times than in the unconstraint dynamics.
The most important effect observed in all non-adiabatic trajectories deactivating via the inter-CT state is the elongation of the C13=C14 bond. The time and ensemble averaged bond lengths of nSBR + in the inter-CT state, and the differences with respect to the intra-CT state are shown in Figure 4. With an average bond length of 1.42 Å, C13=C14 is undoubtedly the weakest double bond in the chromophore and consequently the most likely isomerization site. Note also the large relative change that the C13=C14 bond experiences when going from the intra-CT to the inter-CT state (Figure 4b).
Owing to the very different outcome of the all-trans SBR + photoisomerization reactions in bR and in non-polar solvents, the comparison between the two processes is difficult. We found, however, some interesting similarities in the initial stages of the photoisomerization process. In solution the dissociation of the chromophore-counterion pair and the selective activation of the C13=C14 bond are delayed and take place once the system crosses to the inter-CT state. This is an optically dark state as the oscillator strength decreased by about three orders of magnitude when going from intra-CT to inter-CT. This finding may be related with the observations of Ruhman et al. [22] that the population of the fluorescent state is biologically inactive and that the torsion around C13=C14 occurs only after the decay of the fluorescent state. In addition, from Nogly et al. [17] we know that in bR the isomerization takes place only after the disruption of the counterion network. Thus, one cannot a priori eliminate the possibility that a delocalized inter-molecular charge transfer state may play a role in the photoisomerization.
In the present approach nuclear quantum effects have been neglected. To address this issue, reduced dimensionality models could be constructed, [62] and the importance of quantum effects in the H-bonded nSBR + TFAmoiety could be estimated by comparing approximate results obtained with surface hopping or the density matrix evolution methods to exact quantum results. [63,64] We also note that instanton theory is a very promising approach for computing tunneling rates. [65,66]

SUMMARY
The non-reactive deactivation of the [nSBR + TFA -]TFA complex in non-polar solvents proceeds through several stages: the initial population of the lowest ππ* state, the separation of charge within the retinal chromophore, internal conversion to an inter-CT state from where the dissociation of the H-bonded chromophore-counterion pair takes place.
To study the consequence of the motion on the inter-CT state we performed constrained QM/MM nonadiabatic dynamics simulations of [nSBR + TFA -]TFA. The longer time that the constraint system spends in the inter-CT state allowed us to compute the average length of all double and single bonds in the chromophore with a relatively good statistic. This allowed us to demonstrate that selective elongation of the biologically active C13=C14 takes place in the inter-CT state.
Thus, the population of the inter-CT state of the chromophore-counterion not only provides a framework for understanding the dissociation of the H-bonded ionpair, but also provides a mechanism for explaining the elongation of the biologically active bond in the chromophore. On this ground we speculated whether an intermolecular charge transfer state may also play a role in the isomerization of all-trans SBR + in bR. Future nonadiabatic dynamics simulation of all-trans SBR + isomerization in bR which would account for the relevant counterions at the quantum mechanical level may give an answer to this query.
Notes. The authors declare no competing financial interests.