Microhydration of Protonated Nα-Acetylhistidine: A Theoretical

Aug 7, 2015 - In contrast, the M06/6-311++G(d,p) method appears to slightly overestimate ΔhH° and ΔhG° (MAD = 7.6 kJ mol–1). When looking at the...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCB

Microhydration of Protonated Nα-Acetylhistidine: A Theoretical Approach Vanessa Riffet,† Guy Bouchoux, and Gilles Frison* Laboratoire de Chimie Moléculaire, Ecole Polytechnique and CNRS, 91128 Palaiseau cedex, France S Supporting Information *

ABSTRACT: Extensive exploration of the potential energy surfaces of protonated Nα-acetylhistidine hydrated by 0−3 molecules of water was performed. The methodology combined hierarchical and genealogical (Darwin family tree) approaches using polarizable AMOEBA force field and M06 functional. It is demonstrated that this mixed approach allows recovering a larger number of conformers than the number recovered by using any one of the two methods alone. Hydration enthalpies of protonated Nα-acetylhistidine and of model compounds have been computed using higher theoretical methods, up to the G4MP2 procedure. Excellent agreement with experiment is observed for successive hydration of methylamonium and imidazolium cations using MP2/6-311++G(2d,2p)// M06/6-311++G(d,p) and G4MP2 methods, thereby validating the theory levels used for hydrated protonated Nαacetylhistidine. It is found that the first hydration enthalpy of protonated Nα-acetylhistidine is ca. 10 kJ mol−1 lower than that of imidazolium, a result explained by the local environment of the positively charged imidazolium moiety.



INTRODUCTION Understanding the tridimensional structure of proteins in bulk water is far from being achieved. There is no doubt that protein folding is intimately influenced by intra- and intermolecular interactions, among which are hydrophobic forces and hydrogen bonds.1,2 In this context, studies of solvation by a discrete number of water molecules, i.e., microhydration, allow bridging the gap between isolated and fully solvated molecules. This strategy has been used in a number of experimental and theoretical studies for solutes ranging from inorganic ions to protonated biomolecules such as amino acids or peptides. For instance, complexation thermochemistry for hydration by a discrete number of water molecules of protonated amino acids and peptides has been determined by several mass spectrometry techniques. Indeed, thermal equilibrium, A+(H2O)n−1 + H2O → A+(H2O)n, has been performed into high pressure ion sources3−8 or into water pressurized drift cells9−13 thus allowing the corresponding Gibbs free energy change to be determined. In addition, measuring the equilibrium constants as a function of temperature allowed the determination of the relevant enthalpy and entropy changes via a van’t Hoff type analysis. Blackbody infrared radiative dissociation (BIRD) may provide thermochemical parameters via numerical modeling of dissociation rates of A+(H2O)n complexes.14 Similarly, complexation energies may be deduced from modeling of dissociation efficiencies originating from guided ion beam experiments.15 Beside these experimental source of thermochemical data, infrared multiphoton dissociation (IRMPD) coupled with mass spectrometry brings direct information on the structure of the complexes A+(H2O)n due to typical vibrational signatures of CO, OH, and NH stretching modes.16−18 Comprehensive insight into the solvation process is however possible only if © 2015 American Chemical Society

experimental information is coupled with computation of the corresponding molecular structures. From this point of view, difficulties arise when several sites for hydration are in competition thus leading to an increasing number of isomers and conformers in the hydrated ion population. Therefore, the computational approach used to explore the lowest energy part of the potential energy surface of the system should be both accurate and exhaustive. In the recent years, two theoretical approaches have been proven successful to explore the potential energy surface of microhydrated biological molecules and ions. The first strategy is a hierarchical approach, where as a first step, a large number of starting geometries are generated by molecular dynamics (MD) and/or Monte Carlo (MC) calculations using force field theories. Low energy structures are then reoptimized using quantum chemistry tools at a high, state-of-the-art, level, either directly19,20 or after several intermediate steps21−24 for final selection of the more reliable structures. One crucial point in this strategy is the ability of force fields, and more generally the used intermediate levels, to explore satisfyingly the potential energy surface, both in terms of structures and relative energies. Jensen et al. have shown that for neutral acetyl and N-methyl amino acids25 conventional force fields are of limited accuracy, whereas the AMOEBA polarizable force field26,27 gives significantly more reliable structures and energies. Relevance of the AMOEBA force field, due to an efficient treatment of electrostatic and polarization effects, has been confirmed for charged systems such as Received: June 11, 2015 Revised: August 7, 2015 Published: August 7, 2015 11527

DOI: 10.1021/acs.jpcb.5b05581 J. Phys. Chem. B 2015, 119, 11527−11539

Article

The Journal of Physical Chemistry B protonated and deprotonated methionine,28 protonated Gramicidin S29 and S-nitrosoglutathione30 peptides, or sodiated oligoglycines.31 The second strategy, based on the creation of a Darwinian family tree, considers that each generation of conformations Gn of a molecule microsolvated by n water molecules is constructed from the previous generation Gn−1. This strategy has been developed32 and successfully used for small protonated amino acids such as Gly,33 Ala,34 and Pro.35 Larger systems like deprotonated Asp and Gln behave in the same way,36 as well as other molecular systems like cation-phenol.37 However, recent results may question such an approach when applied to larger systems. For example, if a recent experimental study demonstrates that the linear conformation of diprotonated 1,7-diaminoheptane is kept under hydration up to 10 water molecules, adding more than 12 water molecules induces a more folded conformation.38 This means that the most stable conformers of generation G12 could not be related to that of G10.38 Other examples where hydration induces change (i) in the protonation site,39,40 (ii) between nonzwitterionic and zwiterrionic structures,41−44 or (iii) in the conformational structures45 have been evidenced. In such circumstances, the Darwin tree computational strategy remains to be evaluated. Furthermore, this strategy relies on the knowledge of the first generation G0 which may also be a challenging task for flexible molecules. In the present work, on the basis of the idea that these two strategies turn to be complementary, we define a new way to match these two approaches without significantly increasing the computational time. The two strategies have not been followed in parallel, as done previously,30,46 but sequentially. First, the hierarchical approach is used for each generation to identify the most stable conformers. Then the Darwin tree is built on the basis of all the minima identified. This should prevent gaps in the hierarchical strategy due to the somehow imprecision of the used computational methods. Conversely, the observation of a new conformer at the nth generation, which is not related to any parents, allows including a new family to the tree which would have been missed only following the Darwin tree strategy. Hydrated clusters exhibit various structural patterns characterized by cooperative hydrogen bonding or (and) van der Waals type interactions. These weakly bound species can be accurately described by quantum chemical methods using either high level ab initio molecular orbital theory, or ad hoc densityfunctional theoretical (DFT) models. The latter methods are obviously the best choice in terms of computational efficiency. Recently, versatile classes of functionals that account for dispersive interactions have been developed, which included (i) the nonlocal van der Waals density functionals,47−49 (ii) the local atomic potential approach,50 (iii) semiclassical51−53 or density dependent54 dispersion correction functionals, and (iv) parametrized functionals,55 among them the “meta-hybrid GGA” M06 functional56 used throughout this study. N-Acetylated amino acids are simple model systems that allow investigating the microhydration process of molecules containing a peptide bond. In the gas phase, most of the isolated N-acetylated amino acids protonate on the carbonyl of the acyl moiety.57,58 Exceptions to this rule are amino acids containing a basic side chain such as Arg, Lys, and His.57−59 The present study is focused on microhydration of protonated Nα-acetylhistidine, noted [NAH+H]+. Recently, a detailed study of histidine itself demonstrates the favored protonation

on the imidazole group and the participation of several types of internal hydrogen bonds in the stabilization of neutral, protonated, and deprotonated structures.60 In protonated Nαacetylhistidine, new possibilities of internal hydrogen bond are obviously opened. Furthermore, this small model peptide presents several sites able to interact with water molecules, as summarized in Figure 1. Indeed, this system includes (i) two

Figure 1. Nomenclature used throughout this study to describe hydrogen bond donor sites (a−d, red arrows) and acceptor sites (e−h, blue arrows) of protonated Nα-acetylhistidine, [NAH+H]+.

