Time-Dependent Perspective for the Intramolecular Couplings of the N-H Stretches of Protonated Tryptophan.

Quasi-classical direct dynamics simulations, performed with the B3LYP-D3/cc-pVDZ electronic structure theory, are reported for vibrational relaxation of the three NH stretches of the -NH3+ group of protonated tryptophan (TrpH+), excited to the n = 1 local mode states. The intramolecular vibrational energy relaxation (IVR) rates determined for these states, from the simulations, are in good agreement with the experiment. In accordance with the experiment, IVR for the free NH stretch is slowest, with faster IVR for the remaining two NH stretches which have intermolecular couplings with an O atom and a benzenoid ring. For the free NH and the NH coupled to the benzenoid ring, there are beats (i.e., recurrences) in their relaxations versus time. For the free NH stretch, 50% of the population remained in n = 1 when the trajectories were terminated at 0.4 ps. IVR for the free NH stretch is substantially slower than for the CH stretch in benzene. The agreement found in this study between quasi-classical direct dynamics simulations and experiments indicates the possible applicability of this simulation method to larger biological molecules. Because IVR can drive or inhibit reactions, calculations of IVR time scales are of interest, for example, in unimolecular reactions, mode-specific chemistry, and many photochemical processes.

Introduction

Vibrational spectroscopy is an important experimental means for obtaining structural information for biological molecules.[1−4] Of particular interest are the secondary structures of proteins and peptides, which are affected by their conformations, protonation sites, and hydrogen bonding. In recent work, IR–IR–UV hole-burning,[3] messenger tag vibrational predissociation,[4] and IR–IR hole burning[4] spectroscopies were used to investigate structures of protonated tryptophan (TrpH+) low-energy conformers (see Figure ). The positions and line widths of the N–H stretches of the −NH3+ group provided detailed structural information for the TrpH+ conformers, as shown in Figure . Separation of the conformers by hole-burning spectroscopy[3,4] allowed one to resolve individual contributions from the two conformers A and B being responsible for the splitting of the broad NH3+ lines that were found in earlier infrared multiphoton dissociation experiments.[5,6]

Figure 1

Optimized geometries of conformers A and B of TrpH+ determined from B3LYP-D3/cc-pVDZ calculations.

Figure 2

(a) Predissociation spectrum of singly H2 tagged TrpH+. (b,c) IR–IR hole burning spectra of singly H2 tagged TrpH+ at 3341 cm–1 (red arrow) and 3362 cm–1 (blue arrow), respectively (a–c adapted from ref (4)). (d) Calculated and 0.967 scaled harmonic spectrum of conformer B.

In the previous work,[3] conformer B for TrpH+ (Figure b) was assigned by comparison of experimental and calculated IR spectra. The Greek letters α, β, γ label the three hydrogen atoms in the protonated amino group. These labels indicate the dominant NH stretch in each of the three strong vibrational modes involving the −NH3+ group, denoted by NH(α), NH(β), and NH(γ), as shown in Figure . NH(α) is the mode with the largest amplitude of vibration of H toward the side group, NH(β) is influenced by the electrostatic bond of H to the oxygen atom, and NH(γ) is basically a free NH stretch (displacement vectors are shown in Figure S1 in the Supporting Information).

The experiments measure the absorption band, that is, the absorption probability versus energy, and the width of the absorption band reflects coupling between the excited eigenstate and other modes. This coupling may be considered to be a time-independent picture of intramolecular vibrational energy redistribution (IVR),[7−9] that is, the broader the bandwidth, the faster the IVR. All the measurements have in common that the bandwidth of the NH(γ) mode is smaller, although the NH(γ) mode has the highest energy. Computed with the Beyer–Swinehart algorithm,[10] the vibrational density of states at 3050 cm–1, for the NH(α) absorption region, is dense with ∼8.3 × 106 states/cm–1, while in the NH(γ) region at 3350 cm–1, it increases to ∼26.8 × 106 states/cm–1. However, NH(γ) is a free hydrogen stretch, and it is therefore not unexpected that the coupling strength to the other modes is weaker than in the NH(α) and NH(β) cases, where the H atom interacts with the ring structure and the oxygen atom, respectively (Figure ). The NH(α) and NH(β) modes have comparable bandwidths. Notably, the NH(γ) widths are very similar in all experiments, except for the conformer resolved Yale spectrum with a slightly larger width.

