Anharmonic Vibrational Analyses of Pentapeptide Conformations

Subscriber access provided by EPFL | Scientific Information and Libraries

Article

Anharmonic Vibrational Analyses of Pentapeptide Conformations Explored with Enhanced Sampling Simulations Hiroki Otaki, Kiyoshi Yagi, Shun-ichi Ishiuchi, Masaaki Fujii, and Yuji Sugita J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b06672 • Publication Date (Web): 13 Sep 2016 Downloaded from http://pubs.acs.org on September 16, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Anharmonic Vibrational Analyses of Pentapeptide Conformations Explored with Enhanced Sampling Simulations #

Hiroki Otaki,†, Kiyoshi Yagi,†,¶ Shun-ichi Ishiuchi,‡ Masaaki Fujii,‡and Yuji Sugita*†,¶,§,∥

†RIKEN Theoretical Molecular Science Laboratory, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan. ¶RIKEN iTHES, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan.

‡Laboratory for Chemistry and Life Science, Institute for Innovative Research, Tokyo Institute of Technology, 4259 Nagatsuta-cho, Midori-ku, Yokohama 226-8503, Japan. §RIKEN Advanced Institute for Computational Science, 7-1-26 Minatojima-Minamimachi, Chuo-ku, Kobe, Hyogo 650-0047, Japan. ∥RIKEN Quantitative Biology Center, 1-6-5 Minatojima-Minamimachi, Chuo-ku, Kobe, Hyogo 650-0047, Japan.

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 47

ABSTRACT

An accurate theoretical prediction of the vibrational spectrum of polypeptides remains to be a challenge due to (1) their conformational flexibility and (2) non-negligible anharmonic effects. The former makes difficult the search of conformers that contribute to the spectrum, and the latter requires an expensive, quantum mechanical calculation for both electrons and vibrations. Here, we propose a new theoretical approach, which implements an enhanced conformational sampling by the replica-exchange molecular dynamics (REMD) method, a structural clustering to identify distinct conformations, and a vibrational structure calculation by the second-order vibrational quasi-degenerate perturbation theory (VQDPT2). A systematic mode-selection scheme is developed to reduce the cost of VQDPT2 and the generation of a potential energy surface by the electronic structure calculation. The proposed method is applied to a pentapeptide, SIVSFNH2, for which the infrared (IR) spectrum has been recently measured in the gas-phase with high resolution in the OH and NH stretching region. The theoretical spectrum of the lowest-energy conformer is obtained with a mean absolute deviation of 11.2 cm-1 from the experimental spectrum. Furthermore, the NH stretching frequencies of the five lowest-energy conformers are found to be consistent with the literature values measured for small peptides with a similar secondary structure. Therefore, the proposed method is a promising way to analyze the vibrational spectrum of polypeptides.


2

Page 3 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION Vibrational spectroscopy is a viable tool to elucidate the structure and dynamics of biomolecules. Polypeptides and proteins have been studied by the infrared (IR) spectroscopy1-4 utilizing the amide I band (the C=O stretching mode) that brings about a strong IR signal. The amide A band (the NH stretching mode) is also a sensitive reporter of hydrogen bonds (HBs) and secondary structures, but was hard to measure due to an overlap with a broad OH stretching band of bulk water. Recently, Yan and coworkers5-7 have shown that in situ, real-time detection of the amide A band is feasible for proteins at a lipid-water interface by a chiral sum-frequency generation (SFG) spectroscopy, which makes use of a chiral structure of the peptide backbone and a surface selectivity of the SFG method. The measurement in the gas phase, i.e., in an isolated, solvent-free environment, also offers a highly resolved, conformer-specific spectrum in the amide A region using the double-resonance techniques.8-10 Notably, the laser desorption and the electrospray ionization11-12 have opened the way to produce an abundant amount of peptides in a supersonic jet, which can be subsequently analyzed by the laser spectroscopy. The early works have characterized the HB of amide backbones in short peptides (2 or 3 residues).9, 13-16 The challenge extends to more complex, real systems, e.g., the larger peptides,17-22 the interaction with water or micro-hydration,23-27 and a complex with hydrocarbon molecules.28 Despite the significant advances in the spectroscopic techniques, one of the remained issues is the assignment of the observed vibrational bands. For this purpose, theoretical assistance is highly desirable; however the computation of a vibrational spectrum is a challenge because: (1) polypeptides have large conformational flexibility, which makes it difficult to find the structure that contributes to the spectrum; (2) a quantitative prediction of band positions and intensities re-


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 47

quires a high level of vibrational and electronic structure theories. In particular, high frequency OH/NH stretching modes need a sophisticated treatment due to non-negligible anharmonic effects. In the previous works,29-30 conformational space of peptides has been sampled using molecular dynamics (MD) simulations. After MD simulations to find candidate structures, accurate electronic structure calculations follow to refine the rank of relative energies for peptide structures. Although the scheme has been applied to small peptides so far, the conventional MD can be easily trapped at a local minimum and hardly find the global minimum structure, causing a bias in the initial set of candidate structures. Replica-exchange molecular dynamics31-33 (REMD) method, which is one of the most widely used generalized-ensemble methods, is capable of sampling the wide conformational space of peptides by running multiple copies (or replicas) of MD simulations at different temperatures and by exchanging a pair of replicas so as to achieve a random walk in the temperature space. Recently, Schubert et al.34 have carried out REMD simulation to explore the conformation and ab initio MD simulation to calculate the IR spectrum of 20-residue peptides, Ac-Ala19-Lys+H+ and Ac-Lys-Ala19+H+. They have revealed that, while the former rigidly keeps the helix structure, the latter is more flexible having several conformations in the lowenergy regime. Once a few lowest-energy conformations are obtained from the classical MD simulations, their vibrational spectra are calculated to compare with the experiment and to assign the vibrational band. Most of the previous studies have calculated the spectrum based on the harmonic approximation. The common practice is to scale the harmonic frequency by a factor of 0.95 – 0.99, as the approximation usually overestimates the frequency. Recently, Bouteiller et al.35-36 proposed mode-specific scaling factors to improve the agreement with experiment. However, the neglect


4

Page 5 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of anharmonicity may incur not only a quantitative error, but also a qualitative failure yielding a wrong number of bands. In the presence of resonance, the anharmonic coupling splits the fundamental transition into several bands, which is often observed for hydrogen-bonded NH/OH groups.16, 37 Although the explicit account of anharmonicity used to be costly, recent advances in the vibrational self-consistent field38-39 (VSCF) and post-VSCF40-46 methods have made an application to large molecules feasible. Roy et al.47 have recently calculated the vibrational spectrum of gramicidin S (10 residues) by the second-order perturbative extension to VSCF (VSCFPT2).41 Panek and Jacob,48 and Cheng and Steele49 have proposed an efficient method based on localized vibrational modes and applied to polypeptides. Two of the authors (HO and KY) have recently developed the second-order vibrational quasi-degenerate perturbation theory50-51 (VQDPT2), which is a perturbative extension to VSCF with an account of resonance (or quasidegenerate states). Together with an anharmonic potential generated by the electronic structure calculation, the method is capable of predicting the vibrational spectrum of hydrogen bonded systems with high accuracy.37 The present work is the first report, to the best of our knowledge, that implements both the enhanced conformational search and the anharmonic treatment of vibrational states of polypeptides. In this study, we perform the REMD simulation and the VQDPT2 calculation to predict the vibrational spectrum of a penta-peptide, SIVSF-NH2, shown in Figure 1. The IR spectrum of SIVSF-NH2 has been recently measured by two of us (SI and MF) in a range of 3200 – 3700 cm1

for the OH and NH stretching modes.21 Here, we show that the k-means clustering algorithm52

for the REMD simulation trajectory can efficiently select candidate structures that yield the experimental spectrum. SIVSF-NH2, which has 81 atoms and 237 vibrational degrees of freedom, is a large system for the VQDPT2 calculation. Therefore, we have developed a systematic way to


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 47

construct a reduced dimensional model, which incorporates only the vibrational modes that are strongly coupled to the modes of interest (i.e., OH/NH stretching modes). We show that the method reproduces the experimental spectrum with high accuracy and provides an unambiguous assignment of the vibrational bands.

O H2 N

CH C

O H N

CH C

O H N

CH C

O H N

CH C

CH2

CH CH3

CH CH3

CH2

OH

CH2

CH3

OH

O H N

CH C

NH2

CH2

CH3

Ser

(1)

Ile

Val

Ser