cationic imidazolium hydrogen bond donor sites, NδH (site a) and NεH (site b), (ii) two neutral hydrogen bond donor sites, that is the amide NH (site c) and the carboxylic OH (site d) groups, (iii) two carbonyl oxygen atoms known to be efficient hydrogen bond acceptor sites (sites e and f), and (iv) the hydroxylic oxygen and the nitrogen of the acylamide CH3CONH group which, to a lesser extent, may be also considered as acceptor sites (sites g and h). Note that sites d and e in a syn-HOCO arrangement of the acidic function generally act cooperatively with respect to a water molecule. These hydrogen bond donor and acceptor sites have been quantitatively confirmed by calculation of their NBO charges (vide infra). Last, a water molecule can also insert in an existing hydrogen bond (indicated as i-type conformers) or interact with another water molecule (indicated as w-type conformers). The competition between fixation of a water molecule at one of the active sites presented in Figure 1 or its possible insertion into an intramolecular H-bond is a question explored in this study. Last, we propose a set of evaluated theoretical gas-phase hydration thermochemical parameters based on the most stable identified [NAH+H]+(H2O)n (n = 0−3) species.



METHODOLOGY As indicated in the Introduction, the strategy adopted to explore exhaustively the conformational landscapes of [NAH +H]+(H2O)n (n = 0−3) complexes combines the hierarchical and Darwin tree approaches. (1) For the hierarchical approach, we have used molecular dynamic (MD) and Monte Carlo (MC) simulations together with basin-hopping conformational scanning algorithm as implemented in the Tinker molecular modeling package61 associated with the AMOEBA polarizable force field. The parameters for the neutral COOH termination were adjusted on the basis of those of the carboxylic group of the aspartic acid, which are already available in the AMOEBA parameter set.46 For a given starting structure, MD calculations have been conducted at several different temperatures (300, 500, 750, and 1000 K), and overall dynamic times of 0.16 and 4.0 ns during which structures were sampled every 0.16 and 4.0 ps, respectively, at the AMOEBA level. Considering the basinhopping procedure, the conformational search was monitored 11528

DOI: 10.1021/acs.jpcb.5b05581 J. Phys. Chem. B 2015, 119, 11527−11539

Article

The Journal of Physical Chemistry B

Figure 2. Optimized geometries of the five most stable conformers A, B, C, D, and E of protonated Nα-acetylhistidine [NAH+H]+ (system G0).

via torsional motions using an energy window of 418 kJ mol−1. This exploration has been fully achieved for [NAH+H]+(H2O)n (n = 0−3) systems to generate a first set of structures. This procedure led to 143, 470, 1076, and 3358 minima for n = 0, 1, 2, and 3, respectively. Then, all structures situated in the 0−20 kJ mol−1 (0−60 kJ mol−1 for the first generation G0) potential energy range relative to the most stable conformers given by AMOEBA were optimized at the M06/6-31+G(d,p), thereby generating 59, 37, 44, and 113 conformers for n = 0, 1, 2, and 3, respectively. Further, the low lying structures were reoptimized at the M06/6-311++G(d,p) level. Vibrational frequencies were calculated at this same DFT level to obtain enthalpies and Gibbs free energies at 298 K. (2) All the structures obtained at the M06/6-311++G(d,p) level for each generation G0, G1, G2, and G3 were then used to build a Darwinian tree. Accordingly, examination of the resulting tree allowed us to identify, in a given family, structures that were not generated by the abovedescribed hierarchical approach. These missing structures are built manually, by identifying hydrogen bond donor and acceptor groups available for water molecule fixation, and optimized at the M06/6-311++G(d,p) level. For example, no isomers were found in G2 for the B family after the hierarchical approach although G1 includes Bb and Bd. Therefore, on the basis of Bb and Bd geometries, we have considered new structures Bbx and Bdx including a water molecule at various position x (x = a−h). With consideration of all generations Gn, a total of 209 new structures have been identified by this additional step; among them 29 belong to the 0−15 kJ mol−1 energy range. These conformers are identified in Tables S4−S6 (see Supporting Information), in which the reason for these structures being missing following our first approach is also given. It is of interest to underline that geometry optimizations at the M06/6-311++G(d,p) level were conducted with the SCF convergence criterion equal to 10−8 and the “UltraFine” grid option for integration, as implemented in the Gaussian 09 package.62 This is necessary for systems containing many soft modes as considered here with hydrated species. Moreover, for all the clusters examined in the present work, the lowest frequencies, which correspond to hydrogen bonding network deformations and (or) relative rotation of the water molecules, are essentially lower than 30 cm−1. Obviously, the rigid rotor−

harmonic oscillator approximation implemented in the Gaussian suite of programs is probably not fully adapted to describe these low frequency modes. Thus, even with use of the precautionary “UltraFine” option, care should be taken when handling low frequencies and their closely related thermochemical parameters, particularly the vibrational contribution to entropy. Indeed, in most of the hydrated clusters, frequencies lower than 30 cm−1 account for ca. 15% of the total entropy of the system. For this reason, computed Gibbs free energy estimates should be considered only qualitatively for generations G1−G3. Therefore, throughout this work, emphasis will be essentially given to relative enthalpies. From this latter point of view, higher level calculations were performed on several key structures in order to improve the estimate of the conformer relative energies, and the corresponding hydration energies. A first, economical, approach was to perform single point energy calculation, on the M06/6-311++G(d,p) optimized geometry, at the MP2 level with the 6-311++G(2d,2p) basis set. Thermochemical corrections at zero and 298 K were taken from the M06/6-311++G(d,p) calculation. Calculated data obtained at 298 K has been used to compare with experimental values when available. This procedure has been denoted “MP2//M06” in the following lines. Basis set superposition error (BSSE) has been also considered for selected systems using the counterpoise (CP) approximation. The second, less economical, approach used the recent composite G4MP2 procedure which claims to provide accurate thermochemical parameters of a large variety of compounds or ions. The G4MP2 procedure is based on a series of single point calculations, using B3LYP/6-31G(2df,p) optimized geometries, including CCSD(T) and MP2 levels calculations and HF extrapolation based on triple and quadruple-ζ basis sets.63 NBO charges have been obtained at the M06/6-311++G(d,p) level. For the sake of clarity, and except as otherwise noted, for each generation G1−G3, only minima situated in the 0−10 kJ mol−1 range (0−15 kJ mol−1 range for G0) for ΔH°298 at the M06/6-311++G(d,p) level are described in the following. As indicated above, looser criteria (0−15 kJ mol−1 or larger ranges) were used throughout our research of minima, which guarantee the completeness of our exploration in the more restrictive range. 11529

DOI: 10.1021/acs.jpcb.5b05581 J. Phys. Chem. B 2015, 119, 11527−11539

Article

The Journal of Physical Chemistry B

