Simulation of two-dimensional infrared Raman spectroscopy with application to proteins

Two-dimensional infrared (2DIR) spectroscopy¹ is a powerful technique successfully applied to the amide-I band of proteins to reveal structural information. An alternative technique, which we will here denote two-dimensional infrared Raman (2DIRRaman) spectroscopy,² was first implemented experimentally in 1999.³ The technique has also been denoted Doubly Vibrationally Enhanced (DOVE) four-wave mixing spectroscopy³ or electron–vibration–vibration two-dimensional infrared (EVV 2DIR) spectroscopy.⁴ It has only recently been applied to the amide-I band of proteins,² and the interpretation of these spectra is still under discussion. This paper aims to implement a method to simulate the 2DIRRaman spectrum of the amide-I band and apply it to proteins dominated by α-helical and β-sheet type structures.

The 2DIRRaman technique was originally demonstrated for different solvent mixtures.^3,5 It was later applied to methylene spectra⁶ and to protein sidechain vibrations.⁷ In recent years, the technique has been demonstrated to be powerful in the application of determining drug binding.⁴ Furthermore, the technique was applied to the amide-I band of proteins,² which is well known to be sensitive to the protein secondary structure.⁸ The technique clearly has promising applications, and by combining infrared and Raman excitations, it is uniquely sensitive to couplings between infrared and Raman active modes.

While the theory for 2DIR spectroscopy is quite well developed with contributions from a broad community,^9–23 the theory for 2DIRRaman spectroscopy is still more limited²⁴ with a few applications to extensive systems.^25–28 Many of these approaches utilize techniques originally developed for 2DIR spectroscopy^1,29 or 2DRaman^30–35 spectroscopy. In this paper, we will present an implementation following the Numerical Integration of the Schrödinger Equation (NISE) approach in the NISE program, which was originally developed for 2DIR spectroscopy³⁶ and later extended to two-dimensional electronic (2DES) spectroscopy^37,38 and (two-dimensional) sum-frequency generation [(2D)SFG] spectroscopy.³⁹ The same time-domain approach was also applied by others for FTIR,⁴⁰ Raman,^11,40,41 SFG,^42,43 and 2DIR spectral calculations.^11,14,20,41 We will focus the application to the amide-I band of proteins, which was the focus of many 2DIR studies^1,44–46 and the recent experimental work with the 2DIRRaman technique.² We will show that the predicted spectra reproduce the key features observed in the recent experiments and these features can be understood considering the features observed in the corresponding Raman and infrared absorption spectra.

The remainder of this paper is organized as follows: first, in the section titled Methods, we will describe the fundamental theory behind the 2DIRRaman technique and the combined molecular dynamics–response function protocol used for simulations. We will then present and discuss the simulated spectra in the section titled Results. In this section, we will further discuss the implications of our findings. Finally, we will summarize the conclusions and present an outlook.

We started from the monomer protein data bank⁴⁷ structures 3V03⁴⁸ and 7MG1⁴⁹ for bovine serum albumin (BSE) and concanavalin A (ConA), respectively. The initial structures are shown in Fig. 1. The structures were first solvated with SPC/E water.⁵⁰ Sodium counterions were added to ensure charge neutrality. For BSA, this required 17 sodium ions, and for ConA, six sodium counterions were used. The molecular dynamics simulations were performed using GROMACS, version 2021,⁵¹ and employing the OPLS-AA force field⁵² for proteins. Equilibration was performed by a 100 ps NVT run and a subsequent 100 ps NPT run. The 1 ns production run was performed at 300 K and 1 bar using the velocity rescaling thermostat and the Parrinello–Rahman barostat.⁵³ Long-range interactions were accounted for using the particle-mesh-Ewald scheme⁵⁴ with a cutoff distance for short-range interactions of 1.0 nm. The simulation time step was 2 fs, and frames were stored every 20 fs for the spectral simulations.

For the spectral simulations, the amide-I vibrations are described with a vibrational Hamiltonian of the form¹ 

Here, i labels the amide-I sites, while $Bi†$ and B_i are the Bosonic creation and annihilation operators, respectively. The vibrational frequency of each local amide-I mode i is given by ω_i(t). The coupling between pairs of local modes is given by J_ij(t). The anharmonicity of local mode i is Δ_i(t). These properties are all time-dependent due to the protein dynamics and the interaction with the environment. The anharmonicity is needed for the 2DIR and 2DIRRaman spectra, where double excited states can be reached as a result of two consecutive excitations.