(4)

Phe

Figure 1. The chemical structure of SIVSF-NH2.


6

Page 7 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


METHODS The features of the present work are schematically shown in Figure 2. One is an efficient conformational search based on REMD and the other is a mode-selection scheme that makes VQDPT2 calculations for large molecules feasible.

(a)

(b) Conformational Search Sampling (REMD) 30,000 structures

Selection of Vibrational Modes frozen (157 modes) small coupling

Clustering Analysis (k-means clustering) 144 structures Geometry Optimization (DFT) 5 structures

active (70 modes) large coupling target modes (OH/NH stretching modes)

Vibrational Analysis (VQDPT2)

Figure 2. Schematic illustration of (a) the conformational search and (b) the selection of vibrational modes.

Conformational Search. Conformations of SIVSF-NH2 are searched in a scheme shown in Figure 2 (a). First, we carry out a REMD simulation to sample a large number of peptide con-


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 47

formations in gas. Then, the simulation snapshots are classified into some clusters by the kmeans clustering analysis and a representative structure is taken for each cluster. After density functional theory (DFT) calculations to refine the structure and energy, VQDPT2 calculations are carried out for a few lowest-energy structures to obtain the IR spectrum. The details of each step are described in the following. The system for REMD simulation was a SIVSF-NH2 molecule in vacuum in accordance with the experimental condition.21 A fully extended structure was prepared as an initial structure using the MMTSB Toolset.53 The force-field parameters for the peptide were taken from CHARMM36 with CMAP correction.54 Following the energy minimization using CHARMM,55 the REMD simulation with 12 replicas was carried out in a temperature range distributed between 300 and 1300 K using NAMD (version 2.9).56 The peptide was first equilibrated for 1 ns at each temperature. Then, a REMD simulation was carried out for 60 ns per replica. When the replica exchanges were applied every 2 ps, the acceptance ratio was 28 % on average [see Figure S1 in the supporting information (SI)]. The time step for MD simulation in each replica was set to 1 fs. The non-bonded interaction was calculated every step with no truncation. Langevin thermostat was used to control the target temperatures with the coupling time set to 5 ps-1. MD trajectory snapshots were saved every 2 ps. The 30,000 snapshot structures were used for the analysis. The k-means clustering analysis was performed using the MMTSB Toolset for the REMD trajectory analysis. The root-mean-square deviation (RMSD) between the snapshot structures in terms of Cartesian coordinates of all atoms was used as a measure for the fixed radius clustering. The radius was set to 1.7 Å (see Section 3 for justification) which yielded 144 clusters. For each


8

Page 9 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cluster, the snapshot structure closest to the centroid or the average structure was selected as a representative structure. The 144 structures were used as an initial structure for geometry optimizations at the level of B3LYP exchange-correlation functionals57-58 with a mix usage of 6-31G** and 6-31++G** basis sets59-61 (denoted 6-31(++)G**). The diffuse basis functions were added only to the nitrogen and oxygen atoms and the hydrogen atoms bound to them. The electronic structure calculations were performed using Gaussian09.62 Finally, the anharmonic vibrational analysis was carried out for the five lowest-energy structures. VQDPT2 Calculation with Selected Number of Modes. SIVSF-NH2, which has 81 atoms and 237 vibrational modes, is a large molecule for vibrational structure calculations. In order to reduce the cost, we take advantage that the target modes in the present study are only the ten highest-frequency modes, i.e., the OH and NH stretching modes. As shown in Figure 2 (b), we divide the system into target modes, active modes that strongly couple with the target modes, and others that are kept frozen. The procedure to select the active modes is as follow: 1) Calculate the derivatives of the potential up to the fourth order in terms of normal coordinates (requires 2f +1 points of Hessian calculations, where f is the number of modes)63 and a measure of two- and three-mode coupling strength (MCS) proposed by Seidler et al.64 2) List all combination of two and three modes, in which at least one of the modes is the target mode, in a decreasing order of MCS. 3) Select the modes in the combination from the top of the list, until a maximum number of modes given from the input is reached.


9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 47

In this study, we selected 80 modes for vibrational structure calculations, which incorporated 10 target modes and 70 active modes. oc-VSCF51, 65 was first carried out to optimize the coordinates using the cubic force field (CFF) and a pair selection scheme with a threshold value of 250 cm-1. The normal modes with the frequency lower than 500 cm-1 were kept fixed to avoid divergence of the VSCF calculation. Then, the optimized coordinates were used to generate a hybrid potential energy surface (PES)66 up to the three-mode coupling, which combined the quartic force field (QFF) and the grid potential. QFF was calculated by numerical differentiations of the Hessian matrices. The grid potential was calculated for all the one-mode terms and the two- and three mode coupling terms with MCS larger than 1.0 cm-1. 11 and 7 grid points were employed along the coordinate for the one-mode terms and for the two- and three-mode terms, respectively. The dipole moment surfaces (DMS) were constructed at the same time using the dipole moment obtained at the grid points. Finally, the VQDPT2 calculation was carried out with the parameters to construct the degenerate space set to k = 4 and Ngen = 1. The generation of the anharmonic potential and the vibrational structure calculation were performed using SINDO program developed by us. The bottleneck of the vibrational structure calculation was the generation of the anharmonic PES, which required an enormous number of electronic structure calculations to obtain the energy and Hessian. The generation of the hybrid PES for 80 modes required 160 points of Hessian matrices and ~153,000 points of energies to obtain the QFF and the grid PES, respectively. We used the supercomputer facility (Fujitsu Server PRIMERGY CX400 in Nagoya University). The energy and Hessian calculations took ~ 220 sec and 86 min, respectively, using 24 cores in one node. Therefore, the total cost of PES generation was 9,580 node × hour. Nevertheless, the elec-


10

Page 11 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tronic structure calculations can be performed independently, and thus can be massively parallelized.

RESULTS REMD simulation and clustering analysis. The free-energy landscape of a SIVSF-NH2 molecule in vacuum as a function of backbone dihedral angles φ and ψ of each residue is shown in Figure 3, which is calculated by,

E(φ, ψ ) = −kBT ln(

Nφψ N

),

(1)

where kB and T denote the Boltzmann constant and the temperature, respectively, and N is the total number of snapshots (30,000) and Nφψ is the number of snapshots whose backbone dihedral angles belong to grid points of φ and ψ. The grid point was set every 3.6 degree along φ and ψ. At 1300 K, the backbone dihedral angle distribution of all the residues resembles to each other, indicating that the REMD simulation samples a sufficiently wide range in the backbone conformational space. Note that all the residues show a preference to the α-helix region (φ ~ -80 and ψ ~ -30) and the β-sheet region (φ ~ -100 and ψ ~ 120). We have also confirmed that every replica walks over a wide range in the temperature and potential energy space (see Figures S2-S4). The result substantiates that the structure of SIVSF-NH2 is well sampled by the REMD simulation. In contrast, each residue gives different minima at 392 K and 300 K as shown in Figure 3 (b) and (c), respectively. The landscape at 300 K is found to be similar to that at 392 K, suggesting that these minima represent the essential conformation of SIVSF-NH2 in the low temperature. The position of the minima at 300 K matches well with that of γ and β turns.9, 67 The γ turn is a


11


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 47

structure where the CO of the i-th residue and the NH of the (i+2)-th residue are hydrogenbonded, thereby forming a ring structure with 7 atoms. It is characterized in the Ramachandran

Figure 3. The free-energy landscape as a function of backbone dihedral angles, φ (C-N-Cα-C) –

ψ (N-Cα-C-N), of the four residues, Ile, Val, Ser(4), and Phe, obtained from the REMD trajectory


12

Page 13 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


at (a) 1300 K, (b) 392 K, and (c) 300 K. The positions of γ turns (γL and γD) and type I and II’ of β turns (β-I and β-II’) are marked in (c). plot of the (i+1)-th residue by a minimum around (φ , ψ ) = (75, -64) and (-79, 69) referred to as

γD and γL, respectively, which differ in the conformation of the side chain with respect to the backbone. Ile, Val, and Phe clearly show the minimum corresponding to γD and γL. The β turn is defined as a close contact of the i-th and the (i+3)-th residues, in which the Cα atoms are within a distance of 7 Å. It often accompanies a hydrogen bond between the CO of the i-th residue and the NH of the (i+3)-th residue which forms a turn with 10 atoms. Among the nine types of β turns,67 types I and II’ give rise to minima around (φi+1, ψi+1) = (-60, -30) / (φi+2, ψi+2) = (-90, 0) and (φi+1, ψi+1) = (60, -120) / (φi+2, ψi+2) = (-80, 0), respectively, in the Ramachandran plot of the (i+1)-th and (i+2)-th residues. In Figure 3 (c), type I is observed for Ile-Val, Val-Ser(4), and Ser(4)-Phe, while type II’ is also seen for Val-Ser(4). The backbone dihedral angles in Phe are distributed in a relatively wide range of ψ due to the motion of the terminal NH2 that has less restraint compared to other groups.


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 47

