Analysis of Vibrational Circular Dichroism Spectra of Peptides: A Generalized Coupled Oscillator Approach of a Small Peptide Model Using VCDtools.

Vibrational circular dichroism (VCD) is one of the major spectroscopic tools to study peptides. Nevertheless, a full understanding of what determines the signs and intensities of VCD bands of these compounds in the amide I and amide II spectral regions is still far from complete. In the present work, we study the origin of these VCD signals using the general coupled oscillator (GCO) analysis, a novel approach that has recently been developed. We apply this approach to the ForValNHMe model peptide in both α-helix and β-sheet configurations. We show that the intense VCD signals observed in the amide I and amide II spectral regions essentially have the same underlying mechanism, namely, the through-space coupling of electric dipoles. The crucial role played by intramolecular hydrogen bonds in determining VCD intensities is also illustrated. Moreover, we find that the contributions to the rotational strengths, considered to be insignificant in standard VCD models, may have sizable magnitudes and can thus not always be neglected. In addition, the VCD robustness of the amide I and II modes has been investigated by monitoring the variation of the rotational strength and its contributing terms during linear transit scans and by performing calculations with different computational parameters. From these studies-and in particular, the decomposition of the rotational strength made possible by the GCO analysis-it becomes clear that one should be cautious when employing measures of robustness as proposed previously.

Introduction

The study of peptides using vibrational circular dichroism (VCD) spectroscopy has been the focus of many studies since the introduction of this powerful technique.[1−3] These studies have established that the amide I and amide II bands provide important insight into the dominant and/or fractional average content of the secondary structures of peptides and proteins. Traditionally, the analysis of these bands was done using the coupled oscillator (CO) model.[4] However, after the development of density functional theory (DFT) programs that can accurately and routinely compute VCD spectra, more elaborate computational methods have been developed that now allow one to compute VCD spectra for very large systems.[5−8]

Subsequent computational studies on peptides and proteins[6,9−15] have led to two important conclusions. First, it was concluded that the normal modes associated with the amide I and amide II bands are delocalized modes involving motion of several interacting groups. The general consensus is, therefore, that the observed amide I and amide II VCD signals are caused by distance effects associated with the coupling of groups on different strands. In the case of the amide I VCD bands, a number of computational studies[6,12,13,16] have demonstrated that the very intense VCD bands[17,18] originate from a more general degenerate CO mechanism. The amide II VCD bands, on the other hand, are typically not interpreted in terms of a CO model even though they appear to be governed by a similar mechanism. Second, it has been found that it is more difficult to identify β-sheets than α-helices. This is thought to be a consequence of the fact that in β-sheets, the amide groups that are coupled have almost the same planar arrangement and thus often may not form a chiral arrangement.[14,15]

To investigate this last aspect, Magyarfalvi et al. recently performed a VCD robustness analysis of the amide I and amide II modes for a series of For(AA)NHMe (AA being Val, Asn, Asp, or Cys, n = 1–5 for Val and Asn, and n = 1 for Asp and Cys) model peptides with α-helix and β-sheet backbone conformations.[19] Apart from a few exceptions, it was concluded that for short peptide models, the computed VCD sign of the amide I and amide II modes is most often nonrobust.[20] On the other hand, in longer peptide models, the amide I and amide II modes were found to be robust if they involved couplings between similar amide vibrations localized on the various backbone amide groups. These findings seem to corroborate the experimental observations but do not provide any physical insight into why it is difficult to accurately compute the VCD sign of nonrobust modes. As discussed previously,[21] this is because the performed VCD robustness analysis—which is based on the values computed for the dissymmetry factor[20]—is unable to judge properly the stability of the VCD sign computed for modes in which the dominant contribution to the total VCD intensity is determined by through-space coupling of electric dipoles.

To tackle this issue, we employed, in this work, the general CO (GCO) formalism, a novel approach for analyzing the computed VCD intensities that we have recently introduced.[21,22] The GCO expression is an exact formulation of the VCD CO model which is applicable for all types of normal modes. In addition, the GCO analysis can be used to judge the stability of the computed VCD signs. This makes the GCO analysis ideally suitable for studying the mechanisms that determine the magnitude and sign of the VCD signals observed in the amide I and amide II spectral regions. To make the GCO analysis even more efficient and useful, we have embedded it into the graphical user interface (GUI) of the Amsterdam density functional (ADF) software package.[23,25] This has allowed us to implement a novel normal mode analysis, which makes it possible to monitor on-the-fly the nuclear displacement vectors and the atomic contributions of the vibrational electric DTMs.

In this work, we will use these newly developed tools to analyze ForValNHMe model peptides that have previously been studied by Magyarfalvi et al.[19] The relatively small size of the system allows us to perform a large number of careful calculations on the dependence of the amide I and amide II VCD signals on structural and electronic parameters, as well as detailed normal mode analyses. Clearly, such calculations do not aim to interpret the sign and intensity patterns of VCD spectra of large peptides and proteins. What they do allow us to show is that the newly developed tools enable us (i) to pinpoint the origin of the VCD signal and, at the same time, judge the stability of its sign, (ii) to explain the differences observed between the VCD spectra computed for α-helix and β-sheet conformers, and (iii) to understand why some of the normal modes have a nonrobust VCD sign. Importantly, they also allow us to get a detailed insight into the VCD GCO mechanism. As a result, it will become clear that GCO analyses such as the ones performed in the present work on peptides are much more generally applicable and are in fact an indispensable tool for unraveling the origin and sensitivity of VCD signals in any (supra)molecular system.

Theory