The amide-I Hamiltonian was constructed using AIM, which is a computer program recently developed to read in molecular dynamics trajectories from a broad range of molecular dynamics codes and generate vibrational Hamiltonians.^55,56 The Skinner map was used for the electrostatic frequency shift,⁵⁷ and pre-proline units⁵⁸ were modeled adding a 26 cm⁻¹ red shift to the Skinner map prediction⁵⁵ to reflect that these units are tertiary amide groups. The nearest neighbor couplings were modeled using the Glycine Dipeptide (GLDP) Ramachandran angle map.⁵⁹ Long-range couplings were treated with the dipole–dipole coupling approximation of the TDCTasumi map.⁶⁰ The Tsuboi⁶¹ transition-polarizability map was implemented in the AIM code⁵⁶ and used to describe the Raman transitions. The harmonic rule for the transition-dipole and transition-polarizability for sequence transitions was used, meaning that $μfe=2μeg$ and $αfe=2αeg$⁠. Here, g symbolizes the ground state, e is the single-excited state manifold, and f is the double-excited state manifold. The transitions from the ground state directly to the double excited state manifold were assumed to happen exclusively within the same amide-I unit and with the transition-dipole $μ⃗gf∝μ⃗ge$ based on previous mapping calculations.⁶² The proportionality factor was assumed to be identical for all transitions and set to one in the actual calculations as all signals under these assumptions are proportional to this factor. The anharmonicity was set to 16 cm⁻¹ following previous simulations.^1,63,64

The spectral simulations were performed with the NISE program,^38,65 where the response functions for the Raman spectra and 2DIRRaman spectra were implemented following a procedure similar to that used for infrared absorption and standard 2DIR spectroscopy.²² The response function for the Raman spectra is given by⁶⁶ 

The Raman spectrum is then obtained as the imaginary part of the Fourier transform over time t. The polarized (VV) and depolarized (VH) conditions were obtained by combining the appropriate components of the Raman tensors⁶⁷ and prefactors. The Raman technique was, thus, also implemented in NISE for this study.

The 2DIRRaman spectra were obtained from the response functions connected with the double-sided Feynman diagrams shown in Fig. 2,

The signal of the first diagram is of the rephasing type and can be isolated experimentally. The two other diagrams are of the non-rephasing type and cannot be separated from each other. Here, we will therefore present the rephasing spectra and non-rephasing spectra separately, where diagrams 2 and 3 are added up to give the non-rephasing spectrum. The two-dimensional spectra are obtained by a double Fourier transform over the time delays t₁ = τ₂ − τ₁ and t₃ = τ₃ − τ₂. Here, we will denote the corresponding frequencies ω₁ and ω₃. We note that various names have been used for both the time delays and frequencies in the literature. We based our choice on the ease to connect with traditional 2DIR spectra, which also allows using similar naming conventions in the simulation program. The different polarization components are obtained using the identical polarization averaging to conventional 2DIR spectroscopy.⁶⁸ This results in three polarizations in the isotropic situation corresponding to parallel polarization, perpendicular polarization, and cross-polarization.⁶⁹ With parallel polarization, we understand the situation, where both the polarization of the infrared interactions and the polarization of the visible ones are parallel to each other. With perpendicular polarization, we understand the situation, where the polarization of the infrared interactions is parallel to each other, but perpendicular to the polarization of the optical interactions. Denoting the Cartesian polarization components of the infrared fields in the frame of the MD simulation a and b and the Cartesian polarization components of the optical fields in that frame c and d, the expression for the parallel and perpendicular 2DIRRaman spectra are thus

The values of the weight factors $Aabcd‖$ and $Aabcd⊥$ are given in Ref. 68. The linear absorption spectra and 2DIR spectra were calculated using the same trajectories and settings using the existing implementation in NISE.⁶⁵ For the 2DIR spectra, the waiting time was set to zero, while the 2DIRRaman spectra do not have a waiting time. For all simulations, the coherence times were varied from 0 to 2.56 ps in 20 fs steps.

The linear absorption and Raman spectra for BSA and ConA are presented in Fig. 3. For BSA, all three spectra look very similar. They all contain a single broad feature centered at around 1630 cm⁻¹, which is slightly skewed to the blue side. Experimentally, a Raman peak position was reported at 1654 cm⁻¹ (Ref. 70), essentially coinciding with the 1650 cm⁻¹ reported for the infrared spectra.⁷¹ The ∼20 cm⁻¹ difference between the present calculations and experimental values is in good agreement with the systematic shift reported for the frequency mapping combination used here.⁶³ For ConA, the simulated Raman spectra peak at 1650 cm⁻¹, while in the infrared spectrum, a broad double peak structure is observed with features around 1610 and 1620 cm⁻¹. Experimentally, 1671 cm⁻¹ was reported for the Raman spectra,⁷⁰ while the infrared spectra exhibit a double peak structure with peaks at around 1620 and 1640 cm ⁻¹ (Ref. 72). Apart from a systematic shift, the simulated spectra also reproduce the key spectral features of ConA. The non-coincidence effect,¹¹ which is defined by a large difference between the line shape between the FTIR spectrum and the Raman spectra, is significant in ConA and very small in BSA. Our simulations reproduce this effect. A large non-coincidence effect may be a signature of the significant coupling between the local modes, leading to collective delocalized states^45,73 or a large variation in the transition-dipole moment, but not the transition-polarizability.¹¹ In the present situation, the former is most likely the explanation. The qualitative agreement for the FTIR spectra with experiments is good; however, we here refrain from a detailed quantitative comparison with explicit experiments as, for example, for ConA, significant variations between experimental spectra have been observed depending on the exact experimental conditions.^2,72,74 This may, among other things, depend on the potential amyloid formation.⁷⁵ 