Figure 4. The number of clusters as a function of the clustering radius.

We have carried out the k-means clustering analysis using the trajectory obtained at 300 K. The analysis requires a clustering radius as an input parameter to divide the snapshot structures into a group of clusters. As shown in Figure 4, the number of cluster increases steeply upon the decrease of the radius. From the viewpoint of the cost of subsequent electronic and vibrational structure calculations, it is desirable to take a large radius and to keep the number of clusters as small as possible. However, an excessively large radius fails to distinguish different conformers. Therefore, the radius needs to be carefully determined. To this end, we have examined the distribution of the structures obtained from the clustering analysis (the one closest to the center of a cluster) in terms of the dihedral angles (φ, ψ) of the four residues, superimposing it to the free-energy landscape at 300 K. The result is shown in Figure 5. The sampling is obviously insufficient when the radius is set to 2.5 Å missing several minima; for example, (φ, ψ) = (-120, -60) of Phe and (φ, ψ) = (60, -100) of Val. In contrast, the conformers selected with the use of cluster radii of both 1.4 and 1.7 Å cover all the local minima, though the distribution of the former is more dense. Therefore, 144 structures obtained by using a radius of 1.7 Å have been employed for the subsequent calculation.


14

Page 15 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. The distribution of representative structures of the clusters obtained by the k-means clustering analysis using a radius, (a) 1.4 Å, (b) 1.7 Å, and (c) 2.5 Å, in terms of dihedral angles (φ, ψ) of Ile, Val, Ser(4), and Phe, superimposed on the free-energy landscape at 300 K.


15


(a)

80 w/o ZPE

∆ E / kJ mol-1

60

40

20

0 1

21

41

61

81

101

121

141

Conformer ID

(b)

40 with ZPE 30

∆ E / kJ mol-1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 47

w/o ZPE

20

10

0 1

6

11

16

21

Conformer ID

Figure 6. The relative energy of (a) the 144 conformers without the ZPE correction and (b) the 25 lowest energy conformers with and without the ZPE correction.

The energetics based on DFT. Figure 6 shows the relative energy of the conformers obtained by the geometry optimization at the B3LYP/6-31(++)G** level of theory. The rise of the energy with respect to the number of conformer is rather slow beyond ~ 20 kJ mol-1 with a slope of 0.34 kJ mol-1 conformer-1. In other words, there are about three conformers in an energy range of 1 kJ mol-1. This result indicates that the potential of SIVSF-NH2 is rugged and complex with many


16

Page 17 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


local minima. The relative energy with and without the zero-point energy (ZPE) correction is compared in Figure 6 (b) for the 25 lowest energy conformers. Although the energy is slightly different in number, the overall trend remains to be the same before and after the ZPE correction. Therefore, other conformers that lie higher in energy may be safely discarded. The structure of the five lowest-energy conformers is shown in Figure 7 together with the intramolecular HBs defined by,68 rDA < 3.2 Å and θDHA > 130º ,

(2)

where D and A are the donor (either O or N) and acceptor (i.e., O) atoms, and rDA and θDHA denote the distance between D ··· A and the angle of D-H ··· A. The backbone structure is characterized by the γ and β turns as suggested in the Ramachandran plot, Figure 3 (c). Conformers 1 and 3 have a similar backbone structure with two γ turns and one β turn. Conformers 2 and 5 are also somehow similar forming the γ and β turns in the same place, though the latter is more stretched due to the conformation of Ile. Conformer 4 has a sequence of γ turns and appears rather packed. Interestingly, the stability of the conformer hardly correlates with the number of HBs. For example, both Conformers 1 and 4 have seven HBs in total, yet Conformer 1 is more stable than Conformer 4 by as much as 17.9 kJ mol-1. The stability of Conformers 4 and 5 is different by only 0.2 kJ mol-1, even though the latter has less number of HBs. Therefore, not only the HB but also the distortion of the backbone and/or the interaction between the sidechains play an important role to determine the structural stability. The vibrational analysis has been carried out for these five conformers.


17


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

γL turn

γL turn

(a)

Page 18 of 47

β turn

γL turn

β turn

OH-π γL turn

conformer 1 (b)

Ser(1)

Ile γD turn

γD turn

Val

Ser(4)

Phe

β turn

β turn

OH-π

conformer 2

Ser(1)

γL turn

(c)

Ile γL turn

Val

Ser(4)

β turn

Phe γL turn

β turn

OH-π

γL turn

conformer 3

Ser(1)

γD turn γL turn

(d) γL turn

Ile γL turn

Val

Ser(4)

γD turn γL turn

Phe

γD turn

γD turn

conformer 4

Ser(1)

Ile γL turn

(e)

Val

Ser(4)

β turn

β turn

γL turn

Phe

β turn

OH-π

β turn

conformer 5

Ser(1)

Ile

Val

Ser(4)

Phe

Figure 7. The structure (left) and the HB pattern (right) of the five lowest-energy conformers. Hydrogen atoms bound to carbon atoms are not shown for clarity. Color code: white (hydrogen), blue (nitrogen), red (oxygen), gray (carbon), backbone atoms (green).


18

Page 19 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(sym-NH2)F NHV NHI

NHF

Harmonic (asym-NH2)F

OHS1 OHS4

NHS4

oc-VQDPT2 (30 modes)


COF +NHI-b



3200

3300

3400

3500

wavenumber /

3600

3700

3800

cm-1

Figure 8. The theoretical IR spectrum of Conformer 1 computed by the harmonic approximation and the VQDPT2 method using different numbers of active modes. The OH and NH stretching modes of residue R are labeled as OHR and NHR, respectively. The symmetric and asymmetric NH2 stretching modes of Phe are denoted as (sym-NH2)F and (asym-NH2)F, respectively. “COF + NHI-b” is a combination tone of the CO stretching mode of Phe and the NH bending mode of Ile. The spectrum is constructed using a Lorentzian line shape function with the full width at half maximum of 5 cm-1.


19


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 47

Vibrational Analysis. Figure 8 shows the IR spectrum of the lowest-energy conformer, Conformer 1, in the OH/NH stretching region calculated by the harmonic approximation and the VQDPT2 method. Eight peaks are observed in the harmonic spectrum, the two highest-frequency ones being the OH stretching modes of serine and the six others being the NH stretching modes of the backbone. The amide NH stretching modes of Ile, Val, and Phe are found to be the three lowest-frequency peaks, while that of Ser(4) is observed at a position more than 100 cm-1 higher in frequency. This feature is not surprising because the backbone NH in Ile, Val, and Phe forms HBs, while that of Ser(4) remains free [cf. Figure 7 (a)]. The remaining two peaks are assigned to the symmetric and anti-symmetric NH2 stretching modes of Phe. Those of Ser(1) are not visible due to their small IR intensities. In passing, we note that the NH stretching modes are localized to each residue (> 90 % in weight), whereas those of Ile and Val were heavily mixed in the previous report.21 The difference is caused by the level of electronic structure theory employed. The anharmonicity shifts the spectrum substantially lower in frequency by 150 – 200 cm-1 compared to the harmonic one. The amount of shift is dependent on the character of the mode; the ratio between the harmonic and anharmonic frequency, listed in Table 1, varies in a range of 0.94 – 0.96. It is notable that the band positions of asymmetric NH2 stretching mode of Phe and NH stretching mode of Ser(4) are interchanged, since the former is subjected to stronger anharmonicity than the latter. Therefore, the common practice to scale the harmonic frequency with a constant factor leads to a wrong assignment of the band.


20

Page 21 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48


Table 1. The vibrational frequency of NH and OH stretching modes of Conformer 1 (in cm-1) calculated by the harmonic approximation (ω) and VQDPT2 method (ν) using 30, 40, 60, and 80 active modes and the ratio between the harmonic and anharmonic frequency ( f = ν / ω ).

Assignment

Harmonic

VQDPT2 mode = 30

mode = 40

mode = 60

mode = 80

ω

ν

f

ν

