Letter pubs.acs.org/JPCL
Does Thermal Breathing Affect Collision Cross Sections of Gas-Phase Peptide Ions? An Ab Initio Molecular Dynamics Study Robert Pepin, Alessio Petrone, Kenneth J. Laszlo, Matthew F. Bush, Xiaosong Li, and František Tureček* Department of Chemistry, University of Washington, Bagley Hall, Box 351700, Seattle, Washington 98195-1700, United States S Supporting Information *
ABSTRACT: Ab initio molecular dynamics (AIMD) with density functional theory (DFT) was applied to explore conformational motions and collision cross sections (Ω) of folded (2) and extended (7) conformers of doubly charged peptide ions, (Ala−Ala−Leu−Arg + 2H)2+, in the gas phase at 300 and 473 K. The experimental Ω of (Ala−Ala−Leu−Arg +2H)2+ was measured as 149 ± 1.2 Å2 at 298 K. Thermally distributed mean values of Ω for 2 and 7 at 300 and 473 K were only 0.8−1.1% larger than for the equilibrium 0 K structures. Long (>10 ps) trajectory calculations indicated entropy-driven conformational change of 2 to 7 that occurred at random within a ∼ 4 ps time window. The experimental Ω was found to fit the calculated population averaged values for 2 and 7, indicating a rapid conformer interconversion. Overall, thermal breathing had only a minor effect on the peptide ion collision cross sections.
M
trajectory (TM) model11 that has been refined for different collision gases5 and here is called TM*. The conversion to Ωth of a molecular structure critically depends on the veracity of the calculated structure and the choice of the collision model. Small ions with up to ca. < 200 atoms can be treated by ab initio or density functional theory (DFT) calculations to provide geometries of local energy minima and atomic charges, thus allowing the use of TM methods to calculate Ωth. With larger ions, structures are usually obtained by force-field molecular dynamics calculations and converted to Ωth by a method not using ion−molecule potentials.12−14 This approach has the disadvantage of not using equilibrium geometries for gas-phase ions and not considering intermolecular interaction potentials between the ion and the gas molecules. The DFT-TM approach has been shown to produce Ωth within a few percent of the experimental Ω for several peptide ions.15,16 However, an inherent weakness of DFT-TM is that it uses 0 K equilibrium structures to approximate the properties of gas-phase ions in experiments that are performed at ambient or elevated temperatures12−14 where a number of peptide vibrational modes are active. Vibrational effects on ion drift mobilities have been studied for rigid ion structures such as C60 and found to be weak.11,17,18 However, peptide ions exhibit a range of conformational flexibilities due to low-frequency vibrational modes that may affect the ion collision cross section. Furthermore, applied studies19 often use dynamic methods of ion mobility mass spectrometry, such as the traveling wave technique (TWAVE).5,20 Ion storage and injection into TWAVE have
olecular shape is an important characteristic that affects physical properties, reactivity, and molecular recognition of neutral molecules and ions. Regarding isolated molecules or ions in the gas phase, the shape of an ion can be characterized by its collision cross section (Ω) that is related to the ion velocity when drifting through a gas in a weak electrostatic field.1 With the advent and developments of ion mobility mass spectrometry,2,3 there has been increased interest in the experimental and theoretical determination of collision cross sections for gas-phase ions ranging from ionized drug molecules to protein complexes. Experimental methods for determining Ω are based on drift time measurements of ions accelerated by electrostatic fields, according to the Mason−Schamp equation for ion motion (1),1 where e is the elementary charge, z is the ion charge state, N is the drift gas number density, μ is the reduced mass of the ion-collision gas atom pair, kB is the Boltzmann constant, T is the gas absolute temperature, and K is the ion mobility 1/2 3ez ⎛ 2π ⎞ 1 Ω= ⎜ ⎟ 16N ⎝ μkBT ⎠ K
(1)
Because Ω is a scalar quantity, it is desirable to relate it to the ion 3D structure that can be obtained by electronic structure calculations. To this end, the coordinates and charges of atoms in an ion can be used to calculate its theoretical collision cross section (Ωth) in a given gas. An efficient software (MOBCAL)4 has been introduced and developed5,6 to allow calculations of Ωth at different levels of approximation. The standard models for Ωth calculation, dealing with atom coordinates only, are the elastic hard sphere scattering (EHSS)7 and projection average approximation (PA)8−10 methods. Intermolecular interaction potentials are considered in the more sophisticated ion © XXXX American Chemical Society
Received: May 31, 2016 Accepted: July 7, 2016
2765
DOI: 10.1021/acs.jpclett.6b01187 J. Phys. Chem. Lett. 2016, 7, 2765−2771
Letter
The Journal of Physical Chemistry Letters
change upon adjusting the width of the mobility injection time from 200 to 50 μs. The experimental data indicate that (AALR + 2H)2+ was traveling through the ion mobility drift cell as a single conformer or a mixture of conformers of very similar collision cross sections. The experimental cross section of (AALR + 2H)2+ with He was determined as Ω = 149 ± 1.2 Å2 from three measurements. The fwhm of the ion mobility profile on the Ω scale was 10.8 Å2. Full geometry optimization found nine 0 K equilibrium structures of lowest energy (1−9) within an energy range of ca. 30 kJ mol−1, according to B3LYP calculations (Table 1). The optimized ion geometries substantially differed in their folding patterns, which ranged from tightly folded conformers 1−5 to largely unfolded ones (8, 9) (Figure 2). The folded conformers shared some overall structure features in that the N-terminal ammonium of Ala-1 was hydrogen bonded to the carboxyl C O of Arg-4 and the protonated Arg side chain group was hydrogen bonded to the Leu-3 amide carbonyl. Differences were found in the carboxyl conformation (s-cis in 3 and 4 and strans in 1, 2, and 5) and H-bonding patterns of the Arg NH (1−3) and NH2 (4, 5) groups. The least folded conformers 8 and 9 had the N-terminal ammonium solvated by the proximate Ala-1 amide carbonyl, whereas the inner N−H (9) and outer NH2 (8) guanidinium protons were H-bonded to the Arg-4 carboxyl. Conformers 6 and 7 had the same arrangement of the N-terminal ammonium group as in 8 and 9, but the Arg N−H (7) and NH2 (6) protons made hydrogen bonds to the Leu-3 amide carbonyl. The 0 K relative enthalpies of conformers 1−9 showed a general trend of stabilization of the more folded structures 1, 2, and 5 when based on the M06-2X, ωB97X-D and MP2 calculations. In contrast, B3LYP favored the unfolded structures 6−9 (Table 1). The relative energies resulted from balancing two major terms, one of which was the potential-energyincreasing Coulomb repulsion of the charged groups that increases as the groups are spaced closer to each other because of ion folding. The other factor was the attractive ion-dipole and dipole−dipole interactions of hydrogen bonds that were more abundant in the folded structures. The 0 K relative enthalpies (Table 1) were further substantially modified by including 298 K ion entropies that favored structures 3 and 6− 9, which had a larger number of low-frequency normal modes (ν < 100 cm−1) than the more tightly folded structures 1, 2, 4,
been reported to increase the effective temperature of peptide ions and thus potentially affect their stability, conformations, and cross sections.21,22 Here, we investigate the effect of ion internal energy on the collision cross sections of gas phase peptide ions. To this end, we chose doubly protonated Ala−Ala−Leu−Arg (AALR + 2H)2+ ions that have been previously shown to assume substantially different but nearly isoenergetic conformations.23 We compare the experimental Ω to theoretical values for ion conformers that were calculated for geometries corresponding to local energy minima at 0 K and for structures generated by ab initio molecular dynamics calculations at 300 and 473 K. Details of the computational procedures are given in the Supporting Information. We wish to show that including the thermal motions in gas-phase ions has a only minor effect on the theoretical collision cross sections. Ion mobility measurements of (AALR+2H)2+ in He gas gave the arrival time profile shown in Figure 1. This shows a single
Figure 1. Arrival time profile of (AALR + 2H)2+ ions in the ion mobility drift cell.
Gaussian peak with a full width at half-maximum (fwhm) of 0.16 ms which was obtained from a least-squares Gaussian fit to the experimental data, giving root-mean-square deviation, rmsd = 1.8%. Comparison of the arrival time profile for (AALR + 2H)2+ with that of a calibration peptide ion (GRGDS + H)+ showed no broadening that would indicate the presence of multiple overlapping peaks. The arrival time profile did not Table 1. Relative Energies of (AALR + 2H)2+ Conformers
relative energya,b B3LYP
c
ωB97X-De
M06-2Xd
MP2c
ion
ΔH°0
ΔG°298f
ΔH°0
ΔG°298f
ΔH°0
ΔG°298f
ΔH°0
ΔG°298f
1 2 3 4 5 6 7 8 9
10.2 0.0 −1.9 3.7 3.3 −6.8 −14.7 −2.2 −6.6
11.6 0.0 −16.4 0.5 4.7 −37.9 −41.1 −26.6 −27.7
1.9 0.0 18.9 17.2 1.3 33.8 20.8 35.5 32.7
3.3 0.0 4.5 14.0 2.6 2.7 −5.6 11.1 11.6
−0.5 0.0 10.9 11.2 3.0 25.1 12.7 33.0 26.5
0.8 0.0 −3.6 8.0 4.3 −6.0 −13.7 8.5 5.5
2.8 0.0 12.0 12.4 3.9 27.3 14.7 39.9 33.5
4.2 0.0 −2.5 9.2 5.3 −3.8 −11.7 15.4 12.5
a In units of kJ mol−1. bIncluding zero-point vibrational energy corrections and referring to 0 K. cFrom single-point energy calculations with the 6311++G(2d,p) basis set on B3LYP/6-31+G(d,p) optimized geometries. dFrom single-point energy calculations with the 6-311++G(2d,p) basis set on M06-2X/6-31+G(d,p) optimized geometries. eFrom single-point energy calculations with the 6-311++G(2d,p) basis set on ωB97X-D/631+G(d,p) optimized geometries. fIncluding 298 K enthalpies and entropies with corrections for hindered rotors.
2766
DOI: 10.1021/acs.jpclett.6b01187 J. Phys. Chem. Lett. 2016, 7, 2765−2771
Letter
The Journal of Physical Chemistry Letters
Figure 2. M06-2X/6-31+G(d,p) optimized structures of ions 1−9. The atom color-coding is as follows: cyan = C, red = O, blue = N, gray = H. Major hydrogen bonds are annotated by double-headed green arrows.
and 5. The vibrational entropy effects resulted in a substantial lowering of the ΔG298 of conformers 3 and 6−9. The order of the ΔG298(ion) also depended on the computational method; the trend ΔG298(7) < ΔG298(6) < ΔG298(3) < ΔG298(2) was supported by the majority of methods we used (Table 1). However, the overall differences in the ΔG298(ion) were rather small and because of the large entropy differences among the isomers it was problematic to firmly assign the lowest-energy conformer. The geometries and atomic charges in 1−9 were used to calculate theoretical collision cross sections (Ωth) using the projection average approximation (PA) and ion trajectory (TM and TM*) methods. The TM* calculations employed interaction potentials optimized for small ions by Campuzano et al.5 The calculated Ωth (Table 2) followed an obvious trend of Ωth(PA) < Ωth(TM*) < Ωth(TM). The root-mean square deviations (rmsd) showed the Ωth(PA) being on average 5.3 Å2 (3.6%) smaller than the Ωth(TM*) whereas the Ωth(TM) were on average 4.3 Å2 (2.8%) larger than the Ωth(TM*). In addition, the Ωth that were calculated for the M06-2X optimized ion structures were on average 2.7 Å2 (1.7%) smaller than those for B3LYP structures, indicating that a tighter hydrogen bonding was favored by M06-2X. This is consistent with the ΔH0 energy data (Table 1) where M06-2X favored the more folded ion structures. Similar differences between Ωth from B3LYP and M06-2X calculations have been reported for other peptide ions15,16 indicating that this trend may be of a general nature. Considering the spread of the calculated Ωth,
Table 2. Theoretical Collision Cross Sections of (AALR + 2H)2+ Ion Conformers with He Gas collision cross section (Å2) B3LYP
M06-2X
ion
PA
TM
TM*
PA
TM
TM*
1 2 3 4 5 6 7 8 9
138.7 142.5 143.6 141.8 142.5 154.7 153.3 159.2 158.1
147.0 151.3 154.3 150.7 156.6 164.6 163.4 166.7 165.6
144.5 147.3 150.7 148.1 147.4 158.7 158.3 161.6 161.0
135.9 140.2 140.6 138.8 140.2 152.2 150.5 157.9 156.6
143.8 148.8 150.4 146.9 148.6 161.5 160.7 164.7 163.6
141.5 145.3 147.7 145.1 144.8 155.7 156.2 159.7 159.1
the low-energy conformers 2−7 can be considered as plausible candidates to fit the experimental Ω (149 Å2). For AIMD calculations, we selected conformers 2 and 7 representing, respectively, the folded and unfolded low-energy peptide ions. Because AIMD was run with M06-2X/6-31G(d), we first checked if the smaller basis set had any effect on the starting geometries. Also, the effect of dispersion interactions was examined by including Grimme’s D3 dispersion corrections24 in the geometry optimization procedure. Table S1 and S2 (Supporting Information) show the comparison of dihedral angles and interatomic distances obtained for 2 and 7 by the different gradient optimization procedures. The data showed that the optimized structures were closely similar, 2767
DOI: 10.1021/acs.jpclett.6b01187 J. Phys. Chem. Lett. 2016, 7, 2765−2771
Letter
The Journal of Physical Chemistry Letters indicating that the less computationally demanding M06-2X/631G(d) procedure can adequately capture the ion geometry. In particular, the parameters expressing the distance between the N-terminal N1 and the carboxyl O43 atoms were within 0.02 Å from the optimizations. This parameter, R(N1−O43), along with the maximum nonbonding distance between the remotest atoms (Rmax) are quite distinct for 2 and 7 (Table S1, S2) and therefore were used to gauge the thermally induced conformational motion in the ions in the course of the AIMD run. The normalized distributions of Rmax obtained from 99 to 106 AIMD trajectories for 2 and 7 are shown in Figure 3. After
Figure 3. Normalized distribution of the maximum nonbonding distances evaluated at each step of AIMD trajectories at 300 and 473 K. Black line, ion 2 at 300 K, 106 trajectories; red line, 2 at 473 K, 99 trajectories; blue line, ion 7 at 300 K, 101 trajectories; pink line, 7 at 473 K, 106 trajectories. Black dashed line, 13 ps trajectory for 2 at 300 K; red dashed line, 20 ps trajectory for 2 at 473 K. Arrows show the static Rmax in 2 (12.56 Å) and 7 (14.97 Å).
Figure 4. Thermal distributions and population mean values of TM* calculated collision cross sections in 300 K He drift gas of ions 2 and 7 at (a) 300 K and (b) 473 K ion internal temperature. (c) Collision cross sections in 473 K He drift gas of ions 2 and 7 at 473 K ion internal temperature.
1.1 ps at 300 K, the folded ion 2 shows a distribution of Rmax extending from 12 to14.3 Å with a maximum at 13.35 Å. This corresponds to a 6.3% increase of Rmax relative to the equilibrium 0 K structure. The trajectories run at 473 K showed a distribution of Rmax extending to 14.6 Å with a maximum at 13.5 Å. Extending the running time for the AIMD trajectories resulted in a further minor shift and broadening of the Rmax distribution. The AIMD trajectories for the less folded ion 7 also showed an upward shift of the Rmax distribution. In this case the 300 and 473 K distributions had nearly identical maxima at 15.45 Å. This corresponds to a 3.2% increase relative to the 0 K equilibrium structure. The temperature had an effect on the width of the Rmax distribution that was broader at 473 K, reaching from 13.3 Å through 17.2 Å, whereas the 300 K distribution was between 14 and 17 Å. The Rmax values for conformers 2 and 7 overlapped, and the overlap interval increased with temperature. Overall, the shift of the Rmax with temperature was greater for the folded conformer 2 whereas the distribution width was greater for 7. To analyze the effect of thermal motion on the ion collision cross sections, AIMD snapshots along the trajectories were sampled at regular intervals, and 900−970 extracted structures were used to calculate Ω tr with MOBCAL. The Ω tr distributions calculated by the TM* model are shown in Figure 4, and extracted values for the population maxima and distribution full widths are collated in Table 3. The calculated
data for 2 and 7 indicate very small thermal effects on the calculated Ω. For 2, the Ωtr at 300 K is only 1.6 Å2 (1.1%) larger than the static Ω and does not change at 473 K. Likewise for 7, the 300 K Ωtr is 1.1 Å2 (0.7%) larger than the static value and does not change at 473 K. Temperature effects are seen for the Ωtr distribution widths (fwhm) that for 2 increase from 3.5 to 3.9 Å2 between 300 and 473 K. This change is only marginally greater for 7 where the Ωtr fwhm increases from 4.4 to 5.1 Å2 between 300 and 473 K. Somewhat larger effects are found for the full distribution ranges of Ωtr. With 2 the 300 K range (141−150 Å2) increases to 138−155 Å2 at 473 K and likewise for 7 where the Ωtr range goes from 150 to 162 Å2 to 151−168 Å2 at 300 and 473 K, respectively. The AIMD results for doubly charged ions 2 and 7 are consistent with previous studies of temperature effects on ion mobility of other gas-phase ions. In a study of sodiated polyethylene glycol ions, Bowers and co-workers found very similar Ω for structures that were based on force-field molecular modeling at 0 and 300 K.18 Jarrold and co-workers reported a series of peptide ions to which they assigned globular and αhelix conformations on the basis of collision cross sections.25 They found that the collision cross section of α-helical [Ac(AGG)5K + H]+ ions decreased with the drift gas temperature whereas the globular ions were remarkably insensitive to the experimental temperature in the range of 295−409 K. Russell and co-workers have reported effects of 2768
DOI: 10.1021/acs.jpclett.6b01187 J. Phys. Chem. Lett. 2016, 7, 2765−2771
Letter
The Journal of Physical Chemistry Letters Table 3. Collision Cross Section Means and Ranges for Ions 2 and 7.a Collision Cross Section (Ωth, Å2) Mean Ion 2 7
300 K
a
473 K d
146.7 (145.1) 156.8 (155.6)d
Rangea b
146.8 156.7
473 K
c
137.8 147.5
300 K
473 Kb
473 Kc
141−150 150−162
138−155 151−168
130−145 139−153
a
From TM* calculations on AIMD-M06-2X/6-31G(d) geometries and atomic charges. bIon 473 K rovibrational temperature in 298 K drift gas. cIon and drift gas temperature at 473 K. dΩth of static M06-2X/6-31G(d) optimized structures.
Figure 5. M06-2X/6-31G(d) structure snapshots along the AIMD trajectory starting from 2. Yellow double-ended arrows indicate transient hydrogen bonds developing in the course of ion unfolding. Green double-ended arrows indicate persistent hydrogen bonds. Distances are given in Ångstrøms.
K (Figure 4c). In addition, the difference between the population averaged Ωth for 2 and 7 decreased from 10 Å2 in He at 300 K to 8.7 Å2 in He at 473 Å2. Hence, in addition to the adverse decrease in resolution in ion mobility measurements with increasing temperature due to diffusion,27 the present computational analysis predicts a temperature-dependent compression of the peptide ion cross sections. The temperature effect on ion mobility of C60 ions and the role of intermolecular potentials were analyzed in detail by Jarrold et al.11 In contrast to collision-gas temperature effects, conformational changes caused by temperature and ion motions through gases and resulting in major changes in Ω have been amply documented for gas phase peptide and protein ions.14,29−32 We examined the long trajectories starting from the folded conformer 2, showing the evolution of R(N1−O43) that measures the distance between the carboxyl and N-terminal amine groups in the ion (Figure S1a, Supporting Information). The 300 K trajectory shows an oscillation of R(N1−O43) for up to 13 ps after which the structure suddenly unfolds within ca. 4 ps. This development is equivalently represented by the Ωth calculated with TM* for snapshots sampled at 0.2 ps
lowering the drift gas temperature to 80 K and demonstrated improved ion mobility resolution, although effects on the Ω have not been studied.26,27 Our TM* calculations refer to collisions with thermal (300 K) He atoms whereas the drifting ion is vibrationally excited to the defined effective temperature which is 300 or 473 K. Relating to the experimental conditions (2.5 Torr He pressure at 298 K), ions 2 and 7 are calculated to undergo ca. 337 000 collisions with He atoms during their 2.2 ms drift time. This gives the average time between the collisions as Δtcoll = 0.0022/ 337 000 = 6.5 × 10−9 s, which is substantially longer than the 1.1 ps interval during which the distribution of Ωth develops by ion conformational breathing. This indicates that the ion internal motion is substantially faster than the collision rate and the calculated Ωth must be population averaged28 (Table 3). This averaging involves both the vibrationally stretched and compressed conformers, resulting in only minor changes in the average Ωth relative to those for the 0 K local energy minima. In contrast to the minor ion temperature effects on the Ωth, TM* calculations of collisions with He at the equilibrium gas temperature of 473 K show ca. 6% smaller cross sections for [AALR + 2H]2+ conformers 2 and 7 compared to those at 300 2769
DOI: 10.1021/acs.jpclett.6b01187 J. Phys. Chem. Lett. 2016, 7, 2765−2771
Letter
The Journal of Physical Chemistry Letters intervals (Figure S1c). The Ωth of 2 oscillates about the initial mean value of 145.8 ± 1.5 Å2 for ca. 13 ps and then rapidly increases to 159.9 ± 1.9 Å2 due to peptide ion unfolding. The unfolding steps are visualized with snapshot structures in Figure 5. Starting from 2, the general features of its conformation are retained for at least 12 ps, showing hydrogen bonds between the C-terminal COOH and the N-terminal NH3+ and Ala-2 amide groups, as well as between the Arg guanidinium and the Leu-3 amide in the structure captured at 1 ps. In the critical region starting at 12.7 ps, the Arg H-bonding is disrupted by thermal motion of the Arg side chain, and the other two Hbonds are stretched to 2.51−2.75 Å. The NH3+ carboxyl Hbond is irreversibly disrupted at 12.9 ps to unfold the Nterminal Ala residue. In the further stages captured at 13.1− 14.1 ps, the COOHAla-2 amide H-bond oscillates while the Arg side chain undergoes multiple rotations to finally create a H-bond to the carboxyl group. The structures developing after 16.9 ps are similar to that of the extended conformer 9. Note that 9 is not the lowest ΔG298 conformer according to M06-2X calculations, and so it is not excluded that it will undergo further refolding at longer times. The conspicuous feature of the AIMD trajectory calculations is the fast conversion from conformers resembling 2 to those resembling 9 that occurs within ca. 4 ps. The main stages of the unfolding process occur at even shorter times between 12.7 and 14.1 ps. The unfolding involves several stretching and bending normal vibrational modes in 2. Within 1.4 ps, modes with v > 24 cm−1, which represent 194 out of all 195 modes in 2, undergo one or more vibrational periods. This is consistent with the AIMD analysis that multiple mode coupling leads to conformational changes within 1.4 ps. It should be noted that the above analysis is based on a single AIMD trajectory, and so different initial conditions of vibrational energy distribution were not examined. The main limit to more extensive trajectory analysis is the number of DFT energy and gradient evaluations (100 000 for a 20 ps run), which become prohibitively expensive for running multiple trajectories, as done for the 1.1 ps AIMD runs. Thus, the results from the single run cannot be used to establish the unfolding kinetics. This is illustrated by a single trajectory AIMD run of 2 at 473 K (Figure S1b), which shows fluctuations of the R(N1− O43) distance with occasional peaks stretching to 4.8 Å, but no overall structure unfolding analogous to that indicated by the 300 K trajectory. It is presumed that the probability for unfolding and refolding events will increase during the much longer time period pertinent to the experimental conditions where the ions are trapped for