In Fig. 4, we present the conventional 2DIR spectra of BSA and ConA. The spectra reproduce the experimental BSA⁷⁶ and ConA⁷² 2DIR spectra reasonably well, considering the spectral variation already discussed for the FTIR spectra. The BSA spectrum is dominated by a diagonally elongated diagonal peak and its corresponding overtone peak shifted by the anharmonicity below the diagonal. The ConA spectrum has more features; in particular, a clear cross-peak is observed above the diagonal. This is the tell-tale signature of the β-sheet structure.^44,72,77 The stripes seen in the lowest contour lines result from the ringing from the Fourier transform. This artifact can be reduced by extending the coherence time used and averaging over more samples.

We present the simulated rephasing 2DIRRaman spectra in Fig. 5. All spectra are dominated by a clear rephasing line shape feature. Little to no difference is observed between parallel polarization and perpendicular polarization spectra. The red main peak for BSA is clearly located below the diagonal line. This shape is in good agreement with experimental observations.² The experimental conditions for the spectra in figures 12(a) and 12(b) of that reference resemble the simulated conditions the best. The calculated spectra are all presented in the ideal instantaneous pulse limit, while the experimental spectra have a finite pulse duration and the spectra in Ref. 2 were obtained by varying the pulse frequencies instead of Fourier transforming the time delays.⁷⁸ For ConA, the red main peak clearly goes above the diagonal line and the upper blue sideband is stronger than the lower one. Again, this is in good agreement with experimental observations.² This observed difference between the BSA and ConA spectra can be understood as the last interaction in the 2DIRRaman spectra is a Raman interaction and the states observed on the ω₃-axis are expected to be Raman active, while the states observed along the ω₁-axis are infrared active. For BSA, the infrared and Raman active modes are all around the same frequency as there are essentially no non-coincidence effects in the α-helix system. On the contrary, the difference between the infrared and Raman spectra for ConA leads to a cross-peak above the diagonal. Compared to the conventional 2DIR spectra presented on the same scale in Fig. 4, the 2DIRRaman spectra are much broader. This is because, in 2DIR, the rephasing and non-rephasing contributions are added up, resulting in a purely absorptive spectrum,⁷⁹ which has narrow peaks. This is not possible in 2DIRRaman, where the peaks inherently have dispersive components.

In conventional 2DIR spectroscopy, polarization control is a powerful tool to study the structure^69,80 and dynamics.^81,82 2DIRRaman does not have an intrinsic waiting time during which the system is in a population state, and the method is, thus, not suited for studying dynamics. To examine the use of polarization for studying the structure, we looked at the rephasing 2DIRRaman signal integrated over ω₁ (see Fig. 6). In conventional 2DIR, the intensity ratio between the parallel and perpendicular spectra would be 3-to-1 if a single transition would be present and a value deviating from this would indicate the presence of multiple interacting transitions.⁷⁸ For the 2DIRRaman spectra calculated here, the parallel and perpendicular spectra have essentially identical intensities. It, thus, appears that little additional information will be revealed in the amide-I spectral region by varying the laser polarization. The observed intensity difference is also much smaller than that observed between the Raman VV and Raman VH spectra. This can be easily understood as for both the parallel and perpendicular 2DIRRaman signals, the polarizations of the optical fields are parallel to each other and the signal is sensitive to the diagonal part of the Raman tensor as is the case in the Raman VV spectra. The intensities plotted for BSA and ConA are normalized to the number of amide-I units, and the peak signal for BSA measured per amide-I unit is, thus, a bit higher than for ConA.

The simulated non-rephasing 2DIRRaman spectra are shown in Fig. 7. These spectra have both a blue peak, which originates from the third diagram in Fig. 2, and a red peak originating from the second diagram. Both peaks are elongated perpendicular to the black line, which is drawn as a diagonal line but displaced by 1600 cm⁻¹ along the ω₁-axis. This reflects the non-rephasing origin of the signal. The distance between the oppositely signed features along the ω₃-axis is expected to be given by the anharmonicity. The main difference between the spectra for the two proteins is that for ConA, both the red and the blue features are split in two clear peak features visible along the ω₁-axis, while for BSA, such a substructure is absent. The non-rephasing 2DIRRaman spectra, thus, exhibit a clear fingerprint of the β-sheet structure in ConA. As for the rephasing spectra, the difference between the parallel and perpendicular spectra is minimal.