f

ν

f

ν

f

NHI

3427.8

3238.2

0.945

3249.5

0.948

3245.6

0.947

3261.9

0.952

NHV

3458.7

3292.1

0.952

3292.5

0.952

3289.9

0.951

3288.7

0.951

NHF

3475.0

3284.9

0.945

3306.6

0.952

3308.9

0.952

3314.3

0.954

(sym-NH2)S1

3498.1

3345.9

0.956

3347.7

0.957

3357.4

0.960

3363.3

0.961

(sym-NH2)F

3498.9

3343.0

0.955

3339.9

0.955

3355.6

0.959

3358.8

0.960

(asym-NH2)S1

3587.0

3407.5

0.950

3410.2

0.951

3429.1

0.956

3446.5

0.961

NHS4

3604.7

3443.2

0.955

3448.8

0.957

3455.1

0.958

3460.5

0.960

(asym-NH2)F

3616.6

3424.3

0.947

3432.8

0.949

3432.7

0.949

3436.1

0.950

OHS1

3682.5

3476.5

0.944

3473.8

0.943

3478.0

0.944

3491.6

0.948

OHS4

3773.5

3584.9

0.950

3583.3

0.950

3577.9

0.948

3581.7

0.949


21


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 47

(a) mode 54

(b) mode 78

(c) mode 184

(d) mode 73

(e) mode 139

(f) mode 156

Figure 9. Three representative vibrational modes of Conformer 1 incorporated in the 30- [(a)(c)] and 80-mode [(d)-(f)] models. See also Tables 2 and 3.

Table 2. Analyses of three vibrational modes selected in the 30-mode model.a

mode

ω

assignment

(cm-1)

otherb

weight of dominant groups group

weight

group

weight

54

430.3

out-of-plane

OHS4

0.90

78

748.4

out-of-plane

NHV

0.43

NHI

0.12

184

1565.2

bending

NHI

0.50

NHF

0.20

group

weight 0.06

NHF

a

The information of all modes are listed in Table S1.

b

Weight of the displacement of atoms not in OH, NH, and NH2 groups.

0.10

0.29 0.26


22

Page 23 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3. Analyses of three vibrational modes selected in the 80-mode model but not in the 60mode model.a

mode

ω

assignment

(cm-1)

otherb

weight of dominant groups group

weight

group

weight

group

weight

NHS4

0.07

73

681.9

-

NHI

0.11

NHV

0.09

139

1276.5

-

NHI

0.23

(NH2)S1

0.07

0.64

156

1380.4

bending

OHS1

0.15

(NH2)S1

0.09

0.75

a

The information of all modes are listed in Table S1.

b

Weight of the displacement of atoms not in OH, NH, and NH2 groups.

0.64

VQDPT2 results show a smooth convergence of the spectrum with respect to the increase in the number of active modes. Three representative modes in the 30- and 80-mode models are visualized in Figure 9. Tables 2 and 3 list the weight of chemical groups in those modes calculated by,

wυ(I ) = ∑

2

∑ ( cυξ ) ,

(3)

i

i∈I ξ =x,y,z

υ

where i runs over atoms in a chemical group I (either OH, NH, or NH2) and ciξ denotes the normal displacement vector of mode υ. Interestingly, the vibrational modes selected in the 30-mode model are localized to OH, NH, and/or NH2 groups. The assignment is readily possible from their visual representation [see Figure 9 (a) – (c)]. In contrast, the modes selected in the 80-mode model have relatively minor contribution of the OH, NH, and NH2 groups, though the total weight amounts to as much as ~ 0.4. As shown in Figure 9 (d) – (f), these modes are not local-


23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 47

ized to any part of chemical groups but delocalized over several residues. The result is consistent with Pele and Gerber,69 who suggested that the mode-mode coupling is strong when the modes share the motion of similar groups of atoms. Also, it has been shown that the potential is more compactly represented in the modes localized to each residue than in normal modes delocalized over several residues.48-49 The present scheme based on MCS systematically incorporates the modes which have important contribution to the OH/NH stretching modes. The VQDPT2 result with 80 active modes shows a good convergence of the spectrum. In the 60- and 80-mode models, the NH stretching mode of Ile (NHI) yields two bands due to Fermi resonance with a combination tone. In the 80-mode model, the vibrational states at 3244.0 and 3261.9 cm-1 are, respectively, represented as,

Ψ1 = 0.77 19511841 + 0.57 2281 ,

(4)

Ψ 2 = 0.66 2281 − 0.56 19511841 + 0.34 19511851 ,

(5)

where m1 and m1n1 denote the VSCF configuration function of a fundamental tone of mode m and a combination tone of mode m and n. The mode numbers, 184, 185, 195, and 228, are the NH bending mode of Ile (NHI-b) and Phe (NHF-b), the CO stretching mode of Phe (COF), and NHI, respectively. The two states shares 2281 , which carries the IR intensity, an thus exhibit a noticeable band in the spectrum. Interestingly, the Fermi splitting is manifested through the vibrational modes that are spatially close; the NH bond of Ile and the CO bond of Phe are hydrogen bonded (see Fig. 7). The resonance is absent in the 30- and 40-mode models because COF is not selected in these models (see Table S1).


24

Page 25 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Exp.

NHI

OHS1

NHV NHF (asym-NH2)F (sym-NH2)F NHS4

Conformer 1 (0.0 kJ mol-1) OHS4

Conformer 2 (sym-NH2)F (8.4 kJ mol-1) NHF (f-NH)S1 OHS1OH

S4 NHI NHS4 (asym-NH2)F

(hb-NH)S1 NHV

Conformer 3 (17.8 kJ mol-1)

(hb-NH)S1 NHF NHV

NHV NHS4

(f-NH)S1 OHS1 NHI NH (f-NH)F OHS4 S4

(hb-NH)F

OHS1 Conformer 4 (17.9 kJ mol-1) NHI (hb-NH)F (sym-NH2)S1 NHF (asym-NH2)S1 (f-NH)

OHS4 F Conformer 5 NHI (sym-NH2)F (18.1 kJ mol-1) NHF NH V OHS1 (asym-NH2)F NHS4

3100

3200

3300

3400

3500

wavenumber / cm-1

OHS4

3600

3700

Figure 10. Theoretical IR spectrum of the five lowest-energy conformers compared with the experimental one. The hydrogen bonded and free NH stretching modes of the terminal NH2 group, in which one of the NH bonds is hydrogen bonded, are labeled “hb-NH” and “f-NH”, respectively. Other labels are the same as in Figure 8. The labels colored in green and blue are the NH stretching modes involved in the γ and β turns, respectively. The red labels are the NH stretching modes of Ser(4) free from HB. The calculated spectrum is constructed using a Lorentzian line shape function with the full width at half maximum of 5 cm-1.


25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 47

Figure 10 shows the IR spectrum calculated for the five lowest-energy conformers together with the experimental one. The spectrum of Conformer 1 is in best agreement with the experiment among the five calculated ones. In Conformers 2 and 3, there are more peaks present than in the experimental spectrum due to visible intensities carried by the NH2 stretching modes of Ser(1). Conformer 4 yields a number of peaks below 3250 cm-1, whereas the peak is absent below 3350 cm-1 in Conformer 5. For a quantitative purpose, the frequencies and IR intensities of Conformer 1 obtained by theory and experiment are compared in Table 4. The experimental frequency and intensity are obtained by fitting the observed spectrum to Lorentz functions. Note that the broad spectrum below 3270 cm-1 fits well to two functions centered at 3272.0 and 3269.5 cm-1 augmented with a broad function centered at 3214.7 cm-1 (shown in Figure S5). The frequencies are highly accurate with mean and maximum absolute deviations of 11.2 cm-1 and 16.7 cm-1 (for NHV), respectively. The ratio of the relative intensity (δR in Table 4) is obtained around 0.5 – 0.6 for the NH stretching bands, indicating that the intensity of OHS4, which is taken as a reference, is overestimated. We note that, if the intensity of OHS4 is multiplied by 0.6, the average of δR improves to 1.1. The intensity of OHS1 is further overestimated by 1.5 times relative to OHS4. The relative intensities between the OH and NH stretching modes may be improved by using larger basis sets and more accurate dipole moment. Nonetheless, the agreement in the intensities is also reasonably accurate. Therefore, SIVSF-NH2 observed in the experiment is unambiguously assigned to Conformer 1 from both energetic and vibrational points of view. A notable discrepancy is found in a region around 3250 cm-1, where the theoretical spectrum shows a sharp Fermi splitting of NHI at 3261.9 and 3244.0 cm-1, while the experiment yields only a single shoulder band at 3269.5 cm-1. The reason is likely due to the accuracy of the PES, since the resonance is sensitive to the level of electronic structure and the functional form. In