OHCO van der Waals interactions are also evidenced by NCI calculations in conformers A and C, respectively. Similarly, additional stabilizations are brought by COacid···π (in D and E) and COamide···π (in B) interactions. Estimation of the Grimme’s dispersion energy within these five conformers with the B3LYP, PBE0, and M06 functionals has been performed (Table S1). The calculated empirical D3 dispersion energies within A−E are significant (between −63 and −71 kJ mol−1 at the B3LYP level). This confirms the importance of using a functional such as M06 which includes the treatment of dispersion. The stabilizing NδH···OCNH interaction occurring in conformers A, C, D, and E is particularly efficient because it involves a strong acidic hydrogen (site a) and the most basic carbonyl oxygen (site f). Accordingly, it is well-established that the proton affinity of secondary amides is ca. 100 kJ mol−1 higher than that of aliphatic acids.58 This is confirmed here by the short H···OC distances ranging from 1.61 to 1.71 Å in A, C, D, and E while it attains 1.89 Å in conformer B. In this latter structure, the second low lying conformation in this system, the amide oxygen (site f), interacts strongly with the protonated imidazolium moiety as evidenced in Figure 4. This is reflected by its low ΔH°298 (between 0.2 and 3.2 kJ mol−1, Table S3) and its low third law entropy (ca. 10 J mol−1 K−1 less than structures A, C, D, and E, Table S3). Structures A and C differ only by the orientation of the carboxylic acid moiety. This is confirmed by the torsion angles along the molecular backbone (see Table S2), which are almost equal except for the N−Cα− C−OH dihedral. A similar situation arises for conformers D and E. In both cases, the enthalpy difference between the two conformers (which consists in a subtle balance between several NamideH···OCOH, NamideH···OHCO, and COacid···π van der Waals interactions) amounts to ∼6 kJ mol−1. Comparison between the torsion angles of B and D (Table S2) indicates roughly similar conformation of the peptide backbone, whereas a rotation of about 120° around the Cβ−Cγ bond is observed in the side chain. Finally, it is noteworthy that the low lying [NAH+H]+ structures have a syn-HOCO arrangement. Accordingly, the most stable conformation bearing an anti-HOCO arrangement (a homologue of conformer B presenting an additional Hbonding of the type COOH···OC(CH3)NH) is situated 29 kJ mol−1 in H°298 above A (M06/6-311++G(d,p) calculation). In the same vein, no structure involving syn-CCNC conformation of the amide moiety has been identified in the considered AMOEBA energy range (0−60 kJ mol−1). Using the computed M06/6-311++G(d,p) G°298 values, a population ratio equal to 74/7/13/5/1% is calculated for conformers A, B, C, D, and E. A slightly different result is obtained when estimating the relative free energies by either the MP2 or the G4MP2 methods. Indeed, at this latter level the A/B/C/D/E population ratio becomes 70%/14%/7%/8%/1%, and similar results are obtained at the MP2 level. Thus, it may be concluded that protonated Nα-acetylhistidine, [NAH+H]+, is undoubtedly more stable in its conformation A but that conformers B, C, and D should be considered in the 298 K population and, obviously, as starting points during the Darwin tree approach. Calculation of the NBO charges of A−E (see Figure S1) indicates that sites a−d (e−h, respectively) have highly positive (negative, respectively) values. This confirms their ability to form hydrogen bonds with water molecules as hydrogen bond donor (acceptor, respectively) sites.

The noncovalent interactions (NCI) topological tool is based on the electron density (ρ) and its derivatives according to the formula s (ρ ) =

|∇ρ| C Fρ 4/3

where s(ρ) is the reduced electron density gradient and CF is the Fermi constant.64,65 This tool has been used to study the various types of nonbonding interactions in the most stable conformers of G0. The weak interactions within molecular systems are represented by reduced density gradient isosurfaces using s = 0.5 au. The sign of the second density Hessian matrix eigenvalue (λ2) allows separating the attractive and repulsive contributions. In order to visually differentiate between the different types of interactions, the pictures were colored in the [−0.05, 0.05] au sign(λ2)ρ range according to the following color code: the attractive weak interactions (such as hydrogen bonds) are in blue, the extremely weak interactions (such as dispersive-like van der Waals) are in green, and the repulsive interactions are in red (see Figure 4). SCF densities were considered, and the NCIPLOT software version 1.0 was used for this work.



RESULTS AND DISCUSSION Nonhydrated Nα-Acetylhistidine, [NAH+H]+, Population of Conformers G0. In this system, five conformers were identified in the 0−15 kJ mol−1 range for ΔH°298 (and ΔG°298) at the M06/6-311++G(d,p) level (structures A−E, Figures 2 and 3). B, C, D, and E are, respectively, situated 3.2, 5.4, 6.6,

Figure 3. M06/6-311++G(d,p) computed relative potential energies (ΔE°, kJ mol−1), enthalpies (ΔH°298, kJ mol−1), and Gibbs free energies (ΔG°298, kJ mol−1) of the five most stable conformers of protonated Nα-acetylhistidine [NAH+H]+ (system G0).

and 13.0 kJ mol−1 higher in enthalpies relative to A at the M06/ 6-311++G(d,p) level. Higher levels of calculations, including MP2 single point energy calculations, as well as MP2 and G4MP2 reoptimizations, show comparable relative energies, enthalpies, and free energies (see Table S3). All these structures exhibit intramolecular H-bond type interactions involving the highly polarized NδH hydrogen atom (donor site a, see Figure 1 for the nomenclature used) and a carbonyl oxygen (acceptor site e and f, Figure 1). Only this noncovalent interaction is depicted in Figure 2 for each conformer. Noncovalent interaction (NCI) analysis (Figure 4) confirms the formation of efficient intramolecular hydrogen bonds: NδH···OCNH in conformers A, C, D, and E and NδH···OCOH in conformer B. The presence of minor NamideH···OCOH and NamideH··· 11530

DOI: 10.1021/acs.jpcb.5b05581 J. Phys. Chem. B 2015, 119, 11527−11539

Article

The Journal of Physical Chemistry B

Figure 4. NCI isosurfaces indicating the nature of noncovalent interactions (blue for H-bond, green for van der Waals interactions, and red for repulsive interaction) within the structures A, B, C, D, and E (isosurfaces generated for s = 0.5 au and colored over the range −0.05 < sign(λ2)ρ < 0.05 au).

Monohydrated Nα-Acetylhistidine, [NAH+H]+(H2O), Population of Conformers G1. On the basis of the combination of hierarchical and genealogical approaches, we have found 10 structures in the 0−10 kJ mol−1 H°298 energy range at the M06/6-311++G(d,p) level (Figure 5). These structures (Ab, Ad, Bb, Cb, Ai, Ci1, Ci2, Cd, Db, and Bd) are presented in Figure 6.