A detailed description of the GCO formalism can be found in ref (22). Here, it will be briefly summarized. The IR and VCD intensities of the fundamental vibrational transition of a normal mode j are given by the dipole (D01) and rotational (R01) strengths[24]where E⃗01tot(j) and M⃗10tot(j) are the total electric and magnetic dipole transition moments (DTMs), respectively, associated with the fundamental vibrational transition (01) of mode j while i is the unit imaginary number.

Fragment Decomposition of the Rotational Strength

In the harmonic approximation, E⃗01tot(j) and M⃗10tot(j) can be written as a sum of atomic components.[24] Consequently, by grouping the atoms of a molecule in three molecular fragments (labeled as A, B, and R) one can write

In eqs and 4, E⃗01 and M⃗01 are the EDTM and MDTM associated with a molecular fragment (X is A, B, or R), λ runs over nuclei, N, N, and N are the number of atoms in the molecular fragments A, B, and R, respectively, while N is the total number of atoms in the molecule (N = N + N + N). By inserting eqs and 4 in eq , we obtain[22]

R01 is the generalized CO term. This term is associated with the through-space interaction of fragments A and B and is given by

The R01 term represents the contribution of the individual fragments A and B

Finally, the R01(j) contribution contains the coupling between fragment R and fragments A and B, and the individual contributions from fragment R

The fragment decomposition of the rotational strengths defined by eqs –8 is exact within the harmonic approximation and can be performed for any type of normal mode in both symmetric and asymmetric chiral molecules. The A and B fragments are the important fragments and as such contain atoms that exhibit large EDTMs and MDTMs. Fragment R contains the remaining atoms, that is, atoms with small or zero EDTMs and MDTMs, and as a result, its contribution to the total rotational strength is expected to be small. Indeed, as shown, one can almost always define the A, B, and R fragments such that the R01 term is very small and negligible, irrespective of the type of the normal mode.[22] More importantly, however, we found that large VCD signals are always due to dominant coupling R01 terms.

GCO VCD Analysis

The fragment decomposition of the rotational strength defined in the previous section allows one to interpret VCD signals in terms of the interaction of two molecular fragments and to identify two sources that contribute dominantly to the rotational strength of a particular transition. The first contribution (R01)—which is expected to be small when the two fragments are achiral—is associated with the two isolated fragments and the second (R01) with the through-space coupling of their EDTMs and MDTMs. This perspective offers more insight into the origin of the VCD signal than eq , although a simple and intuitive interpretation for the magnitude and sign of R01 remains practically impossible because the MDTMs are origin-dependent and imaginary. Recently, it has been shown[22] that more insight into the properties of R01 can be gained by casting eq into the “coupled oscillator” form[4]where ν is the frequency of mode j, c the speed of light in vacuum, and Y⃗(j,A,B), the generalized CO vector, which is defined per normal mode and depends on the definition of the GCO fragments A and B. Unlike eqs and 6, which depend explicitly on the MDTM, the second equality in eq allows us to interpret the VCD intensities in terms of couplings between the EDTMs associated with the GCO fragments A and B, which are basically charge movements separated by the Y⃗ distance. As we have shown,[21,22,25] such an analysis provides detailed insight into the magnitude and sign of the calculated VCD intensities.

A key step in the VCD GCO analysis is the identification of the molecular fragments A, B, and R, providing the most simple and intuitive physical insight. The procedure used to identify the GCO fragments maximizes the R01 term while, at the same time, minimizing the R01 term.[25] Formally, the molecular fragments A, B, and R are defined per normal mode. In practice, it is found, however, that for most modes, the atoms involved in polar bonds provide the most important contributions to the VCD intensities. For this reason, one can often group the atoms of a molecule into VCD-active and VCD-inactive sites. This simplifies the interpretation of the VCD spectra even further because according to eq , intense VCD signals arise when the VCD-active sites of a molecule are in a favorable orientation.

It is important to recognize that the standard CO model is a model used to predict VCD intensities, while the VCD GCO analysis is a postanalysis method performed after running a DFT VCD calculation that is aimed to gain physical insight into the origin of the computed VCD intensities. To connect the R01 term with the standard CO model, we define the CO correction vector Y⃗ aswhere Y⃗ is the distance vector determined by the mass centers of the A and B fragments. Combining eq with eq allows one to link the R01 term directly to the standard CO term R01 viawhere

The magnitude of the errors made when using the standard CO model can be assessed by comparing the values computed with DFT for R01, R01, R01, R01, and R01.

Magnitude and Sign of the VCD Intensities

As mentioned above, typically all intense bands that stand out in a VCD spectrum are characterized by large R01 values which are significantly larger than their associated R01 contributions. Small- and medium-sized VCD bands, on the other hand, are most often characterized by comparable R01 and R01 values. On the basis of eq , one could thus conclude that large VCD intensities will typically be observed when the EDTMs associated with fragments A and B are in a favorable orientation.

This interpretation of eq explains in very simple terms why VCD spectra depend so sensitively on the molecular conformation. Besides this, it also provides insight into the stability of the computed VCD sign, that is, the VCD robustness[26] of a mode. In a series of studies,[19,27−29] Magyarfalvi et al. proposed to use the magnitude of the ratio ζ(j) between the rotational and dipole strengths of the mode j—a quantity that is very similar to Kuhn’s dissymmetry factor—as a measure of its VCD robustness. Modes with ζ(j) > 10 ppm are considered to be robust and their computed VCD sign reliable, whereas modes with ζ(j) < 10 ppm should be considered as nonrobust and their computed sign unreliable.

