The influence of solvent on conformational properties of peptides with Aib residue—a DFT study

The conformational propensities of the Aib residue on the example of two model peptides Ac-Aib-NHMe (1) and Ac-Aib-NMe2 (2), were studied by B3LYP and M06-2X functionals, in the gas phase and in the polar solvents. To verify the reliability of selected functionals, we also performed MP2 calculations for the tested molecules in vacuum. Polarizable continuum models (PCM and SMD) were used to estimate the solvent effect. Ramachandran maps were calculated to find all energy minima. Noncovalent intramolecular interactions due to hydrogen-bonds and dipole attractions between carbonyl groups are responsible for the relative stabilities of the conformers. In order to verify the theoretical results, the available conformations of similar X-ray structures from the Cambridge Crystallographic Data Center (CCDC) were analyzed. The results of the calculations show that both derivatives with the Aib residue in the gas phase prefer structures stabilized by intramolecular N–H⋯O hydrogen bonds, i.e., C5 and C7 conformations, while polar solvent promotes helical conformation with φ, ψ values equal to +/−60°, +/−40°. In addition, in the case of molecule 2, the helical conformation is the only one available in the polar environment. This result is fully consistent with the X-ray data. Graphical abstract Effect of solvent on the Ramachandran maps of the model peptides with Aib residue Electronic supplementary material The online version of this article (10.1007/s00894-017-3508-4) contains supplementary material, which is available to authorized users.