26

Page 27 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


fact, the splitting is absent when the dispersion correction to B3LYP is incorporated (see Sec. 4 and Figure 12 for more details). It is also notable that the present scheme selects the active modes by targeting the NH/OH stretching modes, but not the modes that are in resonance (COF and NHI-b). Thus, the unbalanced description between NHI and COF/NHI-b may have caused an accidental degeneracy. In the previous work, two of the authors (SI and MF)21 assigned the observed conformer to Conformer 1 based on a conformational search using CONFLEX and an IR spectrum calculated by the harmonic approximation. The present result reinforces the finding, as the same conclusion is reached by a rather different approach using REMD and VQDPT2.


27


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 47

Table 4. The vibrational frequency (ν in cm-1) and IR intensity (I in km mol-1) of Conformer 1 calculated by VQDPT2 compared with the experiment.

VQDPT2

Exp.

ν

δν

I

Rb

δRc

ν

Rb

OHS4

3581.7

(-6.3)

124.8

1.00

-

3588.0

1.00

OHS1

3491.6

(-14.9)

464.6

3.72

(1.5)

3506.5

2.42

NHS4

3460.5

(14.8)

100.6

0.81

(0.5)

3445.7

1.60

(asym-NH2)F

3436.1

(-3.9)

95.3

0.76

(0.6)

3440.0

1.23

(sym-NH2)F

3358.8

(8.9)

66.7

0.53

(0.3)

3349.9

1.74

NHF

3314.3

(16.6)

280.3

2.25

(0.6)

3297.7

3.67

NHV

3288.7

(16.7)

286.9

2.30

(0.5)

3272.0

4.95

NHI

3261.9

(-7.6)

142.7

1.14

(0.5)

3269.5

2.16

COF+NHI-b

3244.0

-

119.2

0.96

-

-

-

11.2d a

δν = ν(VQDPT2) – ν(Exp.).

b

IR intensities relative to I(OHS4).

c

δR = R(VQDPT2) / R(Exp.).

d

Mean absolute deviations.

e

An average of δR.

0.7e


28

Page 29 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


DISCUSSION The results in the previous section have shown that the proposed method is capable of finding the global minimum and analyzing the vibrational spectrum of SIVSF-NH2. Let us now discuss the generality of the method and possible extensions to larger peptides. In the previous studies, the conformational search has been performed using a conventional search algorithm. For example, the structure of SIVSF-NH2 was sampled using the CONFLEX program with a constraint on a pattern of HBs deduced from the UV and IR spectra.21 Recently, Li et al.70 proposed a method, in which the dihedral angles of each residue (or fragment) are optimized in a step-by-step manner. These methods have been successful in predicting the structure of short peptides. However, the application to large peptides may be limited since the number of minima grows excessively large. On the other hand, REMD has been developed to overcome this issue and has been applied to simulate many biological processes, for instance, protein folding/unfolding, aggregation, and so on. REMD and other generalized-ensemble simulations, which make an efficient search of biomolecular structures feasible, are already available for studying large, complex systems. Nevertheless, the quality of the potential used for the conformational search, i.e., CHARMM36 in this work, is of critical importance to achieve high accuracy. Figure 11 compares the relative energy of the 144 conformers calculated by CHARMM36, B3LYP, B3LYP-D3 [B3LYP with a dispersion correction of Grimme et al.71], and the resolution of identity (RI)-MP2.72 [The RIMP2 calculation was performed using the ORCA program.73] The accuracy of CHARMM36 and B3LYP is comparable with a correlation coefficient of 0.78 and 0.80, respectively. The result may be rather surprising, considering simplified functional forms and numerous approximations


29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 47

∆E / kJ mol-1


∆E(RI-MP2) / kJ mol-1 Figure 11. Plots of the relative energies of the 144 conformers obtained by RI-MP2/def2TZVP(-df) versus those obtained by CHARMM36 (blue circle), B3LYP/6-31(++)G** (orange triangle), and B3LYP-D3/6-31(++)G** (green square). The correlation coefficient is 0.78 and 0.80, and 0.99 for CHARMM36, B3LYP, and B3LYP-D3, respectively. The B3LYP optimized structure is used to evaluate the energy.

made in the parameters of CHARMM36. Although the classical force field has been criticized in the previous works,29-30, 34 the latest CHARMM force field incorporates a two-dimensional grid potential in terms of the dihedral angles of the backbone or the CMAP correction,54, 74 and thus has a significant advantage over the others for a simulation of peptides. However, it is notable in Figure 11 that Conformer 1 is not the global minimum in CHARMM36 but ~ 10 kJ mol-1 higher


30

Page 31 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in energy than another conformer. Therefore, the two step procedure is still necessary, where the candidate structures initially sampled by CHARMM36 are refined by accurate electronic structure calculations. Although B3LYP is successful in predicting the lowest- and the second lowest-energy conformers, the B3LYP energy of other conformers largely deviates from the RI-MP2 energy. For example, the relative energy of Conformer 4 obtained as 16.2 kJ mol-1 in B3LYP is much larger in RI-MP2 as 51.7 kJ mol-1. The drawback of B3LYP is alleviated by the dispersion correction. B3LYP-D3 greatly improves over B3LYP showing an excellent agreement with RI-MP2 with a correlation coefficient of 0.99. The result demonstrates the impact of dispersion effects on energetics. However, the correction may or may not improve the IR spectrum. The same procedure for Conformer 1 using B3LYP-D3 yields an IR spectrum shown in Figure 12. The band positions are overall in good agreement with experiment with a mean absolute deviation of 16.9 cm-1. In particular, the Fermi splitting of NHI is gone in this case. However, it is also observed that the position of OHS1 is too low, overlapping with NHS4 and (asym-NH2)F. OHS1 is deviated by 54.9 cm-1, which is the largest among all bands. A similar statistics has been recently reported by Barone et al.,75 where B3LYP and B3LYP-D performs comparable in anharmonic calculations. The mean absolute error was obtained as 9 and 11 cm-1 for B3LYP and B3LYP-D, respectively, for a series of small to medium size molecules. In our previous studies on small molecules,51, 76 the use of high level electron correlation methods (e.g. coupled-cluster singles and doubles with perturbative triples, CCSD(T)) with large basis sets (at least triple-zeta quality) was required to achieve an accuracy of 10 cm-1. Therefore, it seems that the present result based on B3LYP/631(++)G(d,p) is subject to some lucky cancellations of errors, suggesting that the density functional and basis sets still need to be carefully selected.


31


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

NHI Exp.

(sym-NH2)F (asym-NH2)F NHF

NHV

Page 32 of 47

OHS4

NHS4 OHS1

VQDPT2 (B3LYP-D3)

3100

3200

3300

3400

3500

3600

3700

wavenumber / cm-1 Figure 12. IR spectrum of Conformer 1 calculated by VQDPT2 based on B3LYP-D3 (bottom) compared with the experiment (top).

Next, let us discuss the accuracy of VQDPT2 for predicting the spectrum of peptides in general. To this end, we compare the NH stretching frequency calculated for Conformers 1 – 5 with the experimental data. Chin et al.9 have previously reported a trend of the NH stretching frequencies of the amide backbone and a terminal -CONH2 group according to their HB pattern by compiling the data of short peptides in the gas phase. As shown in Figure 13 (a), the NH bond free from HB was measured around 3480 cm-1, while β and γ turns were red-shifted to a range of 3410 – 3430 cm-1 and 3340 – 3410 cm-1, respectively, due to the formation of HB. The band position showed a further red-shift when the neighboring CO of the amide group accepted a HB.


32

Page 33 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 12 (b) shows a trend of NH stretching frequencies of a -CONH2 group as a function of HB strength defined as,

δν =

(δν NH )

2

− 4A 2 ,

(6)