Using eqs , 2, 5, and 9, the ζ(j) ratio is expressed as

When R01(j) is the dominant contribution to R01(j), eq can be approximated aswhere θ(j) is the angle between the E⃗01(j) and E⃗01(j) and φ(j), the angle between Y⃗ and E⃗01(j)×E⃗01(j). When sin(θ) or cos(φ) is close to zero, small changes in the relative orientation of the fragments A and B can easily induce a sign change and/or a large change in the magnitude of R01(j). Equation , therefore, suggests that when evaluating the robustness of a mode whose VCD intensity is dominantly determined by R01, it is not sufficient to consider only the magnitude of the associated ζ ratio. Instead, it is necessary to identify the GCO fragments that provide the main contribution to the rotational strength and to investigate whether at ambient temperatures the relative orientation of these two fragments can be changed easily. This is because the sign and magnitude of R01 can critically depend on small variations in the angles θ and φ.

Molecular Structure

A schematic representation of the ForValNHMe peptide is shown in Figure . Eight ForValNHMe conformers have been considered of which four are in the α-helix and four in the β-sheet configuration. Six of the conformers have been taken from ref (19). Following the labeling convention used in ref (19), these conformers are labeled as α1, α2, α3 and β1, β2, β3. For a given configuration (i.e., α or β), these conformers differ only in the orientation of the isopropyl group with respect to the backbone. The α conformers taken from ref (19) have been mirrored to assure that all considered conformations have the same absolute configuration.

Figure 1

Schematic representation of the ForValNHMe peptide as well as definition of the fragments used in the GCO analysis of the amide I modes. The CO fragments A and B are highlighted in red and blue, respectively, while the atoms in the fragment R are shown in black.

The linear transit (LT) scans that will be discussed below revealed two additional conformational energy minima, one in the vicinity of the conformer α1 and one in the vicinity of the conformer β1, which have, therefore, also been considered and labeled as α4 and β4 (see Figure S1 for the structures of all conformers). The conformers α2, α3, α4, and β4 exhibit an intramolecular hydrogen bond between the O1 atom and the N–H bond on the C-terminal, which plays a crucial role in the present study. To emphasize this aspect, these conformers will, therefore, be referred to as α2*, α3*, α4*, and β4*. The molecular structures of all conformers have been optimized at various levels of theory with structures determined at the BP86/TZP level of the theory reported in Table S1 of the Supporting Information.

Results and Discussion

Figure shows the amide I and amide II regions of the VCD spectra computed for the various ForValNHMe conformers. Interestingly, in the amide I region, some conformers exhibit the expected bisignate VCD pattern (α1, α4*, β2, and β3) while the other conformers do not show this feature (α2*, α3*, β1, and β4*). Moving to the amide II region, we find that the spectra of all β conformers exhibit very intense VCD bands with magnitudes comparable to, or even larger than, those of the amide I bands, while the α conformers display only very weak signals. To understand these differences, we will analyze in the following the amide I and amide II VCD bands using the GCO analysis introduced above. Apart from revealing the differences observed between the considered conformers, this analysis will also provide insight into the robustness of the amide I and amide II VCD bands. In the following section, we will discuss these results and refer to the rotational strengths and its various contributions in units of 10–44 esu2·cm2 while the EDTMs are referred to in units of 10–20 esu·cm.

Figure 2

Amide I and amide II regions of the VCD spectra computed for the eight ForValNHMe conformers considered in this study.

Amide I Modes

Normal Mode Analysis

The ForValNHMe peptide has two amide moieties (see Figure ). Normal mode analysis finds that in all conformers, the amide I modes (i.e., modes 60 and 61) dominantly consist of stretching motions of the two C=O bonds but also involve C–H and N–H bending motions. Except in the conformer α4* where the two carbonyl bonds move with similar amplitudes, we find that the normal mode motion is not localized equally on the two carbonyl bonds. One can thus distinguish between a dominant C=O bond that moves with large amplitude and a secondary C=O bond that oscillates with significantly smaller amplitude. Furthermore, modes 60 and 61 can be loosely classified—because there is no symmetry—as in- and out-of-phase with respect to the C=O stretch motions. As can be concluded from Table S9 in the Supporting Information, the out-of-phase mode has a lower frequency than the in-phase mode in conformers without an intramolecular hydrogen bond, whereas the reverse holds for the conformers with an intramolecular H bond. However, irrespective of the conformation, the C5=O6 stretch is always dominant in the mode with lower frequency and the C2=O1 stretch is dominant in the mode with higher frequency as exemplified in Figure where the nuclear displacement vectors of modes 60 and 61 of the conformer α3* are shown.

Figure 3

Analysis of the two amide I modes of the conformer α3*. Left: Nuclear displacement vectors. Right: EDTMs (E⃗01 and E⃗01) associated with the two carbonyl groups.

Generalized CO Analysis

A GCO analysis was performed by decomposing R01 into the contributions defined by eqs –8 with molecular fragments as defined in Figure . We first analyze the magnitude and the sign of the R01, R01, R01, and R01 contributions to the rotational strength R01 of modes 60 and 61. Table S9 shows that the magnitude of R01 varies significantly between the conformers considered and that it is not always the dominant contribution to R01. Similarly, we find that the relative magnitude and sign of R01 and R01 associated with a given normal mode depend sensitively on the molecular conformation. Regarding the magnitude of these two terms, R01 is found to be larger than R01 for all the conformers except for α4* and β4*. Further, because modes 60 and 61 are not pure carbonyl stretching modes, the R01 and R01 contributions to the rotational strength can have sizeable magnitudes and cannot always be neglected. Indeed, the analysis in Table S9 shows that these two terms are the reason why the VCD spectra of the conformers α2*, α3*, β1, and β4* do not exhibit a bisignate VCD pattern in the amide I region. To illustrate this, we note that the generalized CO term R01 (the sum of R01 and R01) always has opposite signs for the two amide I modes of a given conformer. VCD spectra that exhibit a bisignate feature are thus associated with conformers in which R01 is the dominant contribution to the rotational strength. In contrast, for conformers where the sum of R01 and R01 is comparable in magnitude with R01 and also has an opposite sign, VCD spectra may not exhibit a bisignate VCD feature, and this is exactly the case for the conformers α2*, α3*, β1, and β4*.