Since the Aib amino acid is not ribosomally encoded, peptides containing this residue are more resistant to proteolytic enzymes than peptides containing protein amino acids only. The Aib residue is used as a modifier of naturally occurring and biologically active peptides [7,8] due to its unique structural features, introduced by the presence of two methyl groups at C α .
Analysis of peptide crystal structures shows that Aib residues favor the formation of 3 10 -or α-helical structures. The type of helix depends strongly on peptide chain length and on the number of the Aib residues in the peptide. So, it is well recognized that tri-, tetra-and pentapeptides containing at least one Aib residue adopt mainly a 3 10 helix conformation. However, longer (6-20 residues) Aib-containing peptides fold predominantly, but not exclusively, into left-or right-handed α-helices [4,5,9]. The vibrational circular dichroism (VCD) and infre-red (IR) methods are especially reliable for discriminating 3 10 -and α-helices [10]. α-Aminoisobutyric acid Electronic supplementary material The online version of this article (https://doi.org/10.1007/s00894-017-3508-4) contains supplementary material, which is available to authorized users. homooligopeptides in the gas phase and solution were recently studied by Barone and coworkers [11][12][13][14] using an improved AMBER force field. In these studies, the solvent effect was shown as the critical factor governing the conformational behavior of a single Aib residue and in homooligopeptides. Molecular dynamics simulations show that the α-helix is the preferred structure in aqueous solution, while in DMSO the 3 10 -helical structure is predominant.
The α,α-dimethylglycine residue also shows a strong tendency, even stronger than that of proline [15], to promote βturn conformations. For Aib residues, β-turn conformations of type I, I′ and III, III' are usually observed when this nonstandard amino acid residue is placed at both corners of turns. However, occurrence of the Aib residue at the i + 2 position results in a type II β-turn [15][16][17][18][19]. The peptide with the Aib-Gly turn-initiating sequence shows a very stable β-hairpin conformation over a wide temperature range, as studied by isotope-edited IR spectroscopy and molecular modeling [20,21].
The conformational properties of the Aib residue have been extensively studied theoretically. The conformational preferences of a model Ac-Aib-NHMe peptide containing the Aib residue were established for the first time in 1972 [22]. According to the latter authors, the α-aminoisobutyryl residue has a strong tendency to adopt helical conformations, and typical torsion angles φ, ψ for the Aib residue are −57°and −47°, respectively. Subsequent theoretical studies have confirmed these reports. Ramachandran maps calculated using the CFF91 force field indicated that this non-standard amino acid adopts an α-helical conformation in model diamide [23]. However, theoretical studies carried out in the gas phase using quantum-mechanical methods (HF, B3LYP and MP2) showed that the Aib residue has a tendency to adopt C 5 and C 7 conformations stabilized by intramolecular hydrogen bonds [24,25].
Similarly, the potential energy surfaces (PES) calculated by the parm96 force field demonstrated that the most preferred structures of Ac-Aib-NHMe are also C 5 and C 7 conformers [26]. PCM/B3LYP/6-31 + G(d,p) calculations in solvent showed that the most stable structure in a water environment is the extended conformer C 5 , but the energy of the γ turn structure is only 1 kcal mol −1 higher [24].
The conformational properties of Ac-Aib-NMe 2 diamide have not been studied as extensively as their non-methylated C-terminal amide bond analog. The potential energy surfaces were calculated using molecular mechanics methods [23]. These calculations show that the most stable conformation of this peptide is the α conformation with torsional angles of φ, ψ = 60°and 60°, respectively.
As mentioned earlier, the use of non-standard residues such as Aib could be beneficial in enhancing the biological effects of natural or modified peptides. Another promising way to improve the pharmacological parameters of peptides is their modification by replacing the hydrogen atom of the amide bond by a methyl group-referred to as N-methylation. Introduction of a tertiary amide bond into the peptide chain results in a reduction in conformational freedom of the peptide due to steric hindrance [27,28]. In peptides modified in this way, the tendency to adop a cis-configuration of the amide bond is frequently observed [29][30][31]. N-Methylation is a powerful means of increasing the proteolytic stability [32,33], membrane permeability (lipophilicity) [34,35], and bioavailability [36] of natural peptides. There are several examples of modified, N-methylated peptides that exhibit much better pharmacokinetic properties [37][38][39]. Moreover, a few Nmethylated peptides are currently being evaluated in clinical trials, displaying the promise of N-methylation in delivering next generation drugs [40].
The conformational preferences of peptides depend on a delicate balance between intramolecular interactions and the impact of the environment. Although H-bonded interactions are mainly electrostatic in nature, the contribution of dispersion forces in computing accurate interaction energies is not negligible. The hybrid functional B3LYP does not describe dispersion forces correctly [41], and often underestimates the energy of the hydrogen bond [38,39]. The meta-GGA M06-2X dispersion corrected functional gives better results as regards the energy and geometry of hydrogen bonds [42][43][44][45]69]; however, in some cases, it overestimates the interaction energies and predicts unreasonable structures of N-H⋯O hydrogen bonds in peptides [46]. In this study, we wanted to see how both these methods model the conformational properties of Aib residue derivatives, where we consider dispersive interactions to play a particularly important role.
In this report, we present the results of density functional theory (DFT) calculations on the conformational properties of two model peptides with Aib residues: Ac-Aib-NHMe (1) and Ac-Aib-NMe 2 (2). Calculations were performed at the M06-2X/6-31++G(d,p) and B3LYP/6-31++G(d,p) levels of theory in the gas phase, chloroform and water, where the effect of solvent was included using polarized continuum methods (PCM or SMD). We are mainly interested in interactions stabilizing the minima of these compounds, and how the solvent affects their conformational preferences. Another problem analyzed in this paper is the impact of N-methylation on the ability of peptides containing Aib residues to adopt typical secondary structure motifs.

Methods
Theoretical calculations of conformational properties were carried for two model peptides: Ac-Aib-NHMe (1) and Ac-Aib-NMe 2 (2) with the trans N-terminal amide group (ω 01 80°). The chemical structures and torsional parametres of the studied models are defined in Fig. 1. All ab initio and DFT calculations were performed using the Gaussian 09 package [47]. The structural preferences of Ac-Aib-NHMe and Ac-Aib-NMe 2 were determined by Ramachandran maps E(ϕ,ψ) showing the dependence of potential energy on torsional angles φ and ψ.
Due to the symmetry of the studied molecules, the entire conformational map could be reproduced by calculating only half of the grid points [because E(ϕ,ψ) = E(−ϕ,-ψ)]. The ϕ, ψ PES of each molecule was generated on the basis of 84 structures, partly optimized at the B3LYP/6-31++G(d,p) level. In each grid point, the geometrical parameters were fully relaxed, except for the constrained torsion angles ϕ and ψ. The values of these angles varied from −180°to 0°, and from −180°to 180°for ϕ and ψ, respectively, and the step size were 30°. The energy surface was created using the Surfer 8 program with the radial basis function as a gridding method [48]. Single-point MP2 and M06-2X calculations with 6-31++ G(d,p) basis set were performed on partly B3LYP-optimized  structures. To estimate the effects of enviroment on the topology of the energy surfaces, single point calculations were conducted for each grid point using two continuum solvent models: PCM [49,50] and SMD [51].
Ramachandran plots are traditionally used as a convenient way to present the conformational properties of small peptide model, and they are accessible by accurate experimental and computational approaches [23,24,26,30,52,53].
All low energy areas of the conformational maps were analyzed, and the minima found were re-optimized using B3LYP, MP2, M06-2X methods with the 6-31++G(d,p) basis set in vacuo. Because we are aware of potential problems of the MP2 approach combined with finite basis sets when applied to conformers of polypeptides [54], for comparison we recalculated the minima using 6-311++G(3df,2pd) basis set. Next, a full geometry optimization in chloroform and water was performed using the PCM and SMD models by B3LYP/ 6-31++G(d,p) and M06-2X/6-31++G(d,p) methods.
For each conformer, we performed vibrational analysis to check the absence of imaginary freuqencies. The abundances p of individual conformers were estimated on the basis of the relative energies [55]. The regions of the Ramachandran map are employed as conformational descriptors for backbone orientations of peptides and are labeled differently by different research groups. In this paper, the energyminimized conformers of the investigated molecules are described by the general short-hand letter notation introduced by Zimmerman [56]. Figure 2 presents the ΔE = f (φ, ψ) PESs for the studied molecules 1 and 2, respectively, calculated using the M06-2X method in vacuum and water modeled with the PCM method. The analogous conformational maps calculated in chloroform are presented as Figs. S1 and S2 in the supplementary materials. Also the ΔE(ϕ, ψ) PESs calculated by the MP2 and B3LYP methods in the gas phase and in solvent environments are shown in Figs. S3 and S4 in the supplementary materials. On each map, the local minima are depicted with their Zimmerman notation. Conformations C, E, A, F and D are equivalent to the γ-turn (C 7 ), extended (C 5 ), α-helical, polyproline-like (P II ) and β 2 structures in the literature, respectively. Table 1 lists the backbone torsion angles (ϕ,ψ), relative Table 1 Selected torsion angles (°), relative energies ΔE (kcal mol −1 ) and theoretical abundances p (%) of local minima for (1) and (2) in vacuo  Tables 2 and 3, respectively. The presence of achiral α-carbon of the studied molecules results in symmetry of their maps with respect to the point (φ, ψ = 0°, 0°). Therefore, in a discussion of the obtained results, only the minima found in the left halves of the maps have been taken into consideration. A detailed conformational analysis of the conformers of studied molecules was performed on an assumption that hydrogen bonds (N-H⋯O, N-H⋯N, C-H⋯O) and dipole-dipole attractions between carbonyl groups are the main stabilizing internal forces. Tables S1 and S2 in supplementary material collect structural parameters of the X-H⋯A interactions and dipole attractions based on Steiner's [57] and Allen's [58] criteria, respectively. Figures 3 and 4 show the local minima with stabilizing interactions for the two studied compounds in the gas phase.

Results and discussion
Ac-Aib-NHMe (1) As shown in Fig. 2, in the gas phase and in water, the M06-2X conformational maps of Ac-Aib-NHMe reveal four lowenergy conformers-C, E, F, A-and their mirror images. Moreover, a high energy minimum D (in the gas phase ΔE ≈ 4 kcal mol −1 as predicted by M06-2X method) also occurs on Ramachandran maps of molecule 1, which, due to its energy, probably has no practical significance. The conformational maps obtained with B3LYP and MP2 methods (Fig.  S3in supplementary material) look very similar in terms of the shape, number and position of the minima. This probably means that the conformational properties of this compound are determined primarily by steric interactions and that only two areas of the map are accessible: an area of extended conformations with minimum E, and a second area closer to the center of the map, with remaining three low energy minima.
Inspection of the results for molecule 1 listed in Table 1 shows that results obtained by all computational methods used are quite similar. All methods predict that, in a vacuum, both Table 2 Selected torsion angles (°), relative energies ΔE (kcal mol −1 ) and theoretical abundances p (%) of local minima for Ac-Aib-NHNe (1)  conformers stabilized by N-H⋯O hydrogen bonds (C and E) have low energies (ΔE < 1.4 kcal mol −1 ) and their theoretically estimated abundance is in the range 85-98%. According to the results of all methods, the global minimum of 1 is the conformer C stabilized by relatively short C 7 N 2 -H⋯O 1 hydrogen bond ( H⋯O distance =1.93 Å) (see Table S1 in supplementary material) and one weak C β -H⋯O 1 interaction (Fig. 3). Moreover, dipole attraction between the two carbonyl groups also stabilizes this conformation. The second low energy conformer, E, shows M06-2X relative energy of about 0.4 kcal mol −1 higher. This conformer is a fully extended structure with torsion angles φ, ψ = −180°, −179°, and is stabilized mainly by the C 5 hydrogen bond N 1 -H⋯O 2 . Additionally, two hydrogen bonds C β -H⋯O 1 seem to play role in the stabilization of this conformation. The relative energies of conformer E calculated by M06-2X and B3LYP methods are essentially identical. However, the relative energy of this structure obtained with MP2 method is much higher (ΔE = 1.31 kcal mol −1 ). Overstated energy for extended structures obtained by teh MP2 method with this double zeta basis set was also reported for Ac-Gly-Phe-NH 2 , Ac-Gly-ΔPhe-NHMe and Ac-Gly-ΔPhe-NMe 2 dipeptides [59,60]. With the use of a larger triple zeta basis set [6-311++G(3df,2pd)] within the MP2 method, this effect was no longer observed.
Moreover, it is worth noting that both DFT methods estimate the stability of helical conformer A differently. This conformer (φ, ψ = −65°, −26°) is non N-H⋯O hydrogen-bonded. The main stabilizing force is a weak N 2 -H⋯N 1 interaction, which is present only in this conformer and C β -H⋯O 1 contact (see Table S1 in supplementary material). It also has a short antiparallel dipole-dipole attraction (see Table S2). It is interesting that the B3LYP results indicate the high energy of conformer A in the gas phase (ΔE ≈ 2.6 kcal mol −1 ). The results obtained using the M06-2X method (ΔE ≈ 1.8 kcal mol −1 ) seem more reliable because MP2/6-311++G(3df,2pd) predicts similar relative energies of the obtained conformers. Furthermore, this is also consistent with the experimentally demonstrated ability of the Aib residue to induce an α-helix conformation of the peptide [4,5,9].
A systematic exploration of the potential energy surface and full geometry optimizations for all local minima obtained in the gas phase were also carried out in chloroform and water Table 3 Selected torsion angles (°), relative energies ΔE (kcal mol −1 ) and theoretical abundances p (%) of local minima for Ac-Aib-NMe 2 (2)  using M06-2X and B3LYP methods, combined with PCM and SMD models ( Table 2). Calculations predict that the solvent effect on the conformational properties of Ac-Aib-NHMe is rather limited. As already mentioned, the topology of maps calculated in a vacuum and taking into account the effect of the solvent are very similar. All local minima found in the gas phase are also present in chloroform and water, and conformer D always has a high energy (ΔE > 3 kcal mol −1 ). However, the transition from the gas phase to chloroform and water causes some shifts in the backbone torsion angles φ, ψ for the local minima. In particular, the largest changes were observed for conformations A and F, the backbone dihedral angles of which correspond to the α-helix and polyproline II helix (PII) conformations, respectively. For example, on going from the gas phase to water, the M06-2X method combined with PCM model provides shifts of −10.2°and +9.6°in ψ of the conformation A and F, respectively. The corresponding values are −14.1°and +13.6°at the SMD/B3LYP/6-31++ G(d,p) level. Chloroform and water significantly reduce the energy differences between the minima. For example, in chloroform the gap energy between the first and the fourth conformer equals 0.43 kcal mol −1 in the case of geometries calculated at M06-2X level with SMD model of solvent, and 0.79 kcal mol −1 with PCM. Furthermore, solvent environment rearranges the order of low-energy minima. Besides, for each method, different global minima were obtained, and their energies differ less than the error of the DFT methods used. However, a general regularity may be noted: both DFT methods and both solvent models predict a significant stabilization of helical conformer A. Calculations with PCM model predict that conformer A, with torsion angles φ, ψ = −60°, −36°, becomes the global minimum in water. Conformer A is also the most preferred structure of the studied molecule in chloroform at the SMD M06-2X/6-31++G(d,p) level of theory. Another conformer strongly stabilized by the solvent is structure F. Thus, the results obtained by the SMD method predict that conformer F is the lowest-energy structure in water.
It is worth noting that these two conformations discussed above (A and F), are the most strongly influenced by the polar environment, and are stabilized mainly by the short and strong dipole C = O attractions between the carbonyls of the amide groups (see Table S2). The theoretically estimated abundances of these two conformers in water using PCM/ M06-2X and SMD/M06-2X models are 86% and 97%, respectively. The above-described effect of solvent on conformational properties of Ac-Aib-NHMe is very similar to that obtained for the most studied Ac-Ala-NHMe model peptide. Various ab initio and DFT calculations indicated that solvent effects stabilize the conformer corresponding to the α-helix secondary structure, and flatten the PES [61][62][63]. The results of explicit water calculations also show that hydration of the peptide backbone critically depends on the backbone conformation, and allowed us to determine that the N-H⋯O H-bond formed by a dipeptide in its extended conformation is weakened by the close proximity of the O atom of the neighboring peptide group to the NH proton donor [46,[63][64][65][66][67][68].
Ac-Aib-NMe 2 (2) Table 1 shows the backbone torsion angles and relative energies of the local minima obtained for Ac-Aib-NMe 2 . The conformational maps of this molecule in the gas phase, chloroform and water reveal five local minima-A, E, F, D 2 , D-and their mirror images (Fig. 2) regardless of calculation method (the conformational maps obtained with B3LYP and MP2 methods are presented at Fig. S4 in supplementary material). The presence of an additional methyl group at the C-terminal amide bond results in increasing the energy of potential surface center, which corresponds to conformation with φ, ψ = 0°,0°torsion angles. This is caused by a steric repulsion between the Cterminal methyl group (−C 5 H 3 ) and the oxygen atom of the C = O 1 group, and instead of the C 7 conformation observed for Ac-Aib-NHMe (1), there are two high-energy local D minima stabilized by two weak C-H⋯O contacts.
In the gas phase, the most preferred conformation found by both DFT (M06-2X, B3LYP) and MP2 methods, is the fully extended structure situated in region E (φ, ψ torsion angles are almost 180°). In this conformation, the short C 5 N 1 -H⋯O 2 hydrogen bond is the main stabilizing factor, and the H⋯O 2 distance, according to the M06-2X functional, is only 1.76 Å. Moreover, there are three C-H⋯O contacts: C 5 -H⋯O 2 , C β1 -H⋯O 1 , C β2 -H⋯O 1 , which additionally stabilize this global minimum (Fig. 4).
All calculations predict that the second in energy order is conformer A, stabilized mainly by dipole-dipole interaction between C = O groups and by one weak C 4 -H⋯O 2 hydrogen bond. Its relative energy depends strongly on the method of calculation. As in the case of molecule (1), B3LYP density functional distinctly overestimates energy of this helical conformation. A more reliable result is obtained with M06-2X (ΔE = 1.05 kcal mol −1 ) because this value is closer to that obtained using the MP2/6-311++G(3df,2pd) method. The conformational preferences of Ac-Aib-NMe 2 (2) were examined also in chloroform and water (Table 3, Fig. 1). For neither solvent, did we observe a significant effect on the conformation of the compound. The shape of the maps and the number and location of the minima are the same as in the gas phase. However, in all cases, we observe a clear stabilization of the helical conformation A. Regardless of the calculation method and model of solvent, the two lowest-energy structures are conformers A and E, and their summary theoretically estimated abundance is 96-98%. Increased solvent polarity stabilizes the helical conformation to a greater extent. All our results indicate that, in water, this conformer is a global minimum, which remains in accordance with experimental data. Figure 5 shows the energy of the interaction with water (estimated within the implicit solvent model PCM) as a function of the torsion angles φ, ψ for both studied molecules. The brightest areas on the maps indicate conformations of molecules 1 and 2 with the most significant solvent influence. As can be seen, for both compounds, the highest solvent stabilization energies are observed for helical conformations. The (1) (2) energy maxima occur at (ϕ,ψ = (−60°, −100°) and (−70°, −110°) for 1 and 2, respectively. These results are similar to those obtained for Ac-Ala-NHMe with a trans N-terminal amide bond [61,68]. The presence of an additional methyl group at the C-terminal part in diamide 2 causes slightly weaker interaction with the solvent, by about 1 kcal mol −1 . Additionally, the results of calculations show that there is a large difference in the energy values obtained by the PCM and SMD models. Within the SMD model, energy values are 5 kcal mol −1 and 6 kcal mol −1 higher for chloroform and water, respectively. Despite these differences, we can conclude that the solvent stabilizes helical conformations of the studied Aib derivatives, and that this is closely related to the dipole moment of those compounds [30].

X-ray structures of Aib derivatives
To verify the obtained theoretical results, the conformations of the peptides containing the Aib residue, gathered in the Cambridge Crystallographic Data Center (CCDC), were analyzed. The database search yielded 1116 peptide structures with the Aib residue in the vicinity of the secondary amide bond, and 68 structures with a tertiary amide bond. On the calculated map of molecules 1 and 2 in water the φ, ψ torsion angles values corresponding to the conformers found in the crystal state of Aib derivatives are marked with red crosses (Fig. 6). It is apparent from the maps presented that practically all the X-ray structures from CCDC are in regions of the calculated minima. For molecule 1, the vast majority of them are structures, where the Aib residues adopt the right or lefthanded helix conformation, 69% and 29%, respectively, while 1.7% of all crystal structures correspond to structures F. For Aib residues in the vicinity of the tertiary amide bond, all crystal structures found are within the helical conformation (62% right-handed helix, 38% left-handed helix). On each map (Fig. 6) the green bold line defines the area of conformations with relative energies < 4 kcal mol −1 . Determined in this way, available areas of Ramachandran diagram for diamides 1 and 2 are 15% and 4%, respectively. This means that the tertiary amide group at the C-terminal side of Aib residue significantly reduces the conformational freedom of the compound.

Conclusions
The results of the theoretical calculations presented in this paper highlight the effect of chloroform and water on the conformational properties of model peptides with Aib residue. The conformational preferences of two Aib derivatives Ac-Aib-NHMe (1) and Ac-Aib-NMe 2 (2) have been explored by the M06-2X/6-31++G(d,p) and B3LYP/6-31++G(d,p) methods in the gas phase and in a solvent environment, and by ab initio calculations at the MP2/6-311++G(3df,2pd) level in the gas phase. The results obtained show that both studied model peptides in the gas phase adopt structures stabilized by N-H⋯O hydrogen bonds, i.e., C 5 or C 7 conformations. However, in the polar environment, helical conformations φ, ψ = (+/− 60°, +/− 40°) are the most stable, especially in the case of Aib residues in the vicinity of the tertiary amide. As a result, in the case of molecule 2, the helical conformation is the only one available in the polar solvent. These conclusions fully agree with the crystallographic data. The CCDC data shows that 98% of structures -Aib-NH-and 100% of -Aib-NMe-are in helical conformation.