where δνNH is the difference in a frequency between two NH stretching modes and 2A is the difference between symmetric and anti-symmetric NH2 stretching modes free from HB. In this definition, δν = 0 when the NH2 group is free from HB, and δν ≅ δνNH when one of the NH bonds forms a strong HB and δνNH >> 2A. δν was in a range of 45 – 70 cm-1 and > 75 cm-1 for β and γ turns, respectively. The calculated spectrum is consistent with these data. In Figure 13 (a), the NH stretching frequencies free from HB are calculated in a range of 3450 - 3470 cm-1 (also shown in red in Figure 10) and the frequency becomes lower in an order of β turns, and β and γ turns with a neighboring CO accepting HB. In Figure 13 (b), the calculated values trace the experimental trend quite well and the HB strength is calculated in the right range, i.e., β turns for Conformer 5 and γ turns for Conformers 3 and 4. These results indicate that VQDPT2 is capable of predicting the NH stretching frequency or the amide A bands of polypeptides with high accuracy.


33


(a) HB donor

C O

H N

Exp.

HB donor

H C N O

β

HB acceptor

free NH

turns

γ turns This work free NH

β turns γ turns

NH stretching frequency / cm-1

(b) NH2 stretching frequency / cm-1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 47

2

5

asym-NH2 sym-NH2 2

Exp. This work 3 4

f-NH

5

β turns

hb-NH

3 4

γ turns HB strength (δν ) / cm-1

Figure 13. (a) Plots of the NH stretching frequency of the amide group free from HB (crosses), in β turns (triangles), and in γ turns (circles). Filled triangles and circles have a neighboring CO free from HB (left in the inset), while open ones have the CO accepting HB (right in the inset). (b) Plots of the NH2 stretching frequency of the -CONH2 group as a function of HB strength. The HB strength was derived by Eq. (6) using 2A of Conformer 2 (= 134 cm-1). The number near the plot is the index of isomer. The experimental data are taken from Ref. [9].


34

Page 35 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In Figure 13 (a), there are some notable differences between theory and experiment. First, the band positions of β turns are calculated to be lower than the experimental one by about 50 cm-1. The calculated bands originate from NHF of Conformers 1, 2, 3, and 5. It is notable that these conformers have a motif, where an amide group that accepts a β turn donates to a γ turn (CO in Ile and NH in Val; see Figure 7). Secondly, two bands of γ turns are found at 3230 and 3204 cm-1 extending beyond the low frequency limit of the experimental range (~3260 cm-1). The bands come from NHS4 and NHV of Conformer 4, which has a consecutive γ turns (see Figure 7 (d)). These results indicate that the large red-shift of calculated bands is caused by a sequence of hydrogen-bonded amide groups. The longer sequence of HBs makes the amide groups more polarized, thereby inducing stronger HBs and sizable red-shifts of the NH stretching frequency. Note that these structural motifs are not included in the experimental data set for short peptides. Therefore, further systematic theoretical and experimental studies are needed for characterizing the amide A bands of long peptides. Given that the conformational sampling and the vibrational structure calculation are both feasible with sufficient accuracy for each conformer, an extension to a system under thermal equilibrium, where multiple conformers contribute to the spectrum simultaneously, can be made by taking an average of the spectrum over thermally accessible conformers as,

I(ν ) = ∑σ i I i (ν ) ,

(7)

i

where σi and Ii are the weight and spectrum of the i-th conformer, respectively. Recently, we have calculated the IR and Raman spectra of a lipid bilayer (the amide bands of sphingomyelin) based on Eq. (7).77 For an application to polypeptides, we may use the weight derived from the


35


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 47

REMD simulation and the spectrum calculated by VQDPT2. Furthermore, by modeling the effect of environment using, for example, the QM/MM method, it may be feasible to calculate the vibrational spectrum of peptides in biological environment. In order to carry out the weight average in Eq. (7), it is necessary to compute the spectrum of a number of conformers, which poses a computational challenge. The bottleneck of the present approach is exclusively in the generation of the anharmonic PES. This step is expensive because the number of anharmonic, mode-mode coupling terms increases steeply with respect to the size and DFT calculations need to be performed for each term on grid. Thus, it is of vital importance to identify and pre-screen unimportant, weakly coupled terms to alleviate the unfavorable scaling. In this work, we have proposed a mode-selection scheme based on MCS. Furthermore, an interesting insight gained from Figure 9 is that the relevant mode-mode interaction arises in a spatially short range. An efficient PES generation that exploits the spatial locality is a promising direction, which will be investigated in the future work.


36

Page 37 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


CONCLUSIONS We have performed REMD and VQDPT2 to calculate the vibrational spectrum of a polypeptide, SIVSF-NH2. The conformation of SIVSF-NH2 is widely sampled by the REMD simulation, and the resulting trajectory analyzed by the k-means clustering method provides a good set of candidate structures for the energetic and vibrational analyses. A simple scheme based on MCS is proposed and implemented, which systematically selects relevant vibrational modes for vibrational structure calculation. VQDPT2 with 80 active modes is performed for the five lowestenergy conformers. The comparison of the theoretical and experimental spectrum has unambiguously assigned the observed conformer to the lowest-energy one (Conformer 1). Furthermore, the calculated NH stretching frequencies for the γ- and β-turns are found to be consistent with the experimental data known in the literature. Therefore, the proposed approach is a powerful tool to analyze the vibrational spectrum of polypeptides, which helps reveal the intra- and intermolecular hydrogen bonds and their secondary structures.


37


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 47

ASSOCIATED CONTENT Supporting Information. The Supporting Information is available free of charge on the ACS publication website at DOI: xxxxx. The acceptance ratio of adjacent replicas (Figure S1). The replica number of the trajectory at 300 K as a function of time (Figure S2). The temperature and potential energy of a replica as a function of time (Figure S3). The canonical probability distribution of the potential energy (Figure S4). The fitting of the observed spectrum to 9 Lorentz functions. (Figure S5). The character of the 80 vibrational modes employed in the vibrational structure calculations (Table S1).

AUTHOR INFORMATION Corresponding Author *Phone: +81-48-462-1407. E-mail: [email protected] (Y.S.) Present Addresses #Department of Molecular Microbiology and Immunology, Nagasaki University Graduate School of Biomedical Sciences, 1-12-4 Sakamoto, Nagasaki 852-8523, Japan. Notes The authors declare no competing financial interests. ACKNOWLEDGMENT HO is supported by the Special Postdoctoral Researchers Program at RIKEN. This research is partially supported by the “Molecular Systems”, "iTHES", and "Integrated Lipidology" projects in RIKEN, the Center of innovation Program from Japan Science and Technology Agency, JST,


38

Page 39 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


JSPS KAKENHI Grant No. JP26220807 and JP26119006 (to Y. S.), JSPS KAKENHI Grant No. JP16702840 (to K. Y.), and JP25104008 (to M. F.). We used computational resources provided by the HPCI System Research Project (Project ID: hp140105, hp150022), the RIKEN Integrated Cluster of Clusters (RICC), and MEXT SPIRE Supercomputational Life Science (SCLS).


39


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 47

REFERENCES 1. Barth, A., Infrared Spectroscopy of Proteins. Biochim. Biophys. Acta 2007, 1767, 1073-1101. 2. Hiramatsu, H.; Kitagawa, T., FT-IR Approaches on Amyloid Fibril Structure. Biochim. Biophys. Acta 2005, 1753, 100-107. 3. Ganim, Z.; Chung, H. S.; Smith, A. W.; Deflores, L. P.; Jones, K. C.; Tokmakoff, A., Amide I Two-Dimensional Infrared Spectroscopy of Proteins. Acc. Chem. Res. 2008, 41, 432-441. 4. Middleton, C. T.; Marek, P.; Cao, P.; Chiu, C.-c.; Singh, S.; Woys, A. M.; de Pablo, J. J.; Raleigh, D. P.; Zanni, M. T., Two-Dimensional Infrared Spectroscopy Reveals the Complex Behaviour of an Amyloid Fibril Inhibitor. Nature Chemistry 2012, 4, 355-360. 5. Fu, L.; Liu, J.; Yan, E. C. Y., Chiral Sum Frequency Generation Spectroscopy for Characterizing Protein Secondary Structures at Interfaces. J. Am. Chem. Soc. 2011, 133, 8094-8097. 6. Yan, E. C. Y.; Fu, L.; Wang, Z.; Liu, W., Biological Macromolecules at Interfaces Probed by Chiral Vibrational Sum Frequency Generation Spectroscopy. Chem. Rev. 2014, 114, 84718498. 7. Fu, L.; Wang, Z.; Psciuk, B. T.; Xiao, D.; Batista, V. S.; Yan, E. C. Y., Characterization of Parallel β-Sheets at Interfaces by Chiral Sum Frequency Generation Spectroscopy. J. Phys. Chem. Lett. 2015, 6, 1310-1315. 8. Topics in Current Chemistry: Gas-Phase IR Spectroscopy and Structure of Biological Molecules. Edited by Rijs, A. M.; Jos, O. Springer: 2015; Vol. 364. 9. Chin, W.; Piuzzi, F.; Dimicoli, I.; Mons, M., Probing the Competition between Secondary Structures and Local Preferences in Gas Phase Isolated Peptide Backbones. Phys. Chem. Chem. Phys. 2006, 8, 1033-1048. 10. de Vries, M. S.; Hobza, P., Gas-Phase Spectroscopy of Biomolecular Building Blocks. Annu. Rev. Phys. Chem. 2007, 58, 585-612. 11. Rizzo, T. R.; Stearns, J. A.; Boyarkin, O. V., Spectroscopic Studies of Cold, Gas-Phase Biomolecular Ions. Int. Rev. Phys. Chem. 2009, 28, 481-515. 12. Boyarkin, O. V.; Mercier, S. R.; Kamariotis, A.; Rizzo, T. R., Electronic Spectroscopy of Cold, Protonated Tryptophan and Tyrosine. J. Am. Chem. Soc. 2006, 128, 2816-2817.