To understand why the R01 term depends so sensitively on the molecular conformation, we consider the quantities that determine its magnitude, that is, the magnitude of the EDTMs associated with the two fragments as well as the angle θ, the angle between the E⃗01 and E⃗01, and ϕ, the angle between Y⃗ and E⃗01(j) × E⃗01(j). We then find (see Table S10) that the EDTM associated with the dominant carbonyl bonds has a similar magnitude in all modes while the EDTM associated with the secondary carbonyl bonds exhibits much larger variations. This is a consequence of the fact that the involvement of the secondary C=O bonds in the mode motion varies from one conformer/mode to the other. The comparison of the |E|·|E|, sin(θ)cos(ϕ), and R01 values moreover shows that the large values observed for R01 in the conformers α1 and β2 are a consequence of the fact that in these conformers E⃗01 and E⃗01 are favorably oriented. Therefore, it should be clear that in the considered conformers, the relative orientation of E⃗01 and E⃗01 is more important than their magnitude. For example, the conformer α1 has the largest R01 values but these values are obtained by combining the smallest |E|·|E| values with the largest sin(θ)cos(ϕ) values. On the other hand, because α4* has the smallest sin(θ)cos(ϕ) value, the conformer α4* has the smallest R01 values even though it shows the largest |E|·|E| value.

The analysis performed in this section has clearly shown why the standard CO model cannot be used to interpret the VCD bands in the amide I region. Furthermore, the GCO terms have the potential to be the dominant contributions to the rotational strength and thus determine the observed VCD pattern. Whether or not that is the case depends, however, on the structural details associated with the various conformations, that is, on whether the E⃗01, E⃗01, and Y⃗ are in favorable orientations. In addition, we have shown that hydrogen bonding may also influence the observed VCD patterns as it affects the relative position of the in-phase and out-of-phase modes.

LT Calculations

To investigate the robustness of the VCD sign computed for modes 60 and 61, we performed LT calculations in which a sequence of related constrained geometry optimizations is performed along a predefined geometrical path. Such calculations are quite useful as they serve as the first approach to account for structural variations as occurring at room temperature instead of performing full molecular dynamics calculations. Furthermore, they enable one to assess how the various factors contributing to the rotational strength of a particular mode depend on the geometry. Both aspects are highly important for the present study as they allow us to determine and understand the robustness of the amide I modes. In this respect, the dihedral angle τ defined by the two carbonyl bonds is a parameter that can be expected to be particularly important.

Starting from the value in the fully optimized structure, τ was varied by taking 20 steps of 2.5° in each direction. For LT structures that fell within an energy window of 1 kcal/mol with respect to the fully relaxed structure, VCD calculations were performed allowing us to monitor how R01 and its contributions as well as ζ change during the LT scans. Although such an analysis was done for all eight conformers, we will in the following focus on the conformer α3* because the GCO analysis of its VCD spectrum is representative for the analysis of the spectra of the other conformers and provides insight into both the stability of the amide I bands and the GCO mechanism. The analyses performed for the other conformers are given in Section 3 of the Supporting Information. These analyses lead to similar observations and conclusions as drawn here for the conformer α3* and fully support the subsequent analysis.

The upper panel of Figure shows the variation of energy during the LT scan. As can be seen, the energy varies smoothly with energy variations of less than 0.1 kcal/mol for angles between 160 and 190°. This is important as it implies that the relative orientation of the two carbonyl bonds may easily be influenced by perturbations such as the interaction with a solvent but will also be critically dependent on the use of different computational parameters. The middle panel of Figure shows the variation of R01, R01, R01, and R01 + R01 associated with modes 60 and 61. We find that the overall change of R01 in the magnitude and the sign during the entire LT scan is quite significant: R01 changes from +96.5 to −303.2 for mode 60, and from −166.0 to +53.8 for mode 61. In contrast, the magnitude of the dipole strength D01 does not change appreciably during the LT scan (see Table S11). Therefore, the large variations observed for R01 are expected to be directly reflected in the ζ ratios as is confirmed by the middle and lower panels of Figure .

Figure 4

Analysis of the LT calculations performed for the α3* conformer. Upper panel: Variation of energy during the LT scan. Middle panel: Variation of the rotational strength R01 and its contributions during the LT scan. Lower panel: Variation of ζ during the LT scan.