conformation are observed due to formation of a H-bond with a water molecule at the b or d sites, as revealed for the torsion angles of A, Ab, and Ad which differ by less than 1° (Table S2). As expected, a larger deviation is observed for Ai, compared to A. Examination of Figure 5 and Table S4 shows that the global minimum in the G1 population is structure Ab (it represent 60−75% of the 298 K population according to MP2//M06 and G4MP2 calculations, respectively). Ad is the second most stable conformer at the M06 level. However, in agreement with the low difference in energy noted previously between A and B (Figure 3 and Table S3), the second stable structure of generation G1 at the MP2//M06 and G4MP2 levels is Bb, both in H°298 and G°298 (22% and 12% of the 298 K population at the MP2//M06 and G4MP2 levels, Table S4). Beside these two major conformers, a non-negligible amount of conformers Ad, Cb, and Db is predicted in the G1 population. Higher enthalpy conformers Ai, Ci1, Ci2, Cd, and Bd are also entropically disfavored as evidenced in Figure 5. From this point of view it may be noted that “inserted” conformers Ai, Ci1, and Ci2 are fully cyclized, constrained structures characterized by low entropies and thus by higher G°298, in particular at the G4MP2 level (Table S4). Similarly, structures Bd, Cd, and, to a lesser extent, Ad seem to be also entropically disfavored because of the double interaction existing between the water molecule and the acidic moiety (i.e., COacid··· HO(H)···HOacid) inducing a supplementary constraint. At this point, emphasis should be given to a methodological aspect of the present work. The hierarchical approach, which constitutes the first step of our method, was able to generate structures Ab, Ad, Bb, Cb, Ai, Ci1, Cd, and Bd. However, if the close proximity in energies between Ab and Bb mimics that of the couple A and B, it was of interest to see if comparable behavior pertains for inserted structures, as no Bi type structure has been obtained. Regarding the low energy structure of D relative to B and C, the absence of Db was also questioned. To examine these possibilities, the Darwinian tree has been constructed starting from G0 and considering all the H-bond acceptors and donors presented in Figure 1 as well as the

Figure 5. M06/6-311++G(d,p) computed relative potential energies (ΔE°, kJ mol−1), enthalpies (ΔH°298, kJ mol−1), and Gibbs free energies (ΔG°298, kJ mol−1) of the 10 most stable conformers of protonated Nα-acetylhistidine [NAH+H]+(H2O) (system G1).

Due to the cationic nature of the protonated imidazole moiety, formation of a H-bond between the water molecule and site b (Nε-H, Figure 1) is expected to give rise to highly stabilized structures. This is indeed observed for most of the identified low energy complexes [NAH+H]+, leading to Ab, Bb, Cb, and Db, whereas Eb is not situated in the 0−10 kJ mol−1 energy range. For a similar reason, the acidic hydrogen of the carboxylic acid group is responsible for the stability of several other structures (Ad, Bd, and Cd). Finally, a water molecule inserted in a cooperative hydrogen bond network into structures A and C, leading to conformers Ai, Ci1, and Ci2. Within a family of structures, only small variations of the 11531

DOI: 10.1021/acs.jpcb.5b05581 J. Phys. Chem. B 2015, 119, 11527−11539

Article

The Journal of Physical Chemistry B

Figure 6. Optimized geometries (M06/6-311++G(d,p)) of the 10 most stable structures of the [NAH+H]+(H2O) system (Ab, Ad, Bb, Bb, Ai, Ci1, Ci2, Cd, Db and Bd).

insertion modes i. Several new structures were indeed generated, among them the new inserted complex Ci2 and the expected Db complex (see Table S4) which are, respectively, situated only 6.8 and 8.6 kJ mol−1 above Ab at the M06/6-311++G(d,p) level. These structures were missing following our first approach either because the corresponding conformer is not a minimum in the potential energy surface of AMOEBA (case of Ci2), or due to an erroneous value of the relative energy at the force field level (case of Db). Another noticeable result is that inserted structure Bi is a high energy species (ΔH°298 = 12.9 kJ mol−1 at the M06/6-311++G(d,p) level, data not shown). This demonstrates the low efficiency of the internal hydrogen bond (H)OH···OC(OH) to stabilize such structures. Dihydrated Nα-Acetylhistidine, [NAH+H]+(H2O)2, Population of Conformers G2. There were 12 minima in the ΔH°298 = 0−10 kJ mol−1 range identified at the M06/6-311++G(d,p) level (Figures 7 and 8, Table S5). Six are related to conformer A, 2 to conformer B, and 4 to conformer C. In terms of H°, five low lying structures are predicted to be in the ΔH°298 < 5 kJ mol−1 range: Abd, Bbf, Abi, Adi, and Cbi1. Three of these structures, namely, Abd, Abi, and Cbi1, are also the lowest in the G°298 scale and are predicted to represent 90− 70% (MP2//M06 or G4MP2 calculations, Table S5) of the 298 K population. In other words, structures Bbf and Adi are entropically disfavored, spectacularly in the case of complex Bbf. A comparable entropy effect should be mentioned in the case of structure Ade. The situation encountered for Bbf and Ade may be easily understood by considering that they are fully cyclized structures in which the two water molecules act cooperatively to establish a bridge between two sites. In this constrained structure, no rotation of the water molecules is possible, and no frequencies are indeed present below 30 cm−1

Figure 7. M06/6-311++G(d,p) computed relative potential energies (ΔE°, kJ mol−1), enthalpies (ΔH°298, kJ mol−1), and Gibbs free energies (ΔG°298, kJ mol−1) of the 12 most stable conformers of the protonated Nα-acetylhistidine [NAH+H]+(H2O)2 (system G2).

in the M06/6-311++G(d,p) calculations. It may be noted that the possibility of a two water molecule bridge is not offered between sites b and d/e for A- or C-type containing conformers of G2 because these two sites are too distant from each other, as illustrated by the 9.73 Å distance, which separates the two water molecules in Abd. The observation that Abd is the most stable conformer in G2 generation is in line with the fact that A and Ab are the global minima in the G0 and G1 generations. It is noteworthy that the relative energy of Cbd (∼6 kJ mol−1) is similar to the relative energy of C in G0 and Cb in G1 indicating that b and d Hbonds hold similar strengths in G1 and G2 generations. By contrast, the relative energies of Bbd (9.6 kJ mol−1) are significantly larger than that observed for B or for Bb (3−4 kJ mol−1). This reveals that the strength of the d type H-bond is 11532

DOI: 10.1021/acs.jpcb.5b05581 J. Phys. Chem. B 2015, 119, 11527−11539

Article

The Journal of Physical Chemistry B

Figure 8. Optimized M06/6-311++G(d,p) geometries of conformers of [NAH+H]+(H2O)2 in the 0−10 kJ mol−1 enthalpy range (Abd, Bbf, Abi, Adi, Cbi1, Cbd, Ci1d, Ci2d, Bbd, Abw1, Abw2, Ade).

significantly reduced in Bbd with respect to Abd, presumably because the acidic group is also involved in a H-bond with the imidazolium NδH group in Bbd (Figure 8). Mixed structures Abi, Adi, Cbi1, and Ci1d may be seen as cyclized forms Ai and Ci1 stabilized by an additional water molecule involved in either b or d H-bond interactions. Comparison of the enthalpies reported in Tables S4 and S5 reveals that Ai is slightly more efficiently stabilized by these hydrogen bonds than Ci1 and that the b H-bond is stronger than d as already noted. Last, the second water molecule can be located in the second solvation shell through formation of an H-bond with the first water molecule. Different relative orientation of these water molecules induces several conformers noted w1 and w2. From this point of view, the second shell isomers Abw1 and Abw2 are the two other structures that may be of interest in describing the G2 generation. These two structures have high relative enthalpies (9.8 and 9.9 kJ mol−1 at the M06/6-311++G(d,p) level, Figure 7) but are entropically favored, leading to rather low relative free energies, in particular at the G4MP2 level (2.1 and 1.6 kJ mol−1, Table S5). Trihydrated Nα-Acetylhistidine, [NAH+H]+(H2O)3, Population of Conformers G3. The combination of genealogical and hierarchical methods leads to 12 minima in the 0−10 kJ mol−1 enthalpy range at the M06/6-311++G(d,p) level (Figures 9 and 10). Due to their low relative free energies, 5 other conformers situated in the 10−15 kJ mol−1 enthalpy range are included in Figure 9. The G3 lowest lying conformer Abdi results from a water molecule insertion into the intramolecular NδH···OCamide H-bond of the lowest G2 conformer Abd. The following G3 low lying structures may be seen as resulting from hydration at site d of Bbf and Cbi1, or at site b of Ci2d, conformers which were

Figure 9. Computed relative potential energies (ΔE°, kJ mol−1), enthalpies (ΔH°298, kJ mol−1), and Gibbs free energies (ΔG°298, kJ mol−1) of the most stable conformers of protonated Nα-acetylhistidine [NAH+H]+(H2O)3 (system G3).

already identified among the most stable in the G2 system. Despite the absence of low lying E-type conformers in G1 and G2, the conformer Ebid, which is a third generation child of the G0 low lying structure E, is located only 5.5 kJ mol−1 higher in enthalpy than Abdi. Structures involving water molecules in the second solvation shell (i.e., Abdw1, Abdw2, Abdw3, Abiw2, Abiw1, Cbi1w1, Cbi1w2, Cbdw1, and Cbdw2) are situated more than 7 kJ mol−1 higher in enthalpy compared to Abdi, and the first stable structure having a water molecule in the third solvation shell (namely Abw1w3, data not shown) has been found to be more than 20 kJ mol−1 above Abdi. Furthermore, we observed that two water molecules can establish a bridge between sites d and e on structure A (Abde1 and Abde2) or between sites b and f on structure B (Bbfd). This situation is 11533

DOI: 10.1021/acs.jpcb.5b05581 J. Phys. Chem. B 2015, 119, 11527−11539

Article

The Journal of Physical Chemistry B

Figure 10. Optimized M06/6-311++G(d,p) geometries of conformers of [NAH+H]+(H2O)3 in the 0−10 kJ mol−1 enthalpy range (Abdi, Bbfd, Cbi1d, Abde1, Abde2, Ebid, Cbdi2, Abdw1, Abdw2, Abdw3, Fbdi, Abiw2).

Figure 11. Comparative view of the Darwinian trees of the major hydrated [NAH+H]+ structures (the most stable species in term of enthalpy are indicated in bold while less stables are noted into parentheses).

cyclized structures Bbfd, Abde1, and Abde2 are disfavored by their low entropies. General Overview. From this study of the four generations G0−G3, we note that the successive addition of water molecules to the protonated Nα-acetylhistidine mainly preserved the backbone of the lowest energy structures. Conformer of type A remains the most stable structure of all generations following the enthalpy and Gibbs free energy scales. It should however be noted that, in the MP2//M06 enthalpy scale, conformer Bbf is slightly more stable than Abd in the G2 generation, and Bbfd is as stable as Abdi in the G3 generation. Last, conformers of type C become competitive with A and B only when two molecules of water are attached (see Cbi1 in Figure 8 and Table S5). A summary of the most significant structures involved in the

comparable to that encountered with Ade and Bbf in the previously examined G2 system. A new conformation of protonated Nα-acetylhistidine, which does not include intramolecular H-bonds, is identified as Fbid (Figure 10). Complementary calculations confirm that structure F (G0), Fb and Fi (G1), and Fbi (G2) are high energy (>15 kJ mol−1) conformers in their respective generations. Considering the ΔH°298 values, two structures are clearly identified as the more stable: Abdi and Bbfd. On the basis of the MP2 Gibbs free energy calculations (Table S6), approximately 80% of the G3 Boltzmann population at 298 K is predicted to be constituted of structures Abdi, Cbi1d, Abdw1, Abiw1, Cbi1w2, and, to a lesser extent, Cbdi2 and Abiw2 since 11534

DOI: 10.1021/acs.jpcb.5b05581 J. Phys. Chem. B 2015, 119, 11527−11539

Article

The Journal of Physical Chemistry B

Table 1. Computed and Experimental Hydration Enthalpies (ΔhH°298, kJ mol−1), Gibbs Free Energies (ΔhG°298, kJ mol−1), and Entropies (ΔhS°298, J K−1 mol−1) of Model Cations [M + H]+(H2O)n−1 ([M + H]+ = CH3NH3+, Imidazolium) [M + H]+(H2O)n−1 CH3NH3+

CH3NH3+(H2O)

CH3NH3+(H2O)2

imidazolium cation

a

method

ΔhH°298

ΔhG°298

ΔhS°298

M06/6-311++G(d,p) MP2//M06 G4MP2 exptl exptl, av M06/6-311++G(d,p) MP2//M06 G4MP2 exptl exptl, av M06/6-311++G(d,p) MP2//M06 G4MP2 exptl exptl, av M06/6-311++G(d,p) MP2//M06 G4MP2 exptl

79.4 71.2 68.6 62,a 70,b 79c 70 ± 8 66.1 59.8 57.9 55,a 61,b 61c 59 ± 4 56.7 51.7 50.1 51,a 51,b 52c 51 ± 3 68.6 65.4 60.7 62 ± 4e

51.1 42.9 41.0 38,a 43,b 46,c 45d 43 ± 4 35.8 29.4 27.3 28,a 31,b 28,c 31d 30 ± 3 26.1 21.2 22.6 21,a 21,b 19,c 22d 21 ± 3 36.0 32.8 27.6 32 ± 4e

95 95 93 83,a 91,b 110c 95 ± 14 102 102 103 90,a 102,b 112c 101 ± 11 103 103 92 103,a 101,b 111c 105 ± 8 109 109 111 100 ± 8e

Reference 4. bReference 66. cReference 67. dReference 68. eReference 69.

imidazolium cation, the gap between hydration enthalpies or Gibbs free energies calculated using the composite G4MP2 method and the M06/6-311++G(d,p) functional is still observed, whereas MP2//M06 results are in between. The G4MP2 method provides the ΔhH° value in best agreement with experiment; however, MP2//M06 also gives a result falling in the experimental window, though in the upper side. Last, only MP2//M06 provides ΔhG° value in agreement with experiment for imidazolium. Concerning hydration entropies, Table 1 indicates that both B3LYP/6-31+G(d,p) (for G4MP2 calculations) and M06/6-311++G(d,p) functionals give comparable results since the differences in ΔhS° never exceed 2 J mol−1 K−1, except for CH3NH3+(H2O)2 where it attains 10 J mol−1 K−1. There are however noticeable differences between calculated and experimental entropy values; it attains for example ca. 10 J mol−1 K−1 for the imidazolium cation. However, it may be observed that these differences are close to the reported experimental errors and should include uncertainties on the calculated entropies values as underlined in the computational section. Nevertheless, due to these uncertainties, only enthalpies will be discussed hereafter for the hydration process of [NAH+H]+ ion. Hydration enthalpies of variously hydrated [NAH+H]+ ions are presented in Table S7. As observed with the reference systems (Table 1), a systematic shift in ΔhH° is associated with the three levels of calculation used here: the G4MP2 method leading to the lowest values whereas M06/6-311++G(d,p) results are situated 9.2 ± 0.9 kJ mol−1 higher, on the average. This latter overestimate is partly compensated by the simple point MP2 calculation as shown in the third column of Table S7. Note that BSSE estimates done on selected structures reveal ΔhH° calculated at the MP2//M06 level closer to the G4MP2 results (Table S7). As expected, the ΔhH° values calculated for Ab, Bb, Cb, and Db are almost similar indicating that the strength of the Hbond through imidazolium site b (NεH group) is independent of the conformation of the side chain. The average ΔhH° values are 57.7 ± 0.9 kJ mol−1 (M06/6-311++G(d,p) level) or 50.5 ±

population of generations G0 to G3 is presented in Figure 11. In this representation, vertical arrows correspond to addition of a water molecule at site d and diagonal arrows to either bonding to the sites b or f, or insertion of a water molecule into a preexisting intramolecular hydrogen bond. Thermochemistry. Enthalpies (Δ h H° T ), entropies (ΔhS°T), and Gibbs free energies (ΔhG°T) of hydration of species [M + H]+(H2O)n−1, in the gas phase, at temperature T, are defined by considering reactions 1: [M + H]+ (H 2O)n → [M + H]+ (H 2O)n − 1 + H 2O (n = 1 − 3)

(1)

These quantities will be estimated by considering 1 mol of reactant and products in the gas phase at the temperature of 298 K. In order to estimate the strength of a given noncovalent interaction, for a given conformer [M + H]+(H2O)n, each mole of reactant and product is considered to contain only one conformation of generation Gn and the related conformation in the parent generation Gn−1. In fact, realistic hydration thermochemical parameters should consider 1 mol of reactant and products in thermal equilibrium at 298 K. This may be calculated, assuming a Boltzmann distribution of conformers, by estimating the molar fractions xi of each conformer from the corresponding Gibbs free energies. However, the difference in the corresponding ΔhG°T (less than 2 kJ mol−1) is close to (or lower than) the possible computational errors due to the uncertainties on the entropy calculations. In order to test the validity of the levels of theory used in the present study, comparison with experimental thermochemical parameters obtained for hydration of protonated methylamine and protonated imidazole is done in Table 1. Clearly, in the case of protonated methylamine, MP2//M06 and G4MP2 provide comparable hydration enthalpies ΔhH° and Gibbs free energies ΔhG°; moreover, these values are very close to the experimental data (MAD = 1.2 kJ mol−1). In contrast, the M06/6-311++G(d,p) method appears to slightly overestimate ΔhH° and ΔhG° (MAD = 7.6 kJ mol−1). When looking at the 11535

DOI: 10.1021/acs.jpcb.5b05581 J. Phys. Chem. B 2015, 119, 11527−11539

Article

The Journal of Physical Chemistry B 0.8 kJ mol−1 (G4MP2 level). It is interesting to note that these values are 10 kJ mol−1 lower than that calculated for hydration of the imidazolium cation (Table 1). Thus, the intramolecular H-bond involving site a (NδH group) in species containing the A, B, C, and D backbone decreases the strength of imidazolium···OH2 H-bond involving site b. To confirm this hypothesis, we calculated the ΔhH° of imidazolium cation complexed with methylacetamide via one of its N−H bonds (complex 1, Figure 12). The monohydration enthalpy

molecule and site d in A and C should certainly be attributed to the presence of the positively charged imidazolium moiety. This conclusion may be extended to Abd and Cbd species (G2 generation) and to Cbi1d specie (G3 generation) since hydration at site d is associated with ΔhH° of ∼55 and ∼52 kJ mol−1 at the M06/6-311++G(d,p) level. Comparison with previous studies suggests that a similar effect takes place with ammonium group in [Gly+H]+, [Ala+H]+, and [Pro+H]+. Indeed, these protonated amino acids show hydration enthalpies for the first water molecule bound to the COOH group at about 46−50 kJ mol−1 at the MP2 level.33−35 On the opposite, hydration at the site d for B-type conformers (52.2 kJ mol−1 for Bd and 49.9 kJ mol−1 for Bbd at the M06/6-311+ +G(d,p) level, Table S7) is reduced by ∼5 kJ mol−1 compared to A- and C-type conformers, probably due to the intramolecular COcarboxylic···NHimidazolium H-bond which slightly reduces the acidity of the OHcarboxylic group. Insertion of a water molecule in the intramolecular H-bonds in A, B, and C (i-type sites) corresponds to hydration enthalpies of 52, 49, and 57 kJ mol−1 at the M06/6-311+ +G(d,p) level (reduced to 43, 40, and 49 kJ mol−1 at the G4MP2 level). These values are also comparable to the ΔhH° of sites b and d described above. As a consequence, this close similitude in ΔhH° values renders problematic the use of experimental hydration enthalpies to locate the added water molecule. Finally, the building of the second hydration shell, in systems including w-type interactions, is associated with a fairly constant ΔhH° value of ∼45 kJ mol−1 (M06/6-311++G(d,p) level) or 36 kJ mol−1 (G4MP2 calculation). The subsequent interaction of a water molecule is thus progressively weaker as already mentioned in the literature for systems containing ammonium groups.4,7−10,70,71

Figure 12. Optimized geometries (M06/6-311++G(d,p) calculations) of model systems (hydration enthalpies are in kJ mol−1).

calculated at the M06/6-311++G(d,p) level for 1 is equal to 57.6 kJ mol−1, i.e., 11 kJ mol−1 below the value calculated for imidazolium at the same level of theory (Table 1). This leads to the conclusions that the amide···imidazolium H-bond reduces the strength of the imidazolium···OH2 H-bond by ∼10 kJ mol−1 and that the Cα-COOH group, not present in 1, does not influence the coordination of the water molecule. Comparison with previous studies indicates that the hydration enthalpies calculated for Ab, Bb, Cb, and Db are in the same range as the one computed for the third water molecule bound to the protonated ammonium group in [Gly+H]+ and [Ala+H]+ (about 50 kJ mol−1 at the MP2 level).33,34 A similar value is also obtained for the second water molecule bound to the secondary ammonium group of [Pro+H]+ (46 kJ mol−1 at the MP2 level).35 The hydration enthalpies mentioned in Table S7 are computed for the water molecule bound to the site corresponding to the last letter of the name (i.e., for example site d for Ad or Abd). From our data, it is nevertheless possible to compute hydration enthalpies for other sites. In particular, for the G2 generation, this can be achieved for site b. Dissociation of the water molecule bound to site b from Abd, Abi, Cbi1, Cbd, and Bbd leads, respectively, to Ad, Ai, Ci1, Cd, and Bd. Following this procedure (data not shown in Table S7), hydration of site b for G2 is associated with ΔhH° very close to that discussed above for the G1 generation. Accordingly, the average ΔhH° values are 56.7 ± 0.9 kJ mol−1 (M06/6-311++G(d,p) level) or 49.5 ± 0.8 kJ mol−1 (G4MP2 level). Surprisingly, examination of Table S7 reveals that the ΔhH° values of conformers Ad (56.8 kJ mol−1 at the M06/6-311+ +G(d,p) level) and Cd (56.1 kJ mol−1 at the M06/6-311+ +G(d,p) level) are similar to those computed for the site b, even though, in these species, the H-bond implies the neutral COOH group. To understand this result we calculated the hydration enthalpy of acetic acid itself (complex 2, Figure 12). The results are ΔhH°(2) = 42.4, 35.7, and 30.6 kJ mol−1 at the M06/6-311++G(d,p), MP2//M06, and G4MP2 levels, respectively. Thus, ΔhH°(2) is 14 kJ mol−1 lower than ΔhH°(Ad) or ΔhH°(Cd). A comparable difference is observed when considering the more realistic Nα-acetyl-acetamide structure 3. The noticeable strengthening of the H-bond between a water



CONCLUSION Using the hierarchical and the genealogical methods, we have exhaustively explored the [NAH+H]+ (H2O)n (n = 0−3) potential energy surfaces. Among the 38 structures identified for the four generations in the ΔH°298 = 0−10 kJ mol−1 range (at the (M06/6-311++G(d,p) level), 27 (i.e., 71%) have been obtained from the hierarchical approach, and the remaining 11 conformers (29%) result from the Darwinian tree method. This observation clearly demonstrates the usefulness of the double hierarchical/genealogical approach developed in the present study. It also confirms the ability of the AMOEBA force field to correctly describe systems dominated by noncovalent interactions. Three [NAH+H]+ structures, namely, A, B, and C, are the major backbone of all the most stable hydrated complexes identified. The first water molecule is preferentially located on the imidazolium NεH group (site b). Two sites are in competition for the second and third hydration: the hydroxyl of the COOH group (site d) and the insertion into the intramolecular H-bonds which already involve the cationic group and one carbonyl group of the side chain. Compared to model systems, the intramolecular H-bond and the long-range influence of the positive charge explained the computed hydration enthalpies values. Indeed, the former effect has an unfavorable effect on hydration site b while the latter improves the accepting ability of site d. Therefore, hydration enthalpies of A, B, and C on site b are 10 kJ mol−1 lower than that associated with isolated imidazolium cation whereas water bonding on site d leads to hydration enthalpies that are 10 kJ 11536

DOI: 10.1021/acs.jpcb.5b05581 J. Phys. Chem. B 2015, 119, 11527−11539

Article

The Journal of Physical Chemistry B mol−1 higher than that associated with acetic acid. Thereby, these two sites become isoenergetic. This highlights that experimental data do not always allow researchers to clearly distinguish the preferential sites for water coordination with the hydration enthalpy criteria alone. Experimental and G4MP2 computed ΔhH° of model cations are in good agreement. Systematic shifts close to 9 kJ mol−1 are observed with the M06/6-311++G(d,p) level. However, MP2//M06 or G4MP2 calculations allow tightening this gap. The entropic contribution plays a noticeable role in distinguishing between cyclized and pseudolinear structures, with the latter being favored in Gibbs free energy but also subject to larger uncertainties due to the difficulties to accurately calculate the low vibrational frequencies.



(8) Wincel, H. Hydration Energies of Protonated Amino Acids. Chem. Phys. Lett. 2007, 439, 157−161. (9) Wyttenbach, T.; Liu, D.; Bowers, M. T. Hydration of Small Peptides. Int. J. Mass Spectrom. 2005, 240, 221−232. (10) Wyttenbach, T.; Bowers, M. T. Hydration of Biomolecules. Chem. Phys. Lett. 2009, 480, 1−16. (11) Liu, D.; Wyttenbach, T.; Barran, P. E.; Bowers, M. T. Sequencial Hydration of Small Protonated Peptides. J. Am. Chem. Soc. 2003, 125, 8458−8464. (12) Gao, B.; Wyttenbach, T.; Bowers, M. T. Hydration of Protonated Aromatic Amino Acids: Phenylalanine, Tryptophan, and Tyrosine. J. Am. Chem. Soc. 2009, 131, 4695−4701. (13) Kohtani, M.; Jarrold, M. F. The Initial Steps in the Hydration of Unsolvated Peptides: Water Molecule Adsorption on Alanine-Based Helices and Globules. J. Am. Chem. Soc. 2002, 124, 11148−11158. (14) Lemoff, A. S.; Bush, M. F.; Wu, C.; Williams, E. R. Binding Energies of Water to Sodiated Valine and Structural Isomers in the Gas Phase: The Effect of Proton Affinity on Zwitterion Stability. J. Am. Chem. Soc. 2003, 125, 13576−13584. (15) Ye, S. J.; Moision, R. M.; Armentrout, P. B. Sequential Bond Energies of Water to Sodium Proline Cation. Int. J. Mass Spectrom. 2006, 253, 288−304. (16) Prell, J. S.; Chang, T. M.; O’Brien, J. T.; Williams, E. R. Hydration Isomers of Protonated Phenylalanine and Derivatives: Relative Stabilities from Infrared Photodissociation. J. Am. Chem. Soc. 2010, 132, 7811−7819. (17) Scuderi, D.; Bakker, J. M.; Durand, S.; Maitre, P.; Sharma, A.; Martens, J. K.; Nicol, E.; Clavaguera, C.; Ohanessian, G. Structure of Singly Hydrated, Protonated Phospho-Tyrosine. Int. J. Mass Spectrom. 2011, 308, 338−347. (18) Atkins, C. G.; Banu, L.; Rowsell, M.; Blagojevic, V.; Bohme, D. K.; Fridgen, T. D. Structure of [Pb(Gly-H)]+ and the Monosolvated Water and Methanol Solvated Species by Infrared Multiple-Photon Dissociation Spectroscopy, Energy-Resolved Collision-Induced Dissociation, and Electronic Structure Calculations. J. Phys. Chem. B 2009, 113, 14457−14464. (19) Yao, Y.; Chen, D.; Zhang, S.; Li, Y.; Tu, P.; Liu, B.; Dong, M. Building the First Hydration Shell of Deprotonated Glycine by the MCMM and ab Initio Methods. J. Phys. Chem. B 2011, 115, 6213− 6221. (20) Wei, Z.; Chen, D.; Zhao, H.; Li, Y.; Zhu, J.; Liu, B. Ab Initio Investigation of the First Hydration Shell of Protonated Glycine. J. Chem. Phys. 2014, 140, 085103. (21) Aikens, C. M.; Gordon, M. S. Incremental Solvation of Nonionized and Zwitterionic Glycine. J. Am. Chem. Soc. 2006, 128, 12835−12850. (22) Calvo, F.; Douady, J. Stepwise Hydration and Evaporation of Adenosine Monophosphate Nucleotide Anions: A Multiscale Theoretical Study. Phys. Chem. Chem. Phys. 2010, 12, 3404−3414. (23) Lemoff, A. S.; Bush, M. F.; Wu, C. C.; Williams, E. R. Structures and Hydration Enthalpies of Cationized Glutamine and Structural Analogues in the Gas Phase. J. Phys. Chem. A 2005, 109, 1903−1910. (24) Kamariotis, A.; Boyarkin, O. V.; Mercier, S. R.; Beck, R. D.; Bush, M. F.; Williams, E. R.; Rizzo, T. R. Infrared Spectroscopy of Hydrated Amino Acids in the Gas Phase: Protonated and Lithiated Valine. J. Am. Chem. Soc. 2006, 128, 905−916. (25) Kaminsky, J.; Jensen, F. Force Field Modeling of Amino Acid Conformational Energies. J. Chem. Theory Comput. 2007, 3, 1774− 1788. (26) Ren, P.; Ponder, J. W. Consistent Treatment of Inter- and Intramolecular Polarization in Molecular Mechanics Calculations. J. Comput. Chem. 2002, 23, 1497−1506. (27) Ponder, J. W.; Wu, C. J.; Ren, P.; Pande, V. S.; Chodera, J. D.; Schnieders, M. J.; Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio, R. A.; et al. Current Status of the AMOEBA Polarizable Force Field. J. Phys. Chem. B 2010, 114, 2549−2564. (28) Riffet, V.; Frison, G.; Bouchoux, G. Acid-base Thermochemistry of Gaseous Oxygen and Sulfur Substituted Amino Acids (Ser, Thr, Cys, Met). Phys. Chem. Chem. Phys. 2011, 13, 18561−18580.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b05581. Complete author lists for refs 27 and 62; calculated NBO charges, D3 dispersion correction, torsion angles, relative potential energies, enthalpies, Gibbs free energies, and hydration enthalpies; [NAH+H]+(H2O)n Cartesian coordinates calculated at the M06/6-311++G(d,p) level and computed absolute energies, enthalpies, and Gibbs free energies (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +33 (0)1 69 33 48 34. Fax: +33 (0)1 69 33 48 03. Present Address

† (V.R.) Laboratoire de Chimie Théorique, Université Pierre et Marie Curie and CNRS, F-75005 Paris, France.

Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was performed using HPC resources from GENCICINES (Grant x2013087105). Dr. E. Rezabal is acknowledged for her help to NCI calculations.



REFERENCES

(1) Pratt, L. R.; Pohorille, A. Hydrophobic Effects and Modeling of Biophysical Aqueous Solution Interfaces. Chem. Rev. 2002, 102, 2671− 2692. (2) Levy, Y.; Onuchic, J. N. Water Mediation in Protein Folding and Molecular Recognition. Annu. Rev. Biophys. Biomol. Struct. 2006, 35, 389−415. (3) Blades, A. T.; Klassen, J. S.; Kebarle, P. Free Energies of Hydration in the Gas Phase of the Anions of Some Oxo Acids of C, N, S, P, Cl, and I. J. Am. Chem. Soc. 1995, 117, 10563−10571. (4) Gilligan, J. J.; Lampe, F. W.; Nguyen, V. Q.; Vieira, N. E.; Yergey, A. L. Hydration of Alkylammonium Ions in the Gas Phase. J. Phys. Chem. A 2003, 107, 3687−3691. (5) Meot-Ner, M. The Ionic Hydrogen Bond. Chem. Rev. 2005, 105, 213−284. (6) Meot-Ner, M. Update 1 of: Strong Ionic Hydrogen Bonds. Chem. Rev. 2012, 112, PR22−PR103. (7) Wincel, H. Hydration of Gas-Phase Protonated Alkylamines, Amino Acids and Dipeptides Produced by Electrospray. Int. J. Mass Spectrom. 2006, 251, 23−31. 11537

DOI: 10.1021/acs.jpcb.5b05581 J. Phys. Chem. B 2015, 119, 11527−11539

Article

The Journal of Physical Chemistry B (29) Joshi, K.; Semrouni, D.; Ohanessian, G.; Clavaguéra, C. Structures and IR Spectra of the Gramicidin S Peptide: Pushing the Quest for Low Energy Conformation. J. Phys. Chem. B 2012, 116, 483−490. (30) Gregori, B.; Guidoni, L.; Chiavarino, B.; Scuderi, D.; Nicol, E.; Frison, G.; Fornarini, S.; Crestoni, M. E. Vibrational Signatures of SNitrosoglutathione as Gaseous, Protonated Species. J. Phys. Chem. B 2014, 118, 12371−12382. (31) Semrouni, D.; Ohanessian, G.; Clavaguéra, C. Structural, Energetics and Dynamical Properties of Sodiated Oligoglycines: Relevance of a Polarizable Force Field. Phys. Chem. Chem. Phys. 2010, 12, 3450−3462. (32) Michaux, C.; Wouters, J.; Jacquemin, D.; Perpète, E. A. A Theoretical Investigation of the Hydrated Glycine Cation Energetics and Structures. Chem. Phys. Lett. 2007, 445, 57−61. (33) Michaux, C.; Wouters, J.; Perpète, E. A.; Jacquemin, D. Microhydration of Protonated Glycine: An ab Initio Family Tree. J. Phys. Chem. B 2008, 112, 2430−2438. (34) Michaux, C.; Wouters, J.; Perpète, E. A.; Jacquemin, D. Modeling the Microhydration of Protonated Alanine. J. Phys. Chem. B 2008, 112, 9896−9902. (35) Michaux, C.; Wouters, J.; Perpète, E. A.; Jacquemin, D. Stepwise Hydration of Protonated Proline. J. Phys. Chem. B 2008, 112, 7702− 7705. (36) Michaux, C.; Wouters, J.; Perpète, E. A.; Jacquemin, D. Ab Initio Investigation of the Hydration of Deprotonated Amino Acids. J. Am. Soc. Mass Spectrom. 2009, 20, 632−638. (37) Campo-Cacharrόn, A.; Cabaleiro-Lago, E. M.; Rodríguez-Otero, J. Effects of Microhydration on the Characteristics of Cation-Phenol Complexes. Theor. Chem. Acc. 2012, 131, 1290−1302. (38) Demireva, M.; O’Brien, J. T.; Williams, E. R. Water-Induced Folding of 1,7-Diammoniumheptane. J. Am. Chem. Soc. 2012, 134, 11216−11224. (39) Rodziewicz, P.; Doltsinis, N. L. Ab Initio Molecular Dynamics Free-Energy Study of Microhydration Effects on the Neutral− Zwitterion Equilibrium of Phenylalanine. ChemPhysChem 2007, 8, 1959−1968. (40) Chang, T. M.; Prell, J. S.; Warrick, E. R.; Williams, E. R. Where’s the Charge? Protonation Sites in Gaseous Ions Change with Hydration. J. Am. Chem. Soc. 2012, 134, 15805−15813. (41) Bush, M. F.; Prell, J. S.; Saykally, R. J.; Williams, E. R. One Water Molecule Stabilizes the Cationized Arginine Zwitterion. J. Am. Chem. Soc. 2007, 129, 13544−13553. (42) Jockusch, R. A.; Lemoff, A. S.; Williams, E. R. Hydration of Valine-Cation Complexes in the Gas Phase: On the Number of Water Molecules Necessary to Form a Zwitterion. J. Phys. Chem. A 2001, 105, 10929−10942. (43) Bachrach, S. M.; Nguyen, T. T.; Demoin, D. W. Microsolvation of Cysteine: A Density Functional Theory Study. J. Phys. Chem. A 2009, 113, 6172−6181. (44) Kim, J. Y.; Ahn, D. S.; Park, S. W.; Lee, S. Gas Phase Hydration of Amino Acids and Dipeptides: Effects on the Relative Stability of Zwitterion vs. Canonical Conformers. RSC Adv. 2014, 4, 16352− 13361. (45) Nagornova, N. S.; Rizzo, T. R.; Boyarkin, O. V. Interplay of Intra- and Intermolecular H-Bonding in a Progressively Solvated Macrocyclic Peptide. Science 2012, 336, 320−323. (46) Jacquemin, D.; Michaux, C.; Perpète, E. A.; Frison, G. Comparison of Microhydration Methods: Protonated Glycine as a Working Example. J. Phys. Chem. B 2011, 115, 3604−3613. (47) Andersson, Y.; Langreth, D. C.; Lundqvist, B. I. Van der Waals Interactions in Density-Functional Theory. Phys. Rev. Lett. 1996, 76, 102−105. (48) Langreth, D. C.; Dion, M.; Rydberg, H.; Schröder, E.; Hyldgaard, P.; Lundqvist, B. I. Van der Waals Density Functional Theory with Applications. Int. J. Quantum Chem. 2005, 101, 599−610. (49) Sato, T.; Tsuneda, T.; Hirao, K. Van der Waals Interactions Studied by Density Functional Theory. Mol. Phys. 2005, 103, 1151− 1164.

(50) Sun, Y. Y.; Kim, Y.-H.; Lee, K.; Zhang, S. B. Accurate and Efficient Calculation of van der Waals Interactions within Density Functional Theory by Local Atomic Potential Approach. J. Chem. Phys. 2008, 129, 154102. (51) Grimme, S. Accurate Description of van der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463−1473. (52) 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 H-Pu. J. Chem. Phys. 2010, 132, 154104. (53) Jurečka, P.; Č erný, J.; Hobza, P.; Salahub, D. R. Density Functional Theory Augmented with an Empirical Dispersion Term. Interaction Energies and Geometries of 80 Noncovalent Complexes Compared with ab Initio Quantum Mechanics Calculations. J. Comput. Chem. 2007, 28, 555−569. (54) Steinmann, S. N.; Corminboeuf, C. Comprehensive Benchmarking of a Density-Dependent Dispersion Correction. J. Chem. Theory Comput. 2011, 7, 3567−3577. (55) Zhao, Y.; Truhlar, D. G. Density Functionals with Broad Applicability in Chemistry. Acc. Chem. Res. 2008, 41, 157−167. (56) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (57) Addario, V.; Guo, Y.; Chu, I. V.; Ling, Y.; Ruggerio, G.; Rodriquez, C. F.; Hopkinson, A. C.; Siu, K. W. M. Proton Affinities of Methyl Esters of N-Acetylated Amino Acids. Int. J. Mass Spectrom. 2002, 219, 101−114. (58) Bouchoux, G. Gas-phase Basicities of Polyfunctional Molecules. Part 4: Carbonyl Groups as Basic Sites. Mass Spectrom. Rev. 2015, 34, 493−534. (59) Bouchoux, G. Gas-phase Basicities of Polyfunctional Molecules. Part 3: Amino Acids. Mass Spectrom. Rev. 2012, 31, 391−435. (60) Riffet, V.; Bouchoux, G. Gas-phase Structures and Thermochemistry of Neutral Histidine and its Conjugated Acid and Base. Phys. Chem. Chem. Phys. 2013, 15, 6097−6106. (61) Tinker, version 5.0; Software Tools for Molecular Design; Jay Ponder Lab, Department of Chemistry, Washington University: Saint Louis, MO, 2009. (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 B.01; Gaussian, Inc.: Wallingford, CT, 2010. (63) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gaussian-4 Theory using Reduced Order Perturbation Theory. J. Chem. Phys. 2007, 127, 124105. (64) Johnson, E. R.; Keinan, S.; Mori-Sánchez, P.; Contreras-García, J.; Cohen, A. J.; Yang, W. Revealing Noncovalent Interactions. J. Am. Chem. Soc. 2010, 132, 6498−6506. (65) Contreras-García, J.; Johnson, E. R.; Keinan, S.; Chaudret, R.; Piquemal, J. P.; Beratan, D. N.; Yang, W. NCIPLOT: A Program for Plotting Noncovalent Interaction Regions. J. Chem. Theory Comput. 2011, 7, 625−632. (66) Meot-Ner, M. The Ionic Hydrogen Bond and Ion Solvation. 2. Solvation of Onium Ions by One to Seven Water Molecules. Relations between Monomolecular, Specific, and Bulk Hydrogen. J. Am. Chem. Soc. 1984, 106, 1265−1272. (67) Lau, Y. K.; Kebarle, P. Hydrogen Bonding Solvent Effect on the Basicity of Primary Amines CH3NH2, C2H5NH2, and CF3CH2NH2. Can. J. Chem. 1981, 59, 151−155. (68) Banic, C. M.; Iribarne, J. V. Equilibrium Constants for Clustering of Neutral Molecules about Gaseous Ions. J. Chem. Phys. 1985, 83, 6432−6448. (69) Meot-Ner, M. Models for Strong Interactions in Proteins and Enzymes. 2. Interactions of Ions with the Peptide Link and with Imidazole. J. Am. Chem. Soc. 1988, 110, 3075−3080. 11538

DOI: 10.1021/acs.jpcb.5b05581 J. Phys. Chem. B 2015, 119, 11527−11539

Article

The Journal of Physical Chemistry B (70) Klassen, J. S.; Blades, A. T.; Kebarle, P. Determinations of Ion− Molecule Equilibria Involving Ions Produced by Electrospray. Hydration of Protonated Amines, Diamines, and some Small Peptides. J. Phys. Chem. 1995, 99, 15509−15517. (71) Liu, D.; Wyttenbach, T.; Bowers, M. Hydration of Protonated Primary Amines: Effects of Intermolecular and Intramolecular Hydrogen Bonds. T. Int. J. Mass Spectrom. 2004, 236, 81−90.

11539

DOI: 10.1021/acs.jpcb.5b05581 J. Phys. Chem. B 2015, 119, 11527−11539