40

Page 41 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


13. Fricke, H.; Gerlach, A.; Unterberg, C.; Rzepecki, P.; Schrader, T.; Gerhards, M., Structure of the Tripeptide Model Ac-Val-Tyr(Me)-NHMe and its Cluster with Water Investigated by IR/UV Double Resonance Spectroscopy. Phys. Chem. Chem. Phys. 2004, 6, 4636. 14. Fricke, H.; Schäfer, G.; Schrader, T.; Gerhards, M., Secondary Structure Binding Motifs of the Jet Cooled Tetrapeptide Model Ac–Leu–Val–Tyr(Me)–NHMe. Phys. Chem. Chem. Phys. 2007, 9, 4592. 15. Chakraborty, S.; Yamada, K.; Ishiuchi, S.; Fujii, M., Gas Phase IR spectra of tri-Peptide ZPro-Leu-Gly: Effect of C-Terminal Amide Capping on Secondary Structure. Chem. Phys. Lett. 2012, 531, 41-45. 16. Ishiuchi, S.; Yamada, K.; Chakraborty, S.; Yagi, K.; Fujii, M., Gas-Phase Spectroscopy and Anharmonic Vibrational Analysis of the 3-Residue Peptide Z-Pro-Leu-Gly-NH2 by the Laser Desorption Supersonic Jet Technique. Chem. Phys. 2013, 419, 145-152. 17. Abo-Riziq, A.; Crews, B. O.; Callahan, M. P.; Grace, L.; de Vries, M. S., Spectroscopy of Isolated Gramicidin Peptides. Angew. Chem. Int. Ed. 2006, 45, 5166-5169. 18. Rijs, A. M.; Kabeláč, M.; Abo-Riziq, A.; Hobza, P.; de Vries, M. S., Isolated Gramicidin Peptides Probed by IR Spectroscopy. ChemPhysChem 2011, 12, 1816-1821. 19. Nagornova, N. S.; Guglielmi, M.; Doemer, M.; Tavernelli, I.; Rothlisberger, U.; Rizzo, T. R.; Boyarkin, O. V., Cold-Ion Spectroscopy Reveals the Intrinsic Structure of a Decapeptide. Angew. Chem. Int. Ed. 2011, 50, 5383-5386. 20. Zabuga, A. V.; Rizzo, T. R., Capping Motif for Peptide Helix Formation. J. Phys. Chem. Lett. 2015, 6, 1504-1508. 21. Ishiuchi, S.; Yamada, K.; Oba, H.; Wako, H.; Fujii, M., Gas Phase Ultraviolet and Infrared Spectroscopy on a Partial Peptide of β2-Adrenoceptor SIVSF-NH2 by Laser Desorption Supersonic Jet Technique. Phys. Chem. Chem. Phys. 2016, in press, DOI:10.1039/C6CP04196E. 22. Burke, N. L.; Deblase, A. F.; Redwine, J. G.; Hopkins, J. R.; Mcluckey, S. A.; Zwier, T. S., Gas-Phase Folding of a Prototypical Protonated Pentapeptide: Spectroscopic Evidence for Formation of a Charge-Stabilized β-Hairpin. J. Am. Chem. Soc. 2016, 138, 2849-2857. 23. Fricke, H.; Schwing, K.; Gerlach, A.; Unterberg, C.; Gerhards, M., Investigations of the Water Clusters of the Protected Amino Acid Ac-Phe-OMe by Applying IR/UV Double


41


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 47

Resonance Spectroscopy: Microsolvation of the Backbone. Phys. Chem. Chem. Phys. 2010, 12, 3511-3521. 24. Schwing, K.; Reyheller, C.; Schaly, A.; Kubik, S.; Gerhards, M., Structural Analysis of an Isolated Cyclic Tetrapeptide and its Monohydrate by Combined IR/UV Spectroscopy. ChemPhysChem 2011, 12, 1981-1988. 25. Biswal, H. S.; Loquais, Y.; Tardivel, B.; Gloaguen, E.; Mons, M., Isolated Monohydrates of a Model Peptide Chain: Effect of a First Water Molecule on the Secondary Structure of a Capped Phenylalanine. J. Am. Chem. Soc. 2011, 133, 3931-3942. 26. Cirtog, M.; Rijs, A. M.; Loquais, Y.; Brenner, V.; Tardivel, B.; Gloaguen, E.; Mons, M., Far/Mid-Infrared Signatures of Solvent−Solute Interactions in a Microhydrated Model Peptide Chain. J. Phys. Chem. Lett. 2012, 3, 3307-3311. 27. Nagornova, N. S.; Rizzo, T. R.; Boyarkin, O. V., Interplay of Intra- and Intermolecular HBonding in a Progressively Solvated Macrocyclic Peptide. Science 2012, 336, 320-323. 28. Fricke, H.; Gerlach, A.; Unterberg, C.; Wehner, M.; Schrader, T.; Gerhards, M., Interactions of Small Protected Peptides with Aminopyrazole Derivatives: The Efficiency of Blocking a β-Sheet Model in the Gas Phase. Angew. Chem. Int. Ed. 2009, 48, 900-904. 29. Řeha, D.; Valdés, H.; Vondrášek, J.; Hobza, P.; Abo-Riziq, A.; Crews, B.; de Vries, M. S., Structure and IR Spectrum of Phenylalanyl–Glycyl–Glycine Tripetide in the Gas-Phase: IR/UV Experiments, Ab Initio Quantum Chemical Calculations, and Molecular Dynamic Simulations. Chem. Eur. J. 2005, 11, 6803-6817. 30. Valdes, H.; Spiwok, V.; Rezac, J.; Řeha, D.; Abo-Riziq, A. G.; de Vries, M. S.; Hobza, P., Potential-Energy and Free-Energy Surfaces of Glycyl-Phenylalanyl-Alanine (GFA) Tripeptide: Experiment and Theory. Chem. Eur. J. 2008, 14, 4886-4898. 31. Sugita, Y.; Okamoto, Y., Replica-Exchange Molecular Dynamics Method for Protein Folding. Chem. Phys. Lett. 1999, 314, 141-151. 32. Sugita, Y.; Okamoto, Y., Replica-Exchange Multicanonical Algorithm and Multicanonical Replica-Exchange Method for Simulating Systems with Rough Energy Landscape. Chem. Phys. Lett. 2000, 329, 261-270. 33. Sugita, Y.; Okamoto, Y., Free-Energy Calculations in Protein Folding by GeneralizedEnsemble Algorithms. Lecture Notes in Computational Science and Engineering 2002, 303331.


42

Page 43 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


34. Schubert, F.; Rossi, M.; Baldauf, C.; Pagel, K.; Warnke, S.; Helden, G. v.; Filsinger, F.; Kupser, P.; Meijer, G.; Salwiczek, M., et al., Exploring the Conformational Preferences of 20-Residue Peptides in Isolation: Ac-Ala19-Lys + H+ vs. Ac-Lys-Ala19 + H+ and the Current Reach of DFT. Phys. Chem. Chem. Phys. 2015, 17, 7373-7385. 35. Bouteiller, Y.; Gillet, J.-C.; Grégoire, G.; Schermann, J. P., Transferable Specific Scaling Factors for Interpretation of Infrared Spectra of Biomolecules from Density Functional Theory. J. Phys. Chem. A 2008, 112, 11656-11660. 36. Bouteiller, Y.; Poully, J. C.; Desfrançois, C.; Grégoire, G., Evaluation of MP2, DFT, and DFT-D Methods for the Prediction of Infrared Spectra of Peptides. J. Phys. Chem. A 2009, 113, 6301-6307. 37. Yagi, K.; Karasawa, H.; Hirata, S.; Hirao, K., First-Principles Quantum Calculations on the Infrared Spectrum and Vibrational Dynamics of the Guanine-Cytosine Base Pair. ChemPhysChem 2009, 10, 1442-1444. 38. Bowman, J. M., Self-Consistent Field Energies and Wavefunctions for Coupled Oscillators. J. Chem. Phys. 1978, 68, 608-610. 39. Hansen, M.; Sparta, M.; Seidler, P.; Toffoli, D.; Christiansen, O., New Formulation and Implementation of Vibrational Self-Consistent Field Theory. J. Chem. Theory Comput. 2010, 6, 235-248. 40. Christoffel, K. M.; Bowman, J. M., Investigations of Self-Consistent Field, SCF CI and Virtual State Configuration-Interaction Vibrational Energies for a Model 3-Mode System. Chem. Phys. Lett. 1982, 85, 220-224. 41. Norris, L. S.; Ratner, M. A.; Roitberg, A. E.; Gerber, R. B., Møller-Plesset Perturbation Theory Applied to Vibrational Problems. J. Chem. Phys. 1996, 105, 11261-11267. 42. Christiansen, O., Vibrational Coupled Cluster Theory. J. Chem. Phys. 2004, 120, 2149-2159. 43. Christiansen, O., Response Theory for Vibrational Wave Functions. J. Chem. Phys. 2005, 122, 194105. 44. Neff, M.; Rauhut, G., Toward Large Scale Vibrational Configuration Interaction Calculations. J. Chem. Phys. 2009, 131, 124129. 45. Pfeiffer, F.; Rauhut, G., Multi-Reference Vibration Correlation Methods. J. Chem. Phys. 2014, 140, 064110.


43


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 47

46. Mizukami, W.; Tew, D. P., A Second-Order Multi-Reference Perturbation Method for Molecular Vibrations. J. Chem. Phys. 2013, 139, 194108. 47. Roy, T. K.; Kopysov, V.; Nagornova, N. S.; Rizzo, T. R.; Boyarkin, O. V.; Gerber, R. B., Conformational Structures of a Decapeptide Validated by First Principles Calculations and Cold Ion Spectroscopy. ChemPhysChem 2015, 16, 1374-1378. 48. Panek, P. T.; Jacob, C. R., Efficient Calculation of Anharmonic Vibrational Spectra of Large Molecules with Localized Modes. ChemPhysChem 2014, 15, 3365-3377. 49. Cheng, X.; Steele, R. P., Efficient Anharmonic Vibrational Spectroscopy for Large Molecules Using Local-Mode Coordinates. J. Chem. Phys. 2014, 141, 104105. 50. Yagi, K.; Hirata, S.; Hirao, K., Vibrational Quasi-Degenerate Perturbation Theory: Applications to Fermi Resonance in CO2, H2CO, and C6H6. Phys. Chem. Chem. Phys. 2008, 10, 1781-1788. 51. Yagi, K.; Otaki, H., Vibrational Quasi-Degenerate Perturbation Theory with Optimized Coordinates: Applications to Ethylene and trans-1,3-Butadiene. J. Chem. Phys. 2014, 140, 084113. 52. MacQueen, J., Some Methods for Classification and Analysis of Multivariate Observations. Edited by University of California Press: Berkeley, California, 1967; p 281-297. 53. Feig, M.; Karanicolas, J.; Brooks, C. L., III, MMTSB Tool Set: Enhanced Sampling and Multiscale Modeling Methods for Applications in Structural Biology. J. Mol. Graphics Modell. 2004, 22, 377-395. 54. Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D., Jr., Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone φ, ψ and Side-Chain χ1 and χ2 Dihedral Angles. J. Chem. Theory Comput. 2012, 8, 3257-3273. 55. Brooks, B. R.; Brooks, C. L., III; MacKerell, A. D., Jr; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S., et al., CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30, 1545-1614. 56. Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K., Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781-1802.


44

Page 45 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


57. Becke, A. D., Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648-5652. 58. Lee, C.; Yang, W.; Parr, R. G., Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785-789. 59. Hariharan, P. C.; Pople, J. A., The Influence of Polarization Functions on Molecular Orbital Hydrogenation Energies. Theor. Chim. Acta 1973, 28, 213-222. 60. Hehre, W. J.; Ditchfield, R.; Pople, J. A., Self-Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian-Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules. J. Chem. Phys. 1972, 56, 2257-2261. 61. Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. V. R., Efficient Diffuse Function-Augmented Basis Sets for Anion Calculations. III. The 3-21+G Basis Set for FirstRow Elements, Li-F. J. Comput. Chem. 1983, 4, 294-301. 62. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A., et al., Gaussian 09 Revision D.01. Gaussian, Inc. Wallingford CT, 2009. 63. Yagi, K.; Hirao, K.; Taketsugu, T.; Schmidt, M. W.; Gordon, M. S., Ab Initio Vibrational State Calculations with a Quartic Force Field: Applications to H2CO, C2H4, CH3OH, CH3CCH, and C6H6. J. Chem. Phys. 2004, 121, 1383-1389. 64. Seidler, P.; Kaga, T.; Yagi, K.; Christiansen, O.; Hirao, K., On the Coupling Strength in Potential Energy Surfaces for Vibrational Calculations. Chem. Phys. Lett. 2009, 483, 138142. 65. Yagi, K.; Keçeli, M.; Hirata, S., Optimized Coordinates for Anharmonic Vibrational Structure Theories. J. Chem. Phys. 2012, 137, 204118. 66. Yagi, K.; Hirata, S.; Hirao, K., Multiresolution Potential Energy Surfaces for Vibrational State Calculations. Theor. Chem. Acc. 2007, 118, 681-691. 67. Chou, K.-C., Prediction of Tight Turns and Their Types in Proteins. Analytical Biochemistry 2000, 286, 1-16. 68. Desiraju, G. R.; Steiner, T., The Weak Hydrogen Bond. Edited by Oxford University Press: New York, 2001.


45


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 47

69. Pele, L.; Gerber, R. B., On the Number of Significant Mode-Mode Anharmonic Couplings in Vibrational Calculations: Correlation-Corrected Vibrational Self-Consistent Field Treatment of di-, tri-, and tetrapeptides. J. Chem. Phys. 2008, 128, 165105. 70. Li, H.; Lin, Z.; Luo, Y., A Fragment Based Step-by-Step Strategy for Determining the Most Stable Conformers of Biomolecules. Chem. Phys. Lett. 2014, 610-611, 303-309. 71. Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H., A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements HPu. J. Chem. Phys. 2010, 132, 154104. 72. Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R., RI-MP2: Optimized Auxiliary Basis Sets and Demonstration of Efficiency. Chem. Phys. Lett. 1998, 294, 143-152. 73. Neese, F., The ORCA Program System. WIREs Comput. Mol. Sci. 2012, 2, 73-78. 74. Mackerell, A. D.; Feig, M.; Brooks, C. L., Extending the Treatment of Backbone Energetics in Protein Force Fields: Limitations of Gas-Phase Quantum Mechanics in Reproducing Protein Conformational Distributions in Molecular Dynamics Simulations. J. Comput. Chem. 2004, 25, 1400-1415. 75. Barone, V.; Biczysko, M.; Bloino, J., Fully anharmonic IR and Raman spectra of mediumsize molecular systems: accuracy and interpretation. Phys. Chem. Chem. Phys. 2014, 16, 1759-1787. 76. Keçeli, M.; Shiozaki, T.; Yagi, K.; Hirata, S., Anharmonic vibrational frequencies and vibrationally-averaged structures of key species in hydrocarbon combustion: HCO+, HCO, HNO, HOO, HOO-, CH3+, and CH3. Mol. Phys. 2009, 107, 1283-1301. 77. Yagi, K.; Li, P.-C.; Shirota, K.; Kobayashi, T.; Sugita, Y., A Weight Averaged Approach for Predicting Amide Vibrational Bands of a Sphingomyelin Bilayer. Phys. Chem. Chem. Phys. 2015, 17, 29113-29123.


46

Page 47 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table of Contents Graphic


47

Anharmonic Vibrational Analyses of Pentapeptide Conformations

Recommend Documents