The abovementioned TrpH+ vibrational spectroscopies may be analyzed by a time-dependent picture of IVR. In this analysis, a superposition φ of the vibrational eigenstates ψ is considered,[9,11−13] often called a zero-order state. This is shown in Figure , where the initially prepared state |s⟩ is coupled to the remaining zero-order states |l⟩. The initially prepared state |s⟩ evolves in time, and the statistical limit for relaxation of |s⟩ occurs when |s⟩ is coupled to a sufficiently large number of |l⟩ that it appears to decay irreversibly and practically to zero, within typical experimental time scales. For this situation, the probability of remaining in |s⟩ versus time, |c(t)|2, is given by

Figure 3

Depiction of the initially excited superposition and zero-order state |s⟩, zero-order states |l⟩ coupled to |s⟩, and the exact eigenstates |n⟩ = C|s⟩ + ∑C|l⟩. The heights of the bars associated with the eigenstates |n⟩ states are equal to |C|2 (adapted from ref (9)).

The rate constant k is related to the vibrational spectrum for |s⟩ bywhere Γ is the full-width half maximum (fwhm) of the absorption envelope[9] and may be considered to be a rate constant for IVR.[9] The inverse of k is an IVR time τ. The vibrational eigenstates are shown by the |n⟩ states on the right side of Figure .

The experimental line width data from the different experiments, for the TrpH+ NH stretches, is summarized in Table together with the lifetime estimates (eq ). Some of the peaks are not conformer resolved (marked with a star in Table ), and therefore, the values for the peak width is larger in these cases. The lifetime for the NH(γ) mode is longer than those for the NH(α) and NH(β) modes. There is a tendency for the NH(β) mode to have a longer lifetime than NH(α).

Table 1

Experimental Line Widths of the Three N–H Stretch Modes of −NH3+ for TrpH+a

NH3+ mode	fwhm (cm–1)b	lifetime τ (fs)b	fwhm (cm–1)c	fwhm (cm–1)d	lifetime τ (fs)d	fwhm (cm–1)e
NH(α)	23.3 (21.1)	227.8 (251.6)	39.95* (38.6*)	18.6 (18.5)	285.4 (287.0)	39.9* (39.8*)
NH(β)	16.96 (13.9)	313.0 (381.9)	40.77* (39.5*)	11.1 (10.9)	478.3 (487.0)	13.7 (13.5)
NH(γ)	10.2 (5.2)	520.5 (1020.9)	7.9 (n/a)	7.4 (7.1)	717.4 (747.7)	7.3 (7.0)

Values in parentheses are deconvoluted assuming a Gaussian for the laser and a Lorentzian for the absorption band. fwhm of the laser was 7 cm–1 (Yale data) and 1.5 cm–1 (Lausanne data). Values with * are not conformer resolved. The lifetimes were determined from eq .

Yale: double resonance depletion spectrum of TrpH+(H2) of conformer B.[4]

Yale: predissociation spectrum of tagged TrpH+(H2) (all conformers).[4]

Lausanne: conformer B-specific IR spectrum of TrpH+.[3]

Lausanne: IR–UV gain spectrum of all conformers of TrpH+.[3]

Table 2

Parameters of the Morse Fits and Vibrational Levels of the Morse Local Modesa

	D (kcal/mol)	β (1/Å)	r0 (Å)	n = 0	n = 1	n = 2	n = 3
α	59.2379	2.58731	1.03998	4.42	12.75	20.40	27.36
β	89.6204	2.21424	1.04163	4.68	13.67	22.16	30.14
γ	114.740	2.10904	1.02518	5.06	14.83	24.15	33.01

1D potential energy curves for the α, β, and γ NH stretches of TrpH+, calculated with B3LYP-D3/cc-pVDZ theory, were fit with the Morse function. Energies for the vibrational levels are in kcal/mol.

For the work presented here, the quasi-classical trajectory (QCT) method[14] was used to calculate IVR rate constants k, eq , for relaxation of the N–H stretch local modes of the −NH3+ group of TrpH+. These rate constants are compared with those given by the absorption envelopes of the N–H stretches. The motivation of the current study is two-fold. One is to obtain information concerning the atomistic dynamics of the N–H stretch relaxations. The other is to assess the ability of the QCT method to represent IVR relaxations of these N–H stretches.

Simulation Method

Following the previous work,[15,16] the NH stretch superposition states for TrpH+ are identified as local modes. Extensive local mode dynamics have been studied for CH overtones of benzene, both experimentally[17−22] and computationally.[22−28] The computational studies have included both classical[23−25] and quantum[26−28] dynamics simulations. CH local modes have been investigated for other molecules[29] and of interest for the NH stretch local modes of TrpH+ is the NH stretch local mode of pyrrole.[30]