These results thus indicate that in order to assess the robustness of the amide I modes, one should consider the profile of the potential energy surface and the position of the relaxed structure with respect to the points where R01 changes sign. For this purpose, we highlight in Figure the position of the relaxed structure (τ = 175°) in blue while the LT structures where R01 sign changes is highlighted in purple (τ = 170°) and teal (τ = 207.5°) for modes 60 and 61, respectively. As can be seen from the figure, the relaxed structure is located quite close to the point where R01 changes sign for mode 60 but far from the point where R01 changes sign for mode 61. Besides explaining the small VCD intensity in mode 60 (−8.5) and the large VCD intensity in mode 61 (−92.6), this analysis also provides insight into the robustness of the two VCD bands. To change the sign of mode 60, it is sufficient to change τ from 175 to 170°. Such a change is accompanied by a very small energy change (0.003 kcal/mol) and the sign computed for this mode is thus likely to depend on the choice of computational parameters. Mode 61, on the other hand, should be assigned as a robust mode because for a sign change to occur, τ should be changed by more than 32.5°. Even though this change is accompanied by a relatively small energy change (0.25 kcal/mol), it is unlikely that the use of different basis sets or functionals will yield such a large variation in τ.

The VCD robustness criterion proposed by Magyarfalvi and Gobi, on the other hand, would lead to the conclusion that both modes should be considered as nonrobust as ζ(60) = −2.6 and ζ(61) = −9.7. To test the validity of our predictions, we have performed VCD calculations for α3* using different computational parameters, that is, vacuum/BP86/OLYP/TZP/TZTP/DZP, COSMO(DMSO)/BP86/TZP, and COSMO(chloroform)/BP86/TZP. These calculations show that the sign changes only occur for mode 60, and even then only between vacuum and COSMO calculations, for example, from −8.5 in the vacuum/BP86/TZP calculations to +23.5 in COSMO(chloroform) and +33.0 in COSMO(DMSO) calculations. These results thus support the robustness assignment based on the analysis shown in Figure . Similar conclusions are reached when the other conformers are considered (see Section S3 of the Supporting Information).

Variation of R01 and Its Contributions during LT Scan

To further elucidate the dependence of the amide I VCD spectrum on the molecular structure, it is of interest to consider how R01 varies during the LT scan from the perspective of the CO model (see Figure ). Because the standard CO term R01 and R01 exhibit similar trends, it is clear that R01 is the driving force behind the changes observed in R01. At the same time, it is important to note that R01 and R01 often have different magnitudes and can also have different signs. This occurs because R01, R01, and R01 are not negligible. As shown in Figure , the two amide I modes involve C=O stretches along with N–H and C–H bending motions. Because the latter contributions are not the same in the two modes, R01 changes sign at different τ values for the two amide I modes, viz. R01 + R01 has a very different value for both modes. This explains why the standard CO model is not able to predict accurately the sign and magnitude of R01.

The R01 term exhibits a behavior during the LT scan that is typical for CO terms. First, it is observed that when the two carbonyl bonds are in a favorable orientation, the magnitude of R01 is very large. As the relative orientation of the C=O bonds changes, R01 becomes gradually smaller and eventually changes sign. After the sign change, the magnitude of R01 increases steadily as the relative orientation of the two bonds improves. Second, for a given LT structure, R01 has different signs in the in-phase and out-of-phase modes. The only exception is found for the LT structure with τ = 192.5°, which is very close to the point where R01 changes sign. This small discrepancy is due to the perturbed in-phase/out-of-phase behavior resulting from the involvement of the C–H and N–H bending motions. It is also important to notice that R01 does not change sign when τ is 180° as the standard CO model would predict. This is a consequence of the fact that the EDTM associated with a C=O bond is not oriented exactly along the bond. As a result, the angle θ between E⃗01 and E⃗01 can be very different from the dihedral angle τ determined by the two C=O bonds. Indeed, as illustrated in Figure , using the VCDtools GUI, the E⃗01 and E⃗01 vectors deviate from the direction of the bonds by roughly 20°.

To explain the sign change observed for R01 during the LT scan, the variation of the angles θ (the angle between E⃗01 and E⃗01) and ϕ (the angle between Y⃗ and E⃗01 × E⃗01) during the LT scan has been plotted and is shown in Figure . Surprisingly, for both amide I modes, θ varies only minimally during the LT scan: variations of 90° in τ lead to variations in θ of less than 10°. The sign changes observed in Figure for R01 are thus not caused by θ but by ϕ. Indeed, ϕ exhibits significant changes during the LT scan and crosses the 90° line at approximately 192.5 and 197.5° for modes 60 and 61, respectively. These angles coincide exactly with the angles τ at which R01 changes sign for both modes.

Figure 5

Dependence of the angles θ and ϕ associated with normal modes 60 and 61 of conformer α3* on the dihedral angle between the carbonyl bonds.

The analysis in Figure also provides further insight into the CO correction term R01. During the LT scan, the R01 term exhibits a typical CO behavior as its magnitude varies gradually, eventually changing sign and also has different signs for modes 60 and 61. Interestingly, it is found that R01 and R01 change sign at different geometries. The latter observation has two important implications. First, the relative magnitudes and signs of R01 and R01 depend on the relative orientation of the C=O bonds. This means that in certain τ ranges, R01 and R01 not only have opposite signs, but also the correction term is larger than the standard term. This situation occurs for mode 60, when τ is in the (187.5, 192.5°) interval and for mode 61, when τ is in the (182.5, 190°) interval. Second, the Y⃗ and Y⃗ vectors are not parallel. In refs (21) and (22), these vectors were considered to be parallel when quantifying the errors made in the standard CO model. While from a mathematical point of view this assumption is correct, the observation made above shows that in reality the two vectors do not need to be parallel. To this end, we note that—because the vector Y⃗ associated with a given mode is not uniquely defined[22]—it is better to consider its projection on the associated E⃗01 × E⃗01 vector, which is fixed according to eq . We have, therefore, monitored in Figure how the magnitude of the E⃗ × E⃗ vector (|E⃗ × E⃗|), and the projections of the associated Y⃗, Y⃗, and Y⃗ vectors on E⃗ × E⃗ (|Y⃗∥|, |Y⃗∥|, and |Y⃗∥|, respectively) vary. Figure shows that for both modes, |E⃗01 × E⃗01| shows only small variations during the LT scans. On the other hand, |Y⃗∥|, |Y⃗∥|, and their sum |Y⃗∥| vary gradually during the LT scan and eventually all change sign. This behavior, which mirrors closely the variations observed in Figure for the angles θ and ϕ, further supports the conclusion that R01 behaves like a typical CO term.

Figure 6

LT dependence of the magnitude of E⃗ × E⃗ for modes 60 and 61 of conformer α3 as well of the projections (∥) of the associate Y⃗, Y⃗, and Y⃗ vectors on E⃗ × E⃗.

The above analysis has clearly shown that it is far from trivial to explain the intensities and signs of the amide I VCD bands in peptides. It is only by careful and detailed calculations that we can get a hold on the origin of the intensities and the robustness of the bands. GCO analyses in which the geometry dependence and relative magnitudes of the R01, R01, R01, R01, and R01 contributes to R01 are unraveled and are in this respect indispensable.

Amide II Modes

As discussed, the comparison of the computed VCD spectra in Figure shows that in contrast to the α conformers, the β conformers display very intense bands in the amide II region. In this section, we aim to explain these differences.

The conformers of the ForValNHMe peptide have eight normal modes in the 1430–1580 cm–1 amide II spectral region. A detailed normal mode analysis (vide infra, see Section S4 in the Supporting Information) reveals that these modes can be grouped into three different categories (A, B, and C). Type A modes are standard amide II modes in which motion is mostly localized on the C–N–H groups and in which coupling between the motions of C–N–H groups of the C- and N-termini is typically very small. The frequencies of these modes show a clear difference between conformers with or without intramolecular H-bonding, the former being blue shifted with respect to the latter. This blue shift is well-known and finds its origin in the stabilization of the double-bond character of the C–N bond.[30] Type B modes are mixed modes not only consisting of C–H bending but also involving motion of the C–N–H groups. Because the magnitude of the atomic EDTMs associated with the C–N–H groups are comparable or even larger than those of the moving C–H bonds, B modes often have a strong amide II character. Two factors appear to determine the number of type B modes in a given conformer. The first one is the configuration of the conformer as β-sheet conformers have more type B modes than α-helix conformers. The second one is the presence or absence of intramolecular H-bonding: conformers without intramolecular H-bonding are expected to have more B-type modes because in these conformers, the amide II modes (the A-type modes) have lower frequencies and thereby mix more with the C–H bending modes. Type C modes are typical C–H bending modes that do not mix appreciably with the amide II modes. For the present study, these modes are less important as they generally only have low VCD intensities. As summarized in Table the 64 modes of the eight conformers in the amide-II spectral region can be grouped into 15 A-type, 19 B-type, and 30 C-type modes.

Table 1

Summary of the Robustness Analysis Performed for the Modes in the Amide II Spectral Region: the Number of Normal Modes (No. NMs) and the Number of VCD Sign Changes (No. Sign Changes) Observed When Comparing Calculations Performed with Six Different Computational Parameters (BP86/DZP, BP86/TZP, BP86/TZ2P, OLYP/DZP, OLYP/TZP, and OLYP/TZ2P)

	A-type modes		B-type modes		C-type modes		all modes
conformer	no. NMs	no. sign changes	no. NMs	no. sign changes	no. NMs	no. sign changes	no. NMs	no. sign changes
α2*	2	1	0	0	6	1	8	2
α3*	2	0	0	0	6	5	8	5
α4*	2	1	1	0	5	1	8	2
β4*	1	1	5	1	2	1	8	3
total H bonding	7	3	6	1	19	8	32	12
α1	2	0	4	2	2	0	8	2
β1	2	1	3	2	3	1	8	4
β2	2	0	4	1	2	2	8	3
β3	2	0	2	0	4	3	8	3
total rest	8	1	13	5	11	6	32	12
total all modes	15	4	19	6	30	14	64	24
total all α modes	8	2	5	2	19	7	32	11
total all β modes	7	2	14	4	11	7	32	13

GCO Analysis

GCO analyses of the amide II modes of the eight conformers show that there are three factors that determine the magnitude of the VCD signals of these modes. The first one is the type of normal motion. Localized normal modes cannot trigger the GCO mechanism and as such most often have small or medium intensities. Mode 59 of α3* is in this regard a typical example. Figure shows that it is localized entirely on the N-terminal and has a medium-sized rotational strength even though it has a very large total EDTM. However, because the rest of the atoms in the molecule are not involved in the normal mode motion, the EDTM associated with the rest of the molecule is for all practical purposes considered to be zero. As a result, the large EDTM associated with the C–N–H atoms does not have the possibility to couple with the EDTM of another moiety (see GCO analysis in Table ). Delocalized modes, on the other hand, have the potential to trigger the GCO mechanism and to have very large VCD intensities.

Figure 7

Nuclear displacement vectors of the normal modes in the amide II spectral region. The left panel shows the modes of conformer α3* while the right panel shows the modes of conformer β2. To illustrate the amide II character of each mode the radius of the atoms has been scaled with the magnitude of the associated atomic EDTMs. The fragments used in the GCO analysis performed in Table for modes 54(α3*), 59(α3*), 53(β2), and 55(β2) are highlighted in red and blue.

Table 2

GCO Analysis Results of Amide II Modes of Conformers α3* and β2a