The presented implementation of the 2DIRRaman technique makes use of the framework available in the NISE program.⁶⁵ The computational effort needed to calculate the 2DIRRaman spectra is dominated by the propagation of the double-excited states. This propagation was the main focus of the algorithm development for 2DIR spectra.^{20–22,83,84} The computational effort for the 2DIRRaman calculations in NISE, therefore, scales as N³ with the number of single excited states N, just as the conventional 2DIR calculations.²² 

In the present simulations, we assumed an isotropic solution. In the experiments,² a film was used. This may result in differences in the spectra. As BSA and ConA are known to form dimers and tetramers in solution, we also examined the rephasing 2DIRRaman spectra for such configurations (not shown), but we found the effect on the simulated spectra to be negligible. It is, however, important to realize that experimental conditions may have a crucial effect on the spectra as previously observed for the 2DIR spectra of ConA.^64,72,74

In the present implementation, we have assumed that the double frequency laser pulse only excites the overtone states and that the strength of this is proportional to the excitation probability from the ground state to the corresponding site basis single excited state. We expect this assumption to hold when the different site basis vibrations are not involving the same atoms. This approximation is, for example, expected to break down if one will treat the amide-I and amide-II vibrations simultaneously. In this case, the amide-I–amide-II combination vibration^85,86 on the same molecule can likely be directly excited. The inclusion of such excitations will require modification of the present code and more specific knowledge about the excitation cross section of the different states in the double excited state manifold.

In the present implementation, the 2DIRRaman spectra are calculated in the perturbative limit through the calculation of the response functions assuming instantaneous interactions with the applied electric fields. In practice,^2,3,6,7 depending on the experimental setup, the applied pulses will all be finite and the experimentally measured spectrum will depend on the details of the pulse shape and duration. This leads to pulse shape effects, which, in principle, can be accounted for in the perturbative limit, by appropriately modifying the Fourier transform used to convert the time-domain response functions to frequency domain spectra.^78,87

In this paper, we have described the implementation of the Raman and 2DIRRaman techniques in the NISE spectral simulation program and demonstrated it with the application to the amide-I vibration of two proteins. The simulated spectra were found to reproduce the key experimental features very well. We could understand these features in the 2DIRRaman spectra using the information on infrared and Raman active states. The key information revealed by the technique is about the coupling between these different states. The infrared and Raman spectra of the α-helix dominated protein essentially coincide, while the non-coincidence effect is substantial in the β-sheet dominated protein. Therefore, the 2DIRRaman spectra are substantially different for these two structural motifs and the technique can be used to distinguish them. In contrast to 2DIR spectroscopy of the amide-I band, we found that the laser polarization reveals no new information. The 2DIRRaman technique probes the double-excited states more directly; this likely makes the signal more sensitive to local effects as double-excited states are typically more localized due to the anharmonicity. This also complicates the interpretation. However, the already reported success in revealing couplings between drug targets and proteins demonstrates the potential of the technique.⁴ Furthermore, the observed couplings are not simply couplings between the fundamental infrared active and fundamental Raman active transitions. This highlights the importance of developing accurate simulation methods as presented here for the interpretation of the signals. Furthermore, 2DIRRaman spectroscopy may be a good additional test of the models already developed to describe vibrational spectra.¹⁶ While conventional 2DIR spectroscopy can be used for microscopy and imaging,^88,89 the Raman detection step in 2DIRRaman spectroscopy should make the technique more suitable for such types of applications.

The present implementations can be applied to study other proteins, but the implementation is, in principle, general and allows application to other systems and spectral regions, provided the underlying frequency mappings are known and the molecular dynamics simulations are performed. The simulation protocol can also be applied to model systems, where the Hamiltonian parameters are either assumed or fitted to a system of interest.

This publication is part of the project “Nanoscale Regulators of Photosynthesis” (with Project No. OCENW.GROOT.2019.086) of the research program NWO Groot, which was (partly) financed by the Dutch Research Council (NWO). The authors are grateful to Paul M. Donaldson for helpful and inspiring discussions and to Steven J. Roeters for useful suggestions.

The authors have no conflicts to disclose.

Carleen D. N. van Hengel: Formal analysis (lead); Investigation (lead); Methodology (equal); Software (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (supporting). Kim E. van Adrichem: Methodology (supporting); Software (supporting); Supervision (equal); Writing – review & editing (supporting). Thomas L. C. Jansen: Conceptualization (lead); Methodology (equal); Project administration (equal); Software (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.