Previous studies of IVR for CH overtones in benzene,[24,27] energy flow in biomolecules,[31] and intramolecular dynamics within the HOSO2 complex[32] have shown that classical dynamics gives results in good agreement with quantum dynamics if quantum-like initial conditions with zero-point energy (zpe) in all modes are used for the classical simulations. Quasi-classical initial conditions[14] were chosen for a classical trajectory simulation of the n = 3 CH overtone of benzene,[24] and excellent agreement with the quantum dynamics[27] was obtained at short time. The classical simulation gave the correct overall absorption envelope but not the detailed structure inside the envelope.

The QCT method,[14] as described in ref (24), was used to choose initial conditions for the simulations reported here. A N–H stretch local mode state was prepared by first adding zpe to all of the TrpH+ vibrational modes, with random phases, and then, additional energy was added to the N–H stretch so that it was in the n = 1 quantum level. The NH(α) local mode thus consists of a background of zpe in all normal modes, on top of which a single quantum of excitation is put into the bond between N and H(α) by modifying the N–H(α) distance and the N–H(α) relative velocity. The NH(β) and NH(γ) local modes are prepared in a similar procedure.

This excitation prepared the initial φ = |ψ(0)⟩ superposition state. The potential energy of the N–H stretch was represented by a Morse function in selecting initial conditions for the local mode state. As described below, the B3LYP-D3/cc-pVDZ theory was used for the direct dynamics simulations, and this level of theory was used to calculate one-dimensional (1D) potential energy curves for the NH(α), NH(β), and NH(γ) stretches. These curves were fit by Morse potentials with the parameters, as listed in Table . The energy levels for the ground state and the first three overtones of the Morse potentials are also listed in Table and are very close to numeric solutions of the 1D Schrödinger equation for the potential energy curves. Energy for the Morse oscillator n quantum level was added with a random phase between the potential and kinetic energy.[24,33] Rotational energy was not added to TrpH+.

Direct dynamics simulations[34] of the relaxation of the n = 1 local mode states of TrpH+, with QCT initial conditions as described above, were performed with B3LYP-D3/cc-pVDZ electronic structure theory.[35−38] The rather small double-zeta basis was chosen to keep the computational costs at bay. In a previous direct molecular dynamics simulation of the TrpH+ IR spectrum,[4] this level of theory gave anharmonic vibrational frequencies, which is in good agreement with the experiment (TrpH+ harmonic vibrational frequencies, calculated at different levels of electronic structure theory, are compared in the Supporting Information). Conformer B was selected in this study because its optimized ground state is 27 meV, more stable than that for conformer A. The trajectories were integrated for 400 fs with a 0.4 fs time step using the velocity-Verlet algorithm.[39] In total, 100 trajectories were calculated for each of the three local modes. Convergence of the trajectory dynamics was checked by comparing the local mode relaxations for 50 and 100 trajectories. The nature of the relaxation dynamics was revealed by 50 trajectories, with the noise in the simulation results only slightly reduced with 100 trajectories. The trajectory simulations were performed with the VENUS chemical dynamics code interfaced with the NWChem electronic structure code.[40−42]

Simulation Results

The property determined from the simulations is the probability P(n,t) that TrpH+ remained in the n = 1 local mode state versus time. P(n,t) is the QCT analogue of |c(t)|2 in eq .[24] The quasi-classical binning model was used to identify if the trajectory was in the local mode state n = 1; that is, the trajectory was assumed to be in n = 1 if the energy for the local mode was in the energy interval [E,E] of the initially prepared local mode. This information was collected versus time for each trajectory and then merged for all the trajectories to form P(n,t).

The P(n,t) determined from the simulations are plotted in Figure . Relaxation is fastest for the NH(β) local mode and slowest for the NH(γ) local mode. Of much interest is that the relaxation dynamics differ for the local modes. As described below, for NH(β), P(n,t) is exponential. In contrast, for the NH(α) and NH(γ) local modes, the P(n,t) are nonexponential with beats (i.e., recurrences), which are more prominent for the longest lived NH(γ) local mode. Each of the local modes has significant population when the trajectories are terminated at 400 fs, which is 0.5, ∼0.3, and ∼0.2 for NH(γ), NH(α), and NH(β), respectively.

Figure 4

Probability vs time, P(n,t), of populating the initially excited n = 1 states of the NH(α), NH(β), and NH(γ) local modes of TrpH+. The 2σ intervals (light colors) of the mean values were obtained from bootstrapping for every time step.

The simulation P(n,t) were fit with the function used previously,[24] that is

This function allows for recurrences with two different frequencies while ensuring an overall exponentially decreasing character by the condition that k > 0 if a ≠ 0. With eight parameters in total, it is flexible enough for our purpose. The absorption spectrum of the exponential a1e– only would give a perfectly Lorentzian line shape.[24] For NH(β), P(n,t) is exponential with fit parameters a1 = 0.856 and k1 = 0.00297 fs–1. The P(n,t) parameters for NH(α) are a1 = 0.907, a2 = 0.048, k1 = 0.00297 fs–1, k2 = 0.00200 fs–1, and ω2 = 0.0378 fs–1. Relaxation of NH(α) is approximately exponential with a single beat. For NH(γ), the P(n,t) parameters are a1 = 0.795, a2 = 0.315, a3 = −0.114, k1 = 0.0010 fs–1, k2 = 0.0108 fs–1, k3 = 0.00206 fs–1, ω2 = 0.0216 fs–1, and ω3 = 0.0180 fs–1. Relaxation of NH(γ) is decidedly nonexponential with two beats.

The width of the absorption envelope for a superposition state is well approximated by the rate constant for the initial exponential decay.[6,9,27] Accordingly, rate constants for the simulations’ initial exponential decays for the NH local mode states were used to determine their IVR lifetimes and compare with the experimental IVR lifetimes, as listed in Table , determined from the experimental line widths. The simulation IVR lifetimes are 211, 337, and 1000 fs for the NH(β), NH(α), and NH(γ) local modes, respectively. There is no quantitative agreement between the simulation and experimental lifetimes, but there is overall good agreement. Except for the experimental IVR lifetime of 487 fs for NH(β),[3] the simulation and experimental lifetimes agree within a factor of two. Both the simulations and experiments indicate that the free NH(γ) local mode has the longest IVR lifetime. However, the simulations find that the IVR lifetime is longer for NH(α), as compared to NH(β), while the experiments find the opposite.[3,4] However, differences in the IVR lifetimes of these two local modes are not substantial.

Conclusions and Future Extensions

There are important conclusions and future extensions for the current simulation study. The QCT simulations reported here is in overall good agreement with the experiment for the IVR lifetimes of the NH stretch local mode states of TrpH+. Such agreement was found previously in simulations of vibrational relaxation of superposition states for benzene,[24,27] biomolecules,[31] and HOSO2.[32] The simulations were performed by direct dynamics,[34,43] with the method of quasi-classical chemical dynamics directly interfaced with electronic structure theory. This simulation methodology is applicable to a broad range of molecules and superposition states. An important component for the simulations is the inclusion of zpe in the trajectory initial conditions to represent the molecule’s density of states for the initial relaxation of the superposition state.[44,45] Without zpe, classical dynamics becomes more regular suppressing IVR. This has been illustrated in detail for relaxation of the CH local mode in benzene.[25]

There is a striking difference between the IVR lifetimes for the n = 3 CH local mode of C6H6 and the free NH(γ) n = 1 local mode of TrpH+. The IVR lifetime for the C6H6 local mode, determined from the width of the absorption envelope from both classical[24] and quantum[27] dynamics, is 62.5 fs. The IVR lifetime for the TrpH+ local mode is 1000 fs from the current simulation and 1020.9, 747.7 fs from the experiment.[3,4] Rapid IVR for the C6H6 local mode is the result of a 2:1 Fermi resonance between the C–H stretch and the adjacent HCC bend.[23,26] It is of interest to determine IVR pathway(s) for the free NH(γ) local mode of TrpH+.

There are many possible extensions of the current simulation study. Of interest are the pathways for energy flow from the excited NH local modes. For example, how do the intermolecule couplings for the NH(α) and NH(β) local modes affect energy flow? This information may be obtained by following energy transfer to the TrpH+ normal modes versus time.[46] If the simulation P(n,t) for the local modes were extended to longer times, absorption spectra for the local modes could be determined.[24] It is clearly of interest to determine the couplings leading to the beats (recurrences) in the P(n,t) for the NH(α) and NH(γ) local modes. The current simulations investigated IVR for the n = 1 NH local modes of TrpH+. The simulations could be extended to study IVR for n > 1. There is always the question regarding how direct dynamics results depend on the level of electronic structure theory used for the simulations, and different electronic structure theories could be considered to study the relaxation of the NH local modes of TrpH+.

This work presented here may be the first quasi-classical direct dynamics simulation of IVR for a molecule as large as TrpH+. The agreement found in the experiment holds the promise of using this simulation approach to study IVR in other biological molecules. In future work, it would be of interest to consider other simulation methods based on classical dynamics. Fixed energy Wigner sampling[47] could be used to sample the zpe level of TrpH+ instead of the QCT method used here. Time-dependent semiclassical dynamics[48,49] would introduce quantum effects into the simulations, but because of their computation cost, they are not possible by direct dynamics. One approach would be to use direct dynamics to construct an analytic and/or interpolated potential energy function for TrpH+ and use it for the semiclassical dynamics. This would be a challenging study.