conf.	NM	freq.	D01	R01	R01GCO	R01IF	R01R
β2	53	1445.6	416.7	–142.4	–145.8	+16.1	–12.6
β2	55	1454.2	279.6	+230.9	+185.4	+49.1	–3.6
α3*	54	1449.3	4.3	+3.7	+0.2	+4.0	–0.5
α3*	59	1565.6	387.8	–24.5	+2.8	–21.5	–5.8

GCO fragments are defined in Figure . Frequencies are given in cm–1, dipole strength D01 in 10–40 esu2·cm2, and rotational strength R01 and its contributions R01, R01, and R01 in 10–44 esu2·cm2.

The second factor is the involvement of the C–N–H group in the normal mode motion. Typical amide II modes involve stretches of the C–N bond along with bending of the N–H bond. Because of the polar character of these two types of bonds, the associated atomic EDTMs are significantly larger than those associated with the C–H bending motions. As indicated previously, even small motions of the atoms in the C–N–H groups induce atomic EDTMs that are comparable to—or larger than—those associated with C–H bending motions with large amplitudes. Because large R01 values require at least one of the interacting EDTMs to be large, intense VCD bands are thus associated with normal modes that are localized on at least one of the C–N–H groups. This point is further illustrated by the GCO analyses reported in Table . Modes 53 and 55 of conformer β2 exhibit very intense VCD signals. In line with the previous reasoning, we find that in both cases, R01 is significantly larger than the other two contributions to the rotational strengths. We can, therefore, conclude that the large VCD intensities are the result of the GCO coupling between the EDTMs associated with the C–N–H groups. In mode 54 of the conformer α3*, on the other hand, we have only C–H bending movements which generate weak EDTMs that are unable to couple. As a result, the mode is accompanied by only a small R01 value and a small total VCD intensity.

The final factor concerns the relative orientation of the interacting EDTMs. This is a consequence of eqs , 12 and 13. As can be seen in Figure , not all modes involving movements of the C–N–H group along with C–H bending motions have large rotational strengths. Instead, large rotational strengths are only observed for modes where the interacting EDTMs have a favorable orientation with respect to each other.

The normal mode and GCO analyses discussed above clearly explain why the β conformers exhibit very intense bands in the amide II region while α conformers do not. The intense VCD signals observed in this region are determined by large GCO terms which in turn are determined by the through-space interaction of large EDTMs. As a result, intense VCD signals are more likely to be observed in mixed modes that are delocalized and characterized by large EDTMs, that is, in B-type modes. The A-type modes—which are localized—and the C-type modes—which are delocalized but have small EDTMs—most often have small VCD intensities.

VCD Robustness in the Amide II Region

Unlike the amide I modes, the stability of most of the intense VCD bands in the amide II spectral region cannot be investigated consistently by LT calculations because the delocalized nature of the GCO fragments makes it difficult to find LT parameters that allow one to change the computed VCD sign without changing the conformation. For this reason, the robustness of the computed VCD sign was investigated by performing geometry optimizations and VCD calculations with six different computational parameters (BP86/DZP, BP86/TZP, BP86/TZ2P, OLYP/DZP, OLYP/TZP, and OLYP/TZ2P). To identify VCD sign changes induced by the use of different computational parameters, the normal modes predicted by the six sets of calculations have been mapped by computing normal mode overlaps.[26]

The results of these analyses are summarized in Table , which shows that 24 of the 64 modes in the amide II spectral region change sign upon changing computational parameters. Interestingly, no significant differences are observed between α-helix and β-sheet conformers, nor between structures with and without intramolecular hydrogen bonding. When grouping the modes into A, B, and C types, VCD sign changes are observed for 4 out of the 15 A-type modes (27%), 7 of the 19 B-type modes (37%), and 14 of the 30 C-type modes (47%). It is important to realize that the number of B and C modes found in a given conformer depends on the computational parameters used in the calculations, in contrast to the number of A modes, which is not affected by the change of computational parameters.

Type B modes are predominantly C–H bending modes that also involve small movements of the C–N–H group(s). This weak coupling between the C–H bending motions and the movements localized on the C–N–H group is very sensitive to the computational parameters, and this has far-reaching consequences. To illustrate this point, we consider mode 56 of the conformer β3 whose nuclear displacement vectors calculated at the BP86/DZP and BP86/TZP levels are compared in Figure . As can be seen, both calculations predict the mode to be a C–H bending mode localized on the C-terminal methyl group and involving also small N–H bending contributions. For the H atoms moving with large amplitudes, the two sets of calculations make similar predictions for the nuclear displacement vectors. In contrast, the two sets of calculations make quite different predictions for the small-amplitude C–H and N–H bending motions. Not only is localization of the motion on the N–H bond different (4.4% in BP86/TZP vs 11.1% in BP86/DZP), but the two calculations also predict a different phase for the coupling between the C–H and N–H bending motions (see Figure ). The consequence of this in-phase/out-of-phase behavior is quite dramatic as it leads to opposite VCD signs for the two modes.

Figure 8

Normal mode 56 of ForValNHMe-β3 calculated at the BP86/TZP (top) and BP86/DZP (bottom) levels. The green arrows indicate the direction and amplitude of the vibration. The size of the atoms corresponds to the magnitude of the associated atomic EDTM.

From the perspective of the GCO analysis, the two modes in Figure are quite different. Both calculations predict that the R01 term is the largest contribution to the rotational strength, but its sign and magnitude are different in the two calculations (see Table ). The difference in the sign is a consequence of the different phase observed in the mode coupling, while the difference in magnitude is caused by the different involvement of the N–H bending motion in the normal mode motion. Further, we note that both sets of calculations predict large absolute ζ values which provides further evidence that the magnitude of ζ does not reflect the stability of the computed VCD sign when the VCD intensity of a mode is determined by the GCO mechanism.

Table 3

Normal Mode and GCO Analyses of Mode 56 of Conformer β3 Predicted at the BP86/TZP and BP86/DZP Levelsa

calc.	conf.	mode	overlap	freq.	ζ	D01	R01	R01GCO	R01IF	R01R
BP86/TZP	β3	56	1.00	1457.5	–117.7	7.2	–8.5	–5.33	–4.42	+1.24
BP86/DZP	β3	56	0.87	1457.3	+44.6	93.3	+41.6	+41.58	–0.47	+0.49

C-terminal methyl group was chosen as GCO fragment A and the N-terminal C–N–H group as GCO fragment B (see Figure ). Frequency is given in cm–1, dipole strength D01 in 10–40 esu2·cm2, rotational strength R01 and its contributions R01, R01, and R01 in 10–44 esu2·cm2, and the dissymmetry factor ζ in ppm. The normal mode overlap is dimensionless.

Mode 56 of the conformer β3 is representative for most modes in the amide II spectral region of all conformers. Analysis of the data shows that 15 of the 24 observed VCD sign changes are caused by computational uncertainties related to the localization of the normal mode motion on the C–N–H group while the remaining sign changes are associated with modes that have very weak VCD intensities (e.g. |R01| < 10). Similarly, there is hardly any correlation between the magnitude of ζ and the observed VCD sign changes which occur for both small and high ζ values. In fact, in some conformers (α3*, β1, β3, and β4*) the modes with a high ζ value have a nonrobust VCD sign. From the perspective of the GCO mechanism, this is not surprising. Large ζ values are often correlated with large R01 contributions but R01 depends sensitively on the computational parameters used in calculations as has been shown above.

Finally, we note that the present study has also shown that differences between the relaxed structures predicted by various computational parameters may induce VCD sign changes. For example, the optimized structures with the DZP basis set do not contain the H-bond in the conformers α2* and α3*. As a result, some of the normal modes show VCD bands with opposite signs. To avoid such errors, we have employed, in the present study, the relaxed TZP structures as the starting point for the geometry optimizations performed using DZP. With this approach, we find very similar relaxed structures predicted by the six sets of calculations.

Conclusions

Using the ForValNHMe peptide as a prototypical peptide, we have presented in this work a thorough analysis of the amide I and amide II VCD bands of peptides. By combining the recently developed VCD GCO formalism with a novel normal modes analysis that monitors the magnitude of the nuclear displacement together with the magnitude of the vibrational DTMs, we have been able to explain all the differences observed for the α and β conformers in these spectral regions. The performed analyses have revealed which electronic and spatial parameters determine the intensity of the amide I and amide II VCD bands, and have highlighted the critical role played by intramolecular hydrogen-bonding. As shown, the GCO mechanism is responsible for all VCD bands that stand out in both spectral regions. Essentially, this means that the same underlying mechanism, that is, the through-space coupling of electric dipoles, is responsible for the intense VCD bands observed in the amide I and amide II spectral regions. As a result, both amide I and amide II VCD bands are found to be very sensitive to the relative orientation of the C- and N-termini and on the localization of the normal mode motion on these two termini. This, in turn, implies that one should be cautious to determine the stability of computed VCD signs by evaluating the magnitude of the ζ ratio, as proposed previously.

The present study has also provided additional insight into the VCD GCO mechanism. By studying the trends of the R01 and R01 terms during LT scans, new insight was gained regarding the CO correction term R01. It was shown that R01 and the standard R01 CO term behave in a very similar fashion during the LT scans. However, the behavior of the R01 term is shifted with respect to R01. Therefore, R01 can be larger than R01 while at the same time have a different sign. In other words, terms that are neglected in the standard CO model can not only be larger than what this model considers to be important but also have the opposite sign.

The detailed insight that has been obtained in the VCD properties of the ForValNHMe peptide is essential to further our understanding of the VCD spectra of peptides in general, and how these spectra reflect key details of their electronic and spatial structure. At an even more general level, it is clear that supported by theoretical GCO analyses VCD spectroscopy is coming of age and becoming an indispensable tool for studies of the stereochemistry of chiral systems.

Methods

All geometry optimizations,[31] VCD calculations,[32] and LT scans were performed at the DFT level of theory using the ADF program package.[23,33] Two exchange–correlation functionals (BP86[34,35] and OLYP[36,37]) and three different basis sets (the ADF DZP, TZP, and TZ2P basis sets[38]) were used to investigate the effects of different parameters. For the larger part of the calculations, BP86/TZP was used, and these as such served as a reference when comparing between different sets of parameters. In order to investigate solvent effects, the COSMO-solvation model was used as well.[39,40]

To determine the dependence of the calculated VCD sign on the relative orientation of the coupled molecular fragments, LT calculations[41,42] were performed. In these calculations, the molecule structure was changed systematically starting from the fully relaxed structure by varying a LT parameter (e.g., bond length, angle, or dihedral angle) in predefined steps. After each LT step, a constrained geometry optimization was performed by relaxing all structural parameters except the LT parameter which was kept at the predefined value.

The VCD spectra, normal modes, and EDTMs were analyzed and visualized using a new implementation of the VCDtools program[43] in the ADF GUI program.[25] VCDtools computes the rotational strength R01 and its various contributions, for example R, R, R, R, and R using atomic polar tensors, atomic axial tensors, and nuclear displacement vectors obtained from DFT calculations.