Kinetically Trapped Liquid-State Conformers of a Sodiated Model Peptide Observed in the Gas Phase

We investigate the peptide AcPheAla5LysH+, a model system for studying helix formation in the gas phase, in order to fully understand the forces that stabilize the helical structure. In particular, we address the question of whether the local fixation of the positive charge at the peptide's C-terminus is a prerequisite for forming helices by replacing the protonated C-terminal Lys residue by Ala and a sodium cation. The combination of gas-phase vibrational spectroscopy of cryogenically cooled ions with molecular simulations based on density-functional theory (DFT) allows for detailed structure elucidation. For sodiated AcPheAla6, we find globular rather than helical structures, as the mobile positive charge strongly interacts with the peptide backbone and disrupts secondary structure formation. Interestingly, the global minimum structure from simulation is not present in the experiment. We interpret that this is due to high barriers involved in re-arranging the peptide-cation interaction that ultimately result in kinetically trapped structures being observed in the experiment.


Introduction
Helical secondary structural motifs, such as α and 3 10 , are common in proteins. 1 In solution, helix propensity is determined both by intramolecular interactions and protein-solvent interaction. Gas-phase systems offer the opportunity to study the "undamped" intramolecular interactions that shape peptides, thereby shedding light on intrinsic helix propensities and bonding interactions. Gas-phase helices have been investigated using ion mobility spectrometry 2-4 and vibrational spectroscopy. [5][6][7][8][9][10][11][12] The combination of these experimental techniques with molecular simulations based on density-functional theory (DFT) allows for structure elucidation, as it helps to interpret experimentally obtained spectra. Moreover, a rigorous experiment-theory comparison allows for the assessment of the accuracy and predictive power of simulation approaches. 13 Pioneering ion-mobility experiments in the group of Jarrold 2,3 examined the role of N-and C-terminal residues on gas-phase helix formation for the sequences Ala n H + , AcLysAla n H + , and AcAla n LysH + . They concluded that Ala n H + and AcLysAla n H + adopt globular conformations in the gas phase independent of the length of the amino-acid chain while AcAla n LysH + is helical for n > 8. 14 The identities of these structures were confirmed by theoretical and experimental vibrational spectroscopy in the work of Rossi et al. 10 and Schubert et al., 12 respectively. Similar studies focused on peptides of the form AcPheAla n LysH + with n = 1-5, 10, where phenylalanine provides a UV chromophore, which allows for conformer-specific IR-UV double resonance spectroscopy. [6][7][8][9] Vibrational signatures of individual conformers add a new dimension to peptide structural analysis beyond the orientationally averaged collisional cross section provided by ion mobility. In these experiments, the number of residues necessary to form a helix was found to be six, 8,14 but much of the hydrogen bonding pattern responsible for the formation of this motif is already present even with only three residues. 9,15 In conjunction with computational vibrational spectroscopy based on DFT, [8][9][10]16 such spectra allowed for determining detailed molecular structures and critically examining evidence for helix formation of peptides in isolation.  Figure 1: Illustration of helix-stabilizing factors for peptides in the gas phase (a) and structural formulas of (b) AcPheAla 5 LysH + and (c) AcPheAla 6 + Na + . Figure 1(a) illustrates the helix-stabilizing factors in polyalanine peptides, shown for the specific case of AcPheAla 5 LysH + . Work by the groups of Jarrold, 2,3 Rizzo, [6][7][8][9] and Blum 12 showed that intramolecular hydrogen bonds play an important role and that the design concept can even be transferred to non-natural peptides. 11 Hoffmann et al. 17 could show that deleting a single hydrogen bond had little impact on the overall helix stability. In addition to their energetic stability, hydrogen bonds are aligned in helices, and the resulting macro-dipole favorably interacts with the positive charge of the protonated lysine (Lys) sidechain at the C-terminus. Moreover, the capping of the "dangling" carbonyl groups near the C-terminus by the Lys side-chain provides additional stability.
To investigate the importance of the charge fixed at the C-terminus, we focus on the wellstudied system 8,16 of AcPheAla 5 LysH + and compare it to AcPheAla 6 + Na + (Figures 1(b) and 1(c), respectively). In the latter, Lys is formally replaced by alanine (Ala) and a sodium cation (Na + ) in order to introduce a freely movable positive charge. The resulting rich possibilities for electrostatic interaction can locally disrupt hydrogen-bonding networks and induce unconventional backbone conformations. [18][19][20][21] Consequently, the cation-binding site, and hence the conformation as a whole, is not a priori obvious. Ion mobility studies on metallated peptides (e.g. sodiated species of Ala n + M + 22 ) suggest that the cation plays the same role as the charged Lys side-chain in AcAla n LysH + for peptides with n > 12. For shorter peptides, calculated collisional cross sections (CCS) for globular and helical structures are both in agreement with the experimental CCS, preventing a definitive structural assignment.
In the present work, we couple IR-UV double resonance spectroscopy and theory in order to unravel the structure of the system of AcPheAla 6 + Na + with the aim of understanding whether a freely movable cation is sufficient to stabilize helix formation or if the C-terminal localization is a prerequisite for that.

Experimental Setup
The experimental setup has been described in detail elsewhere. 23 In brief, a nano-electrospray ion source is combined with a cooled ion trap (4 K) for spectroscopic studies of gas-phase ions. Conformer-selective IR spectra are recorded by applying IR-UV double resonance.
A measurement is performed by fixing the wavenumber of the UV laser to a line in the electronic spectrum and scanning the wavenumber of an infrared laser. When the IR pulse is in resonance with a vibrational transition of the ion, part of the population is removed from the ground state, leading to a decrease in UV-induced fragmentation. Scanning the IR wavenumber, one obtains a conformer-specific vibrational spectrum. Performing the same experiment on each line of the electronic spectrum allows for assignment of each UV spectral feature to a particular conformer.

Computational Methods
The applied conformational search algorithm is similar to the one used by Rossi et al. 16 First, a global conformational search is performed on the force field (FF) level using CHARMM22 24 and OPLS-AA, 25,26 separately. To that end, a basin-hopping approach 27 was applied using the scan program of the TINKER molecular modeling package. 28,29 For the system of AcPheAla 5 LysH + (AcPheAla 6 +Na + ) 603 280 (626 829) conformers were found using CHARMM22 and 643 938 (635 120) conformers were found using OPLS-AA. Single-point energy calculations at the generalized-gradient approximated (GGA) DFT level of theory have been performed for all these FF conformers. All DFT calculations were done using the all-electron/fullpotential electronic structure code package FHI-aims. [30][31][32] To be more precise, energies were computed at the PBE+vdW level, i.e. using the PBE 33 functional and a pair-wise van der Waals correction scheme (vdW). 34 Furthermore, FHI-aims-specific tier 1 basis sets and light settings have been used that are provided out-of-the-box to control the computational accuracy intended to give reliable energies energy for screening purposes. 30 For the two FFs individually, the 500 conformers with the lowest FF energy and the 500 conformers with the lowest DFT energy, i.e. a grand total of 2000 conformers, have been selected. The 2000 selected conformers were then geometry optimized at the PBE+vdW level using tier 1 basis sets and light settings. A hierarchical clustering scheme was applied in order to rule out duplicates. Further relaxation was then accomplished at the PBE+vdW level using FHI-aims-specific tier 2 basis sets and tight settings that are intended to provide meVlevel accurate energy differences, 30 i.e. within 0.02 kcal/mol. After clustering, this resulted in 324 (159) conformers for the system of AcPheAla 5 LysH + (AcPheAla 6 + Na + ) in the lowenergy region, i.e. within 6 kcal/mol from the global minimum. These conformers were then again locally refined at the PBE0+MBD level, i.e. using the hybrid exchange-correlation (xc) functional PBE0 35 augmented by a many-body dispersion (MBD) correction, 36 using tier 2 basis sets and tight settings which resulted in 52 (23) conformers in the low-energy region, i.e. within 3 kcal/mol from the global minimum.

Results and Discussion
AcPheAla 5 LysH + For our comparative study, a firm assignment of measured conformer-selective IR spectra to their calculated counterparts is of paramount importance. To that end, we first re-assess the peptide AcPheAla 5 LysH + and demonstrate that the applied conformational search technique completely grasps the conformational space energetically close to the global minimum, and that the applied level of theory is capable of reproducing the energetics as well as the vibrational properties of the conformers. For this we compare our results to previous work on AcPheAla 5 LysH + by Stearns et al., 8 where the 45 lowest-energy structures were selected out of a set of 1,000 force-field minima and subsequently optimized using DFT with a hybrid exchange-correlation (xc) functional. Even though four structures were successfully assigned to the experimental spectra, the question whether the search was complete and the whether these conformers are located in the global minimum region remained open. This did, in part, motivate an exhaustive conformational search by Rossi et al., 16 in which 7 conformers were found within 1 kcal/mol of the global minimum on the potential-energy surface (PES).
The authors were able to assign the experimentally observed structures to the global minima populated at low temperature by using the hybrid xc-functional PBE0, 35 augmented by a many-body dispersion (MBD) correction, 36 and including zero-point energy corrections.
The latter were computed with the generalized-gradient approximation functional PBE 33 and a pair-wise van der Waals correction (vdW), 34 which proved however unsatisfying for the prediction of vibrational spectra. It was suggested that using a hybrid xc-functional was necessary, which was a natural assumption since this level of theory was necessary for a correct conformational energy prediction in the first place. Furthermore, it was assumed that an anharmonic treatment was needed to yield improved spectra.
The conformational search strategy has already been laid out in detail in the previous section, including numbers illustrating the exhaustiveness of the search. The fact that we find two additional conformers within 1 kcal/mol from the lowest-energy conformer gives us confidence in the conformational search. The corresponding hierarchy of the relative DFT energy ∆E on the PES is shown in Figure 2. Nine conformers were found within 1 kcal/mol from the global minimum.
Since the experimental measurement takes place on cold ions in the gas phase, the PES High computational costs prohibited the systematic use of hybrid xc-functionals for the  calculation of harmonic vibrations in the previous study by Rossi et al. 16 To complete the picture, we repeat the harmonic vibrational free energy calculations at the PBE0+MBD level, confirming the already obtained result. Figure 3(a) shows the energy hierarchies for ∆E, ∆F (10 K), and ∆F (300 K) for the four lowest-energy conformers illustrated in Figure 3(b).
Conformers A and B are virtually identical near the C-terminus, but differ near the Nterminus by a tilted Phe side chain. The difference between conformers C and D is similar.
All four conformers show helical structure motifs: conformer C possesses one 3 10 -and two α-helical turns, conformer D features one 3 10 -and one α-helical turn, and conformers A and B each possess two 3 10 -and one α turn.
For this work, the original IR-UV double resonance experiment by Stearns et al. 8 has been repeated to allow conformer-selective IR spectra to be compared to their theoretical counterparts calculated at the PBE0+MBD level. The affiliated UV spectrum including peak assignments to their corresponding conformers is provided in the supporting information.
The conformer-selective IR spectra are shown in Figure 3(c). Conformers A and B could be attributed to their corresponding observed IR spectra. While the agreement is very good, the match between experimental and theoretical IR spectra is not perfect. This discrepancy is commonly attributed to two factors: (i) The effect of a possible incomplete characterization of electron exchange and correlation, despite the use of the hybrid functional PBE0, and (ii) the treatment of anharmonic vibrations and nuclear quantum effects. 37 Both of these effects are corrected for solely by applying a scale factor to the vibrational frequencies.
The assumption of a uniform overestimation of the harmonic vibrational modes with respect to experiment is debatable as they depend on the theoretical method, the used basis set, and the system itself. 38,39 In this work, we focus on the frequency region of 3200 cm −1 to 3500 cm −1 which is sensitive to N−H· · · O hydrogen bonding, where a uniform scaling factor of 0.948 yields very good agreement.
The exhaustive conformational search presented here for AcPheAla 5 LysH + , and the rigorous treatment of harmonic vibrations at the hybrid xc level allowed for (i) reproducing the known energy hierarchy and finding additional conformers in the low-energy region and (ii) calculating well-fitting harmonic IR spectra for the conformers in the low-energy region.
In this way we confirm the conformers predicted by Stearns et al. 8 and Rossi et al. 16 and can rule out any other competing conformers. This also shows that calculating computationally costly anharmonic IR spectra is not required in this case. Now that we have confirmed the accuracy of our simulation approach, we tackle AcPheAla 6 + Na + , a more challenging system because of the additional conformational degrees of freedom due to the "unfixed" cation.
AcPheAla 6 + Na + Figure 4 shows the energy hierarchies of the relative PBE0+MBD energies ∆E as well as the relative Helmholtz free energies ∆F at 10 K and 300 K with harmonic vibrational free energy contributions at the PBE+vdW level that were obtained for AcPheAla 6 +Na + . The four presumably dominant conformers are presented in Figure 5(b). Of the four conformer-selective IR spectra that were recorded, two of them correspond to conformers with particularly high intensity in the UV spectrum (see Figure S2, supporting information). The measured IR spectra of these two conformers, IIa and IIb, show very good agreement with the IR spec-  Relative energy [kcal/mol] Figure 4: Energy hierarchies of conformers of AcPheAla 6 + Na + at the PBE0+MBD energy ∆E as well as the relative Helmholtz free energy ∆F at 10 K and 300 K with harmonic vibrational free energy contributions calculated at the PBE+vdW level.
An obvious observation is the outstanding global minimum (conformer I in Figure 5 Ⅰ Ⅲ Ⅱa Ⅱb Exp. Exp.
Ⅰ Ⅱa Ⅲ Ⅱb Ⅰ Ⅱa Ⅲ Ⅱb  structure I does not seem to be observed in the experiment -none of the conformer-selective spectra fit the calculated vibrational signatures (see Figure 5(c)). The structure representing the global minimum is globular and features a cation-π interaction between the Na + and the Phe side chain. If that conformer were present in experiment, one would expect broad features in the UV spectrum due to charge-transfer between Na + and the aromatic ring.
However, no such features have been observed. The reason behind the kinetic trapping of conformers IIa and IIb has to be sought in the experimental procedure in which the molecules are electrosprayed into the apparatus from a solution at room temperature while the actual measurements are taken on isolated molecules at 10 K.
The energy landscape of the system may differ significantly between the system in solution at room temperature and in the gas phase at 10 K. This is particularly true for the global minima. While the global minimum in the gas phase may be energetically favored in comparison to the other structures, kinetic constrains, i.e. high energy barriers between minima, may hinder proper folding while transitioning from solution to the gas phase. Thus, experimentally observed local minima in the gas phase higher in energy are yielded due to their structural bias from aqueous solution at room temperature, resulting in kinetically trapped structures unable to transition into the global minimum.
It is obvious from comparing the ∆F (10 K) and ∆F (300 K) hierarchies (see Figure 5(a)) that the temperature difference does not contribute to a possible kinetic trapping effect.
In fact, the energy gap between the global and the next minimum even increases from 1.6 kcal/mol at 10 K to 2.0 kcal/mol at 300 K. Therefore, kinetic trapping must be caused by solvation effects alone. In order to estimate the magnitude of such an effect, the four lowest-energy conformers presented in Figure 5 have been geometrically optimized with . While in the gas phase conformer I is 1.6 kcal/mol lower in DFT energy than the next minima (conformers IIa and IIb), the situation is reversed when including implicit aqueous solution; conformer I is now 0.9 kcal/mol higher in energy. This suggests that they carry a structural bias from aqueous solution, i.e. the barriers are sufficiently high to kinetically trap them during the electrospray process.
A similar scenario can be seen for conformer III, which is of comparable energy as conformers IIa and IIb on the ∆F (10 K) scale, but the calculated IR spectrum, presented in Figure 5(c), does not match any experimentally observed one. Consulting the ∆F (300 K) scale (see Figure 5(a)) shows that conformer III is 0.9 kcal/mol higher in energy than conformer IIb at room temperature. When re-relaxing the structures to the nearest minimum on the potential energy surface at the PBE0+MBD level including implicit aqueous solvation effects as described above, conformer III becomes further energetically penalized -it is then more than 5.0 kcal/mol higher in energy compared to the other conformers.  There remain two conformers, IV and V, for which the UV spectral signatures have lower intensity (see Figure S2), suggesting that they have smaller populations. The corresponding IR spectra, shown in Figure 6(a), could not be assigned to their calculated counterparts for any structure within 6 kcal/mol from the global minimum on the ∆F (10 K) scale. Similarly, as for IIa and IIb, we assume that these conformers are kinetically trapped, which also renders their assignment difficult in assigning them as these conformers might be higher in energy, and thus no energy criterion can be applied for finding them. Instead we follow an approach 43 where we make use of information from the experiment in order to select from the overall pool of structures for calculation of spectra. Candidates were picked if they feature a free carboxylic acid OH stretch, since the experimental IR spectra show a peak at 3578 cm −1 (see Figure 6(a)). Due to the absence of broad features in the UV spectrum, only structures were considered where the Na + cation was not in close proximity to the phenyl ring. In total, vibrational spectra for 126 conformers have been calculated. In addition to that, local refinement on the PBE0+MBD level for all 52 found minima structures within 3 kcal/mol from the global minimum for the system of AcPheAla 5 LysH + has been laid out after formally replacing Lys with Ala + Na + , with the sodium cation being placed at the position of the amino group nitrogen. Vibrational spectra for the resulting 28 conformers (after clustering) have been calculated as well. As explained above, computationally-costly hybrid xc-functionals are required in order to gain enough accuracy. Only conformer IV (see Figure 6(b)), lying 13.6 kcal/mol higher in energy than the global minimum on the ∆F (10 K) scale, could be assigned to one of the less populated conformers. However, one peak in the simulated vibrational spectrum is blue shifted by 80 cm −1 with respect to the nearest experimental peak, and the corresponding vibrational mode is indicated in Figure 6(b) with a green arrow. Conformer IV is a candidate for the kinetically trapped structure only because of the (partially) matching IR spectra. Taking into account the large computational effort taken, a more appropriate and computationally affordable technique for finding kinetically trapped conformers would be certainly desirable.

Conclusion
Our data indicate that the fixed location of the charge at the C-terminus is imperative for helix formation in peptides of this length in isolation, as this stabilizes the structure through a cation-helix dipole interaction. In the case of the freely-movable sodium cation, the cation-backbone and cation-π interactions seem to be stronger, leading to local distortions of peptide structure, preventing helix stabilization. It is interesting to note the high barriers that seem to be involved in interconverting one structure to another. Even though the cationπ interaction is energetically favored for the AcPheAla 6 + Na + in the gas phase, the system remains kinetically trapped in a structural state that is characterized by cation-backbone interactions and that is energetically preferred in polar solvent.

Acknowledgement
The authors thank the joint Max-Planck-EPFL Center for Molecular Nanoscience and Technology for financial support. CB thanks Dr. Mariana Rossi for sharing her knowledge about theoretical vibrational spectroscopy and Prof. Matthias Scheffler for his continuous support.
The experimental work was supported by the EPFL as well as the Swiss National Science Foundation through grant 200020_165908.

Supporting information available:
• A detailed and technical description of the experiment and of the applied simulations methods as well as additional results can be found in the following sections.
• Numerical details for all shown energy hierarchies as well as experimental and calculated IR spectra, i.e. for Figures 2-6 and Figure 10, as well as corresponding xyz-files of conformers are provided under this link: http://w0.rz-berlin.mpg.de/user/baldauf/carsten_pdf/Schneider_2017_Helices_SI.zip

Experimental setup in detail
The experimental setup has been described in detail elsewhere. 23  Here, they are cooled down by collisions with cold He gas that is pulsed in before their arrival. The He pressure is between 6 · 10 −6 and 10 −5 mbar.

Conformational search approach in detail
The applied conformational search algorithm is similar to the one used by Rossi et al. 16  After clustering, this resulted in 324 (159) conformers for the system of AcPheAla 5 LysH + (AcPheAla 6 + Na + ) in the low-energy region, i.e. within 6 kcal/mol from the global minimum. These conformers were then again locally refined on the PBE0+MBD level using tier 1 basis sets and light settings. After clustering, further geometry relaxation on the PBE0+MBD level using tier 2 basis sets and tight settings resulted in 52 (23) Figure 7: Comparison of energy hierarchies on the PBE0+MBD level (tier 2 basis sets and tight settings) between our own conformational search and the search performed by Rossi et al. 16 Two additional conformers are found in the low-energy region, i.e. within 1 kcal/mol from the global minimum. Conformers with the same energy in both hierarchies correspond to virtually identical structures.
For the system of AcPheAla 5 LysH + , nine conformers were found in the low-energy region, i.e. within 1 kcal/mol from the global minimum, two more than Rossi et al. 16 Figure 7 compares the two hierarchies on the potential energy surface, i.e. on the PBE0+MBD level using tier 2 basis sets and tight settings. Conformers with the same energy in both hierarchies correspond to virtually identical structures.

Helmholtz free energy
The Helmholtz free energy F per molecule in the gas phase is given by with E denoting the DFT energy on the potential energy surface (PES) and F int denoting the free energy contribution due to the internal degrees of freedom, consisting of vibrations and rotations. Assuming harmonic approximation for the intramolecular PES and neglecting any rotational-vibrational coupling, the internal free energy is given by and where N denotes the number of atoms, ω i denotes the vibrational frequency of normal mode i, and I x , I y , and I z denote the moments of inertia along the three axes. The case of T = 0 defines the zero-point energy correction where F int = 3N −6 ih ω i 2 . The Helmholtz free energy F is formally given by F = U −T S, with U , T , and S denoting the internal energy, the temperature, and the entropy of the system, respectively. It is related to the Gibbs free energy G through G = U − T S + pV = F + pV , with p and V denoting the pressure and the volume of the system, respectively. We use Helmholtz free energies in this work because the experiment is essentially done at zero pressure. Furthermore, throughout this work we are exclusively treating relative energies, i.e. comparing energy differences between different conformers (usually with respect to the global minimum) of the same system. Hence, the the pV term cancels. In other words, ∆G = ∆F .  Figure 8: Measured UV spectra for the systems of (a) AcPheAla 5 LysH + and (b) AcPheAla 6 + Na + . In both cases, peaks have been assigned to their identified conformers shown in Figures 3, 5, and 6.
Measured UV spectra for the systems of (a) AcPheAla 5 LysH + and (b) AcPheAla 6 + Na + are presented in Figure 8. In both cases, peaks have been assigned to their identified conformers shown in Figures. 3, 5, and 6 in the paper.
AcPheAla 6 + Na + : Attempting assignemt of IR spectra of kinetically trapped conformers IV and V As stated in the paper, for the system of AcPheAla 6 + Na + obtained IR spectra of kinetically trapped conformers IV and V could not be confidently assigned to their calculated counterparts. For the sake of completeness, Figure 9 shows the corresponding experimentally obtained IR spectra along with several calculated IR spectra that share similarities.
Calculations have been done on the PBE0+MBD level using tier 1 basis sets and light settings. Conformer charmm22Best500FF155865 corresponds to conformer IV illustrated in Figure 6(b) in the paper.  In both cases, several calculated IR spectra that share similarities are shown. Conformer charmm22Best500FF155865 corresponds to conformer IV illustrated in Fioure 6(b) in the paper. xyz-files of the conformers are provided within the SI_data.zip of this SI.

DFT calculations including implicit solvation effects
For the system of AcPheAla 6 + Na + and the four lowest-energy conformers on the ∆F (10 K) scale, re-relaxation was applied on the potential energy surface on the PBE0+MBD level