Ab initio investigation of the first hydration shell of protonated glycine.

Ab initio investigation of the first hydration shell of protonated glycine Zhichao Wei, Dong Chen, Huiling Zhao, Yinli Li, Jichun Zhu, and Bo Liu Citation: The Journal of Chemical Physics 140, 085103 (2014); doi: 10.1063/1.4862985 View online: http://dx.doi.org/10.1063/1.4862985 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/8?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Glycine in aqueous solution: solvation shells, interfacial water, and vibrational spectroscopy from ab initio molecular dynamics J. Chem. Phys. 133, 114508 (2010); 10.1063/1.3481576 Folding processes of the B domain of protein A to the native state observed in all-atom ab initio folding simulations J. Chem. Phys. 128, 235105 (2008); 10.1063/1.2937135 Hybrid ab initio Kohn-Sham density functional theory/frozen-density orbital-free density functional theory simulation method suitable for biological systems J. Chem. Phys. 128, 014101 (2008); 10.1063/1.2814165 Hydration of ion-biomolecule complexes: Ab initio calculations and gas-phase vibrational spectroscopy of K + ( indole ) m ( H 2 O ) n J. Chem. Phys. 124, 184301 (2006); 10.1063/1.2191047 Ab initio molecular dynamics study of glycine intramolecular proton transfer in water J. Chem. Phys. 122, 184506 (2005); 10.1063/1.1885445

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 138.102.119.226 On: Tue, 29 Apr 2014 10:10:54

THE JOURNAL OF CHEMICAL PHYSICS 140, 085103 (2014)

Ab initio investigation of the first hydration shell of protonated glycine Zhichao Wei, Dong Chen,a) Huiling Zhao, Yinli Li, Jichun Zhu, and Bo Liua) Institute of Photo-Biophysics, Physics and Electronics Department, Henan University, 475004 Kaifeng, China

(Received 24 August 2013; accepted 10 January 2014; published online 28 February 2014) The first hydration shell of the protonated glycine is built up using Monte Carlo multiple minimum conformational search analysis with the MMFFs force field. The potential energy surfaces of the protonated glycine and its hydration complexes with up to eight water molecules have been scanned and the energy-minimized structures are predicted using the ab initio calculations. First, three favorable structures of protonated glycine were determined, and the micro-hydration processes showed that water can significantly stabilize the unstable conformers, and then their first hydration shells were established. Finally, we found that seven water molecules are required to fully hydrate the first hydration shell for the most stable conformer of protonated glycine. In order to analyse the hydration process, the dominant hydration sites located around the ammonium and carboxyl groups are studied carefully and systemically. The results indicate that, water molecules hydrate the protonated glycine in an alternative dynamic hydration process which is driven by the competition between different hydration sites. The first three water molecules are strongly attached by the ammonium group, while only the fourth water molecule is attached by the carboxyl group in the ultimate first hydration shell of the protonated glycine. In addition, the first hydration shell model has predicted most identical structures and a reasonable accord in hydration energy and vibrational frequencies of the most stable conformer with the conductor-like polarizable continuum model. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4862985] I. INTRODUCTION

Water can form small active clusters or macro-scale assemblies involved in protonation and micro-hydration,1, 2 which are fundamental chemical reactions and the most important factors in determining structures, properties, and functions of biological molecules.3–11 For example, water can act as proton donor and acceptor in peptide solvation, which is suggested as the substantial driving force in peptide folding. During those folding procedures, amino acids in peptide can form hydrogen bonds with the surrounding water molecules.12, 13 Hydrogen bonds further act as one of the most important noncovalent interactions and provide a key to determine the structures and properties of macromolecules in living systems.14–21 Therefore, investigating the interactions between amino acids and water is an essential step to understand the solvation process and the biological activity of protein, the structure of polymers, proteins, and nucleic acids, which allows opening insights into diverse fields, such as protein engineering, molecule recognition, and drug design. Experimental contributions have assessed the properties of hydrated complexes of charged amino acids and small peptides. Strategies such as electrospray high-pressure mass spectrometer (MS)22–26 and spectroscopy techniques27–34 are able to study multiple hydrations of amino acids and hydration dynamics. Gao et al. studied the hydration of the pro-

a) Authors to whom correspondence should be addressed. Electronic

addresses: [email protected] and [email protected].

0021-9606/2014/140(8)/085103/10/$30.00

tonated aromatic amino acid phenylalanine, tryptophan, and tyrosine with up to five water molecules, in which they found that both the ammonium and carboxyl groups offered good water-binding sites.35 Kamariotis et al. studied the infrared spectra of protonated and lithiated valine in hydration processes including up to four water molecules.36 These studies have investigated the so-called specific hydration of the favorable water binding sites. However, due to experimental limitations such as the pumping requirements of the MS vacuum system, the degree of hydration (the upper bound of water uptake) is limited. As a result, the total water solution effect is not able to be observed. In order to verify and supplement the experiment results, extensive theoretical methods have been applied to simulate the aqueous solution environment. The continuum models such as Onsager and the polarized continuum model (PCM) are generally used to treat the environment, but their accuracy strongly depends on the empirical parameters. For example, the radius for the “united atom” which is used to create the cavity for the PCM model has a great bearing on its calculation. Particularly, the continuum model is shown to be inadequate when describing the strongly interacting Hbonding solvent such as water.37–39 Even, some first principle molecular dynamic simulations such as Born-Oppenheimer molecular dynamics (BOMD)40 and Carr-Parinello molecular dynamics (CPMD)41, 42 have been employed to model the structural properties and determine the hydration shell of biological molecules in aqueous solution. However, their requirement for huge computation ability always should be considered when involving these first principle calculations. Moreover, the molecular dynamics simulation relies on

140, 085103-1

© 2014 AIP Publishing LLC


085103-2

Wei et al.

empirically parameterized potential functions and a molecular mechanics force field for determining the structure.43, 44 Recently studies of the micro-hydration of biological molecules based on suitable building rules and quantum mechanics has been performed focusing on water solution effect on the structure and the properties of biological molecules. Usually, the micro-hydration study of biological molecules relies on chemical intuition or a conformational search to initialize structures and the building rules are crucial to the validity of a hydration shell. Michaux et al.45–49 described the micro-hydration of protonated amino acids by a family tree, which assessed the existence of a Darwinian-like logic to build protonated glycine complexes. In their study, a conformer with up to four surrounding water molecules was finally built; however, it was not adequate to fully saturate the protonated amino acid. With the purpose of simulating the total water solution effect, the first hydration shell structure of GlyH+ · (H2 O)1-N were built in this paper. Here the N is the number of water molecules saturating the first hydration shell of GlyH+ . In our work, the first hydration shell only takes into account water molecules directly interacting with the GlyH+ , by which, the so-called specific hydration process can also be studied. In the real aqueous environment all the hydration sites are hydrated by the surrounding water molecules. Thus such a saturated first hydration shell is proposed to focus the computational resources on the most interesting structures. In addition, a systematic comparison of the conformational structure, the hydration energy, and the vibrational frequency determined by the first hydration shell model and the conductor-like polarizable continuum model (CPCM) was made to assess the accuracy of this method to mimic the water effect on the protonated glycine. The CPCM model was developed from PCM model in which the cavity is created via a series of overlapping spheres. Particular attention is devoted to large systems requiring suitable iterative algorithms to compute the solvation charges. Though, a simulation of the first hydration shell of the deprotonated glycine had been made in our previous work,50 the different charge states of the GlyH+ can lead to different water binding sites, which may result in distinctive hydration shell structure. Moreover, some new physical insights have been made in this paper. This present work combined with our previous work on the deprotonated glycine could shed light on the understanding the water effect on the zwitterionic glycine.

II. COMPUTATIONAL METHODS

The initial structures of the protonated glycine (GlyH+ ) were generated by a 1000-step conformational search (CSearch) which was performed by the Mote Carlo multiple minimum (MCMM) method and MMFFs force fields, as implemented via the MacroModel 9.5 software (Schröinger, LLC, New York, 2007). During the Csearch, all the rotatable bonds (i.e., the C–N bond of the ammonium group, the C–O bond of the carboxyl group, and the C–C bond) were included and all the torsion rotation values were set

J. Chem. Phys. 140, 085103 (2014)

as 0◦ –180◦ , 100 steps per rotatable bond. The energy window for saving structures was set as 50 kJ/mol. The resulting structures were examined, from which five representative structures were chosen for high-level geometry optimization and frequency analysis. The calculations were performed by the Gaussian-03 program at the second-order MøllerPlesset (MP2) level. The extended 6-311++G(d,p) atomic basis set including both polarization and diffuse functions was applied. Three global minimum conformers were generated and they are shown in Figure 1 and named as GA, GB, and GC. The computational details for searching the lowestenergy conformers of GlyH+ · (H2 O)n are as follows. As shown in Scheme 1, the computational procedure consists of 4 stages and 3 operations. In stage 1, one starting structure of GlyH+ · (H2 O)n is generated by randomly adding n water molecules to GA, GB, and GC. Stage 1 was followed by operation (i), in which a 10 000-step Csearch was performed. Hundred steps per rotatable bond and energy window for saving structures of 50 kJ/mol were also defined. During the Csearch, all water molecules were allowed to move freely (translation of 0–5 Å and rotation of 0◦ –180◦ ). In order to prevent the conformers GA, GB, and GC from transition to each other, a constraint of the dihedral angle and torsion rotation of GlyH+ was defined, which allows the geometries of GlyH+ to adjust under MMFFs force field without the transition. In this way, we can avoid repetitive calculations. After operation (i), enormous resulting conformers are located (shown in stage 2). Since in operation (i) similar structures and mirror structures can be generated, operation (ii) was introduced to screen them and then the conformers with water molecules in the first hydration shell (the first shell structure) were identified. Here the first-hydration-shell water molecules are defined as those interacting with GlyH+ directly through hydrogen bonds (threshold value: maximum distance 2.5 Å, minimum donor angle 120◦ , mini acceptor angle 90◦ ). As shown in stage 3, conformers with all water molecules located in the first hydration shell were chosen after operation (ii). Then operation (iii) was applied to these conformers, in which the MP2/6-311++G(d,p) level geometry optimization and harmonic vibrational frequency analysis were performed to characterize the stationary points. All thermodynamic functions were computed under T = 298.15 K and P = 1012.95 mbar using the standard thermochemistry model implemented in Gaussian 03. All energies were corrected by the zero point error. The basis set superposition errors (BSSEs) were corrected using the wellestablished counterpoise approach. In quantum chemistry, calculations using finite basis sets are sensitive to BSSE. In the counterpoise method (CP), the BSSE is calculated by reperforming all the calculations using the mixed basis sets, and the error is then subtracted a posteriori from the uncorrected energy. The BSSE corrections were applied not only to the internal energy but also to the geometry optimization and vibrational frequency calculations.51 If not noted, all the values in this paper are BSSE-free, which is often required when describing hydrogen bond complexes.45–49, 52 The optimized conformers with similar energy and geometry were classified and one of the conformers in one group was


085103-3

Wei et al.

J. Chem. Phys. 140, 085103 (2014)

FIG. 1. Part of the representative conformers and their relative energies (kJ/mol) of GlyH+ · (H2 O)1–4 . The dashed lines indicate the possible hydrogen bonds or non-covalent interactions.

chosen as the representative conformers shown in the stage 4. For example, the starting structure for the GA · (H2 O)2 was achieved by randomly adding two water molecules to conformer GA. After operation (i), 210 conformers were generated. Then 16 conformers with first-hydration-shell water molecules were picked up and 5 representative structures were chosen after the geometry optimization and classifica-

tion. It should be noted that the resulting conformers could be increased substantially when more water was attached to the complexes after operation (i) and therefore the high level MP2 calculations in operation (iii) could be tremendously time consuming. This work also attempts to find a way to save computational resources based for studying the hydration process.


085103-4

Wei et al.

J. Chem. Phys. 140, 085103 (2014)

1. GA · (H2 O)

SCHEME 1. The computational procedure for searching the lowest-energy conformers of GlyH+ · (H2 O)n . (i) MCMM conformational search, (ii) identifying the first hydration shell water-GlyH+ complexes, and (iii) MP2 optimization and classification.

III. RESULTS AND DISCUSSION A. Structures of protonated glycine

The three conformers of GA, GB, and GC are included in Figure 1. Clearly, the different spatial structure and Cartesian coordinates of the carboxyl are the main difference among these three conformers. In GA and GB, the hydrogen atom of the carboxyl group interacts with the oxygen atom of the carbonyl group. While in GC, without such an interaction the energy of GC increases to 43.5 kJ/mol, which indicates that this noncovalent interaction chiefly contributed to the stability of GA and GB. For them all, the hydrogen of the ammonium group forms a strong hydrogen bond with the oxygen of the carboxyl group, with the bond distance of 2.06, 2.28, and 2.45 Å, respectively. The hydrogen bond distance of GB is longer than that of GA, which indicates a weaker hydrogen bond in GB. Consequently, GB is less favorable in energy than GA by 12.3 kJ/mol.

B. The hydrated complexes

To help the readers distinguish the hydrated complexes, a systematic nomenclature composed of three parts is applied. Taking GA3a as an example, the capital letter GA designates the type of GlyH+ and the following number indicates the number of water molecules associated with the protonated glycine; the conformational energies are ordered by appending the lowercase letters of a, b, c, etc., for example, the letter a in GA3a stands for the lowest-energy conformer of GA trihydrate. The building process of GA hydrated complexes is shown as follows. Parts of the conformers of GA complexes are shown in Figure 1 in ball and bond type. All the representative conformers for GA, GB, and GC complexes and their Cartesian coordinates are shown in the supplementary material.59

In Figure 1, the water molecule in GA1a and GA1b are both attached to the ammonium group. As for GA1c, the sole water molecule interacts with the carboxyl group. The energy of GA1b is only 0.6 kJ/mol higher than GA1a, while the energy of GA1c is 10.3 kJ/mol over GA1a. This significant distinction of the relative energy indicates that the first water molecule prefers to hydrate the charged ammonium group. Moreover, the water molecule in GA1a forms a strong hydrogen bond with the ammonium group with the length of 1.66 Å and a weak interaction with the oxygen of the carboxyl group. This water molecule provides a bridge between carboxyl and amino groups, which is defined as bridging water. In the deprotonated glycine monohydrate, bridging water between the carboxylate group and amino group is also found in the conformer G1c.50 On the contrary, G1c is a very unstable conformer whose relative energy is 12.3 kJ/mol. Clearly, a different charge state is usually accompanied by remarkably distinctive water binding sites. 2. GA · (H2 O)2

The most-stable conformer GA2a has one bridging water molecule and one water molecule in front of the ammonium group. The second-stable conformer GA2b has two bridging water molecules and it is 0.2 kJ/mol higher than GA2a in energy. Both conformers GA2c and GA2d have a water molecule binding to the hydrogen of the hydroxyl group and their relative energies are higher than GA2a and GA2b by about 6.9 kJ/mol and 6.7 kJ/mol, respectively. Obviously, GA2c and GA2d are more unstable than GA2a and GA2b. It is indicated that even with two water molecules, the hydration site of the ammonium group still has a strong competition against that of the carboxylic group. However, this value has decreased from 9.7 kJ/mol for the monohydrate to 6.7 kJ/mol for the dihydrate, which indicates a changing competition among different hydration sites. It is interesting to note that the conformer GA2e has a bridge made of two water molecules which are coplanar due to a weak interaction between the water molecule and the oxygen of the carbonyl group. This conformer is nearly a second shell structure. Conformer GA2e reported here corresponds to a bridging minima structure, which is dominant in zwitterionic glycine.53 Completely contrary to the zwitterionic glycine, the protonated glycine-water complexes with a bridge consisting of two water molecules is significantly higher in energy. Conformers GA3f (13.9 kJ/mol higher than GA3a), GA3g (14.7 kJ/mol higher than GA3b), and GA4d (10.0 kJ/mol higher than GA4a) presented in this article are also structures with a bridge consisting of two water molecules. 3. GA · (H2 O)3 and GA · (H2 O)4

In GA3a, the amino group is bonded with three water molecules and two of them are bridging water molecules. Both GA3b and GA3c have a water molecule attached to the carboxylic moiety and their relative energies are 3.7 and


085103-5

Wei et al.

4.7 kJ/mol higher than GA3a, respectively. The results also indicate that the ammonium group is not saturated by two water molecules. The third water molecule still prefers to bind with the cationic group. It still has a competition against the site of the carboxylic group. Conformer GA3e attracts our attention with a water cluster that connects the ammonium group and the carboxylic group. It is encouraged to be a potential energy-minimum structure because water molecules prefer to form cluster via its excellent talent in constructing hydrogen bonding network. For the tetrahydrate, the most stable conformer GA4a has three water molecules attached to the ammonium group and one water molecule attached to the carboxylic group. As the energy will increase if the fourth water locates at other sites (see GA4b and GA4c), it causes that the fourth water molecule prefers to interact with the hydrogen of the carboxylic group. The results suggest that, with three water molecules binding to the strong positive charged ammonium group, the charge seems to be blocked and thus the carboxyl group becomes competitive to attach the fourth water molecule. The ammonium group can strongly attract three water molecules. The experiment operated by Wincel54 showed that the first three water molecules bond to the hydrogen atoms of the ammonium group. Similarly, Kamariotis et al.36 reported that the first three water molecules seem to preferentially bind to the protonated ammonium group in protonated valine. Using infrared photodissociation (IRPD) spectroscopy and computational chemistry, Prell et al.55 found that a substantial population of protonated phenylalanine without water molecules attached to the carboxylic group for ions of PheH+ (H2 O)1–3 . Our calculations are reasonably consistent with these experimental results. Since energies of all the conformers in this paper are calculated at the temperature of zero, it is necessary to compare them with the Gibbs free energies in Michaux’s work45 using the standard thermochemistry model (T = 298.15 K). All conformers of GlyH+ · (H2 O)n = 0–4 predicted in our work, including their structures, hydrogen bond networks and the corresponding relative G values shown in Table I are in good agreement with those determined by Michaux’s work. In Table I, the orders of relative energy are slightly different from that of Gibbs free energy. For example, GA1a is the lowestenergy conformer, but it is the second lowest in Gibbs free energy; GA3e, which is the 5th of seven low-energy conformers, turns to be the highest conformer in Gibbs free energy; GA4b, the second lowest conformer in energy, ranks 4th lowest according to Gibbs free energy. In other words, at high temperature these conformers become a bit unstable. This could be attributed to the effect of the bridging water and the secondshell solvation structure. At higher temperature, the experiments in Michaux’s work show a low ratio of the bridging water structure and the conformers with the second-shell waters become unstable in Gibbs energy. From analyzing the conformers of GA1a, GA3e, and GA4b, GA1a has the structure of the bridging water, while GA3e and GA4b belong to the conformers with second-shell water in the definition of Michaux’s work, which is different from our definition of the second shell in this work. Therefore, the calculations show

J. Chem. Phys. 140, 085103 (2014) TABLE I. Relative Gibbs (G) values in kcal/mol of GlyH+ · (H2 O)n = 0–4 taken from this work and Ref. 45. Labela

Ga

Gb

Labelb

n = 0 GA GB GC

0.00 0.00 A 3.58 3.58 B 9.84 8.37 C

n = 1 GA1b GA1a GA1c

0.00 0.00 A1 0.39 0.36 A2 2.36 2.02 A3

n = 2 GA2a GA2b GA2c GA2d GA2e

0.00 1.46 1.73 1.76 5.18

0.00 1.39 1.69 1.74 5.13

A32 A22 A21 A31 A22b

Labela n=3

Ga

Gb

Labelb

GA3a 0.00 0.00 A322 GA3b 0.07 0.01 A321 GA3c 1.95 0.91 A221 GA3d 3.19 4.26 A222b GA3f 4.28 3.16 A313b GA3g 5.56 5.53 A212b

GA3e n = 4 GA4a GA4c GA4d GA4b GA4e GA4f

5.91 0.00 2.56 2.84 3.73 5.05 8.09

5.91 0.00 2.55 2.82 3.72 5.04 4.78

A222b A3221 A3222b A3212b A3222b A2212b A2212b

a

The label and G values are taken from the frequency calculation in this paper. The label and G values in Ref. 45. The italics are for second (and third)-shell solvation in Ref. 45.

b

the Gibbs free energies for the three specific conformers are higher and also confirm the conclusions drawn from the results in Michaux’s work. Let us now pay attention to the structures and the relative energies of conformers GA · (H2 O)1–4 that will be helpful in investigating the building process of the first hydration shell. First, we can conclude that there are three kinds of “water bridges.” The first one is the bridge made of one water molecule, such as GA1a; the second one is the bridge consisting of two water molecules, such as GA2e; the third one is the bridge including three water molecules, such as GA3e. The aforesaid three kinds of bridge including one, two, and three water are named B1, B2, and B3, respectively. In B1, the water molecule forms a strong hydrogen bond with the ammonium group but a weak interaction with the oxygen of the carbonyl group. In B2, one water molecule is bonded to the ammonium group tightly, just like the water molecule of B1, while the other one locates between the former water molecule and the oxygen of the carbonyl group. Similarly, water molecules network connects the carbonyl group and the ammonium group in B3 (GA3e). The water molecules of B1, B2, and B3 insert into GA between the ammonium group and the carboxyl group to create the co-operatively linked hydrogen bonds, which helps to stabilize the whole complexes. By analyzing the distance of hydrogen bonds among the bridging water, the ammonium and the carbonyl groups, we find that the weak interaction (bond L) between the bridging water and the oxygen of the carbonyl becomes gradually stronger. The bond distance (Table SIII in the supplementary material)59 of bond L continuously decreases from 3.2 to 2.0 Å and the bond angle changes from 92◦ to 142◦ , which indicates that the hydrogen bonding network is enhanced with the increasing hydration water molecules. Especially, the hydrogen bonding becomes remarkably stronger when more than five water molecules are attached. The bond angle of bond L increases sharply from 105◦ to 141◦ . This is due to the formation of the bridge B3 in these conformers, in which a strong hydrogen bond is formed between the bridge water and the oxygen of


085103-6

Wei et al.

FIG. 2. Optimized structure of protonated glycine conformer GA with possible hydration sites. Site I consists of the bridging site I1 and the site in front of the ammonium group I2.

the carbonyl group. Though B3 is not a dominant structure when the number of hydrate water is not more than four, with the hydration process B3 has grown to be a dominant structure for protonated glycine existing in the lowest-energy conformers of GA5a, GA6a, and GA7a. The relevant role water showing is to promote association between two groups through its excellent hydrogen donor and acceptor. In fact, water is directly involved in protein-DNA recognition where are highly solvated. Bridging water often plays a vital role in molecular recognition.56 Second, it is meaningful to investigate the incremental hydration process and figure out the mechanics of how are water molecules arranged their positions. As can be seen in Figure 1, the water molecules in the lowest-energy conformer GA · (H2 O)1–3 locate at the bridging site or in front of the ammonium group. For monohydrate of GA, GA1a is more stable than GA1b by energy of 0.6 kJ/mol which means the bridging site is competitive than the site in front of the ammonium group. While for the dihydrate of GA, GA2a, who has a water

J. Chem. Phys. 140, 085103 (2014)

molecule at the bridging site and a water molecule in front of the ammonium group, is more stable than GA2b which is attached by two bridging water molecules. So, the sequentially added water molecule binds to the alternative site in front of the ammonium group and for the trihydrate of GA the third water molecule (GA3a) bonds to the bridging site again. We can conclude that the bridging site and the site in front of the ammonium group are alternatively hydrated by the first three water molecules. In case of GA4a, the carboxyl group turns to a new competitive water-binding site and the fourth water molecule was attached to the carboxyl group, which also supports our idea. Similarly, for the deprotonated glycine50 when water molecules are attached to glycine anion, the bridging position appears increasingly competitive in hydration process. In a word, the competition capability of different hydration sites is alternatively changing with the sequential hydration process and the hydration dynamics is driven by saturation of the competitive hydration sites. The above analysis shows that the hydrate sites of GA fall into two typical categories: site I (around the ammonium group) and site II (the carboxyl group), as shown in Figure 2, where the bridging site is labelled as site I1 and the site in front of the ammonium group is defined as site I2. As we have discussed, water molecules prefer to occupy the hydrate site I rather than the hydrate site II. Here, a deep insight into the dynamic competition process has been investigated on basis of the hydration energies and interaction energies of site I and site II. The hydration energies are calculated by E hyd = EGA·(H 2 O)n − EGA·(H2 O)n−1 − EH2 O whereEGA·(H2 O)n is the energy of the is calcuoptimized structures of GA · (H2 O)n . EGA·(H 2 O)n−1 lated by removing one water molecule from the optimized conformer of GA · (H2 O)n . The most stable conformers with water molecules at the corresponding hydrate site are used to calculate the hydration energy. Plots of the hydration energies of GA · (H2 O)n=1–4 as a function of the number of water molecules are shown in Figure 3. Initially, the hydration energies of the conformers with water molecules at hydrate site I1 and site I2 are higher than that of the hydrate site II.

FIG. 3. Plots of hydration energy and interaction energy of hydration complex of GA have been addressed as a function of the number of adding water molecules.


085103-7

Wei et al.

J. Chem. Phys. 140, 085103 (2014)

FIG. 4. Low-energy conformers of GA · (H2 O)5 and GA · (H2 O)6 and their relative energies (kJ/mol).

However, the hydration energies of site I1 and site I2 decrease more rapidly than hydrate site II with the increase of the hydrate water. The three lines show a noticeable trend to approach with each other in the process of sequentially adding water molecules to GA. The interaction energies of the three hydrate sites have also been calculated using the counterpoise approach. The interaction energy is caused by the water molecules at one hydration site with its remaining parts. The equation Einteraction = Ecorrected − (Ewater + Eremain ) is introduced to calculate the interaction energy. Ecorrected is the BSSE corrected total energy of the GA-water complex. Ewater and Eremain are the energies of the water molecules and the remaining parts. The corresponding conformers for site I1 are those with all the water molecules being attached to site I1, so conformers GA1a, GA2b, GA3e, and GA4b are chosen for calculation. For site I2, conformers GA1b, GA2a, GA3a, and GA4a are chosen. For site II, GA1c, GA2c, GA3b, and GA4a are chosen. As is shown in Figure 3, the plots of the interaction energies per water molecule of GA · (H2 O)n=1–4 as a function of the number of water molecules are presented. It is clear that the lines of site I1, I2 and site II also have the trend to approach to lower region. For the tetrahydrate, the interaction energies of site I1 and site II are nearly the same. For the first hydration shell of the deprotonated glycine, the hydration sites also fall into two typical categories.50 Obviously, water is preferential to hydrate the charged part rather than the other water-binding sites at the beginning of the hydration process. For the deep hydration, different hydration sites trend to equally bind water molecules.

to the –OH of the carboxyl group and other water molecules are attached to the hydration site I. It suggests that site II is saturated by one water molecule and the sequentially added water molecule prefers to stay at site I1. For instance, in conformer GA6a, the ammonium hydrogen atom H(2) attaches two water molecules. Only one favorable representative structure of GA · (H2 O)7 is found after the Csearch and the optimized structure is shown in Figure 5. In GA · (H2 O)7 , the site II is hydrated by one water molecule and the bridge B3 forms at the bridging site. The ammonium hydrogen atoms H(2) and H(3) both attach two water molecules while H(1) attaches only one. Such a configuration of first-shell water around hydrate site I is ever found in GA6c, which is 6.8 kJ/mol higher than GA6a. Generally speaking, the cost to calculate the energies and structures of glycine-water complexes increases extensively when more water molecules are attached to glycine. Especially, all water molecules are allowed to move freely during the Csearch in our computational method. In Table II, the quantities of conformers GA · (H2 O)1–7 in each stage of the building procedure are given. For GA · (H2 O)1 the number of resulting structures by Csearch is only 43 while for GA · (H2 O)6 it is 6927. This number rises rapidly

4. GA · (H2 O)5–7

Figure 4 shows the possible conformers of GA · (H2 O)5 and GA · (H2 O)6 in which only one water molecule is bonded

FIG. 5. Structure of the GA · (H2 O)7 . The first hydration shell of GA is saturated by seven water molecules.


085103-8

Wei et al.

J. Chem. Phys. 140, 085103 (2014)

TABLE II. The number of conformers in each stage.

GA · (H2 O)n n=1 n=2 n=3 n=4 n=5 n=6 n=7 a b

Resulting structures

Initial structuresa

Possible conformationsb

43 210 563 2655 4609 6927 7946

10 16 29 15 6 5 1

3 5 7 6 3 3 1

Initial structures for MP2 optimization, i.e., the GA-first hydration shell complexes. The representative conformations after high level optimization.

with the increase of the hydrate water molecules shown in Table II. Here we focus on simplifying the calculation process. Compared GA2a to GA1a, it seems that GA2a is built by just adding one water molecule to GA1a. Following this pattern, GA · (H2 O)n≤7 can be considered by adding one water to the lowest conformer of GA · (H2 O)n–1 . Moreover, the “added” water molecule always binds to the competitive sites of GA · (H2 O)n–1 , which could be understood by alternative hydration of competitive sites discussed above. For example, in GA6a, the ammonium hydrogen atoms H(3) attaches two water molecules and H(2) is clear to be an alternative competitive site for the next added water. So it seems that GA · (H2 O)7 can be obtained by adding one water molecule to the competitive hydration site H(2) of GA · (H2 O)6 . 5. GA · (H2 O)8

Thus, a simplified Csearch method can be proposed. First, the starting structure of GA · (H2 O)n can be generated by adding one water molecule to lowest conformer of GA · (H2 O)n–1 . Second, during the Csearch process, the core geometry of GA · (H2 O)n–1 is unchanged, while only the one added water is allowed to move freely to perform the Csearch, resulting in less initial structures and time costs. In this simplified Csearch method, no stable structures with water molecules in the first hydration shell were found on the PES of GA · (H2 O)8 . Moreover, according to the alternative hydration process, we established several starting structures for GA · (H2 O)8 by adding one water molecule to GA · (H2 O)7 at different hydrate sites. After the structural optimization, the added water molecules moved to the second hydration shell. It is deduced that the first hydration shell of GA is formed by seven water molecules.

ber of water molecules in the saturated first hydration shell is seven and the bridging water, bridges B1, B2, and B3, are all found in the hydrated complexes of GC. But for GB, six water molecules are enough to saturate the first hydration shell, and there are no bridging water molecules and bridges such as B1, B2, and B3. It could be understood by the different spatial structures of the ammonium and carbonyl group in GA, GB, and GC. In GA and GB, the ammonium group faces the oxygen atom in the carbonyl group of –COOH (the “trans” form), while in GB, the ammonium group faces the oxygen atom on –OH group of COOH (the “cis” form). The oxygen of the hydroxyl is not as a good hydrogen-bond acceptor as the oxygen in the carbonyl group of COOH. As a consequence, for GB, water molecules binding to the ammonium group have no interaction with the oxygen of the hydroxyl group, which results in a lack of bridging water molecules. Even though, there are obvious similarities in the hydration processes of GA, GB, and GC. For example, for GB and GC, water molecules also prefer to hydrate the ammonium group rather than the carboxyl group at the beginning of the hydration process. Moreover, the hydration processes of GB and GC are driven by the competition between different hydration sites. 2. Water effect

Though GA is the lowest-energy conformer in the three possible structures of protonated glycine, the energy difference among the complexes of the three conformers may change with the increase of hydrated water molecules. In order to clearly illustrate the role of water in the micro-hydration process, the relative energies between hydrate complexes of GB, GC with GA are studied. As shown in Figure 6, the energies of GB and GC are about 15 and 39 kJ/mol higher than that of GA, respectively. With one water molecule hydrated, the relative energies of complexes for GB and GC decrease to 14.1 and 34.8 kJ/mol, respectively, and the values continuously decrease with the hydration process. When GA, GB, and GC are fully saturated by the first hydration shell, the values decrease to 6.2 and 17.9 kJ/mol, respectively. Especially, the relative energies of the hydrate complexes of GC

C. Comparison between GA, GB, and GC

1. Hydration shell

The first hydration shell of GB and GC has been built up by performing the same computational methods as GA. Finally, GB and GC are fully hydrated by six and seven water molecules, respectively. Until now, the first hydration shell of the protonated glycine has been built on the base of the three initial geometries of GA, GB, and GC. The hydration process of the GC is similar with that of GA. For GC, the num-

FIG. 6. Relative energies between hydrate complexes of GB, GC with GA.


085103-9

Wei et al.

decrease sharply. Such a result shows that water has more influence on stabilizing the more unstable molecules. Indeed, water is an excellent component at the interface of biomolecular complexes. It can stabilize and facilitate biomolecules via its hydrogen-bond donor and acceptor. In addition, the relative energies of GA, GB, and GC are calculated using the CPCM model. The energies, which are 6.7 and 16.6 kJ/mol for GB and GC, respectively, are in a good agreement with the results of the first hydration shell model. The saturated first hydration shell model has an approximately similar effect with the continuum dielectric medium model in stabilizing the molecule. D. Comparison of conformational structure, hydration energy, and vibrational frequency determined by the first hydration shell model and CPCM model

1. Conformational structure

We have investigated the structure of GAcpcm and the structure of GlyH+ in GA · (H2 O)7 . The structural parameters are presented in Table SI in the supplementary material showing that the structures of them are nearly identical.59 However the differences of the hydration energy and the vibrational frequency are obvious. 2. Hydration energy

We have to mention that the BSSE correction is not performed using the CPCM model. The difference of the total Gibbs energies between MP2 and CPCM-MP2 calculations for the gas-phase optimized structures was used to estimate the CPCM hydration energy in water. The value for GA is 79 kcal/mol. Such a huge amount is due to the cationic nature of the protonated glycine which is strongly stabilized in the continuum dielectric medium. The hydration energies for each additional water molecules in our first hydration shell model are calculated. The data are provided in Table SII in the supplementary material.59 By summing them, the hydration energy is 33.7 kcal/mol.

J. Chem. Phys. 140, 085103 (2014) TABLE III. Asymmetric stretching vibration of –NH3 and stretching vibration of C=O obtained by the saturated first hydration shell model and the CPCM model in the gas phase. The values of the un-solvated GA are also given to compare. All calculated frequencies were scaled by a factor of 0.95.

vas (NH3 ) vs (C=O)

GAunsolvated

GACPCM

GA · (H2 O)7

3313 1737

3180 1703

3258 1725

shifted by 55 cm−1 with respect to un-solvated GA, but the value by the CPCM model is much more strongly redshifted to lower wave numbers by 133 cm−1 . The redshifted effect on the C=O vibration is 34 cm−1 by the CPCM model and only 12 cm−1 shifted to the lower wavenumber by the saturated first hydration shell model. Third, two types of asymmetric stretching modes of N–H are separated by 56 and 27 cm−1 in un-solvated GA and GA · (H2 O)7 , respectively. While in the solvated GA by the CPCM model, these two types asymmetric stretching mode of NH3 both locate at 3180 cm−1 . This may suggest that the CPCM model is not enough to describe the strongly specific interaction between water in the first hydration shell and GA. The IR spectrums calculated from structure of solvated GA by the CPCM model and the saturated first hydration shell GA · (H2 O)7 are presented in Figure SI in the supplementary material where they can compare to the un-solvated GA.59 From the above, the structures of GA determined by the saturated first hydration shell model and the CPCM model are nearly same while the hydration energy and the vibrational frequency are different. Our first hydration shell model predicts lower hydration energy and smaller vibrational redshift values by water effect. This might result from the lack of the long range interaction of surrounding water in our first-shell water model since the second and the higher hydration shells are not included in our model.

IV. CONCLUSIONS

3. Vibrational property frequency

To our knowledge, the vibrational properties of amino acids are highly sensitive to the rearrangement of the neighboring water molecules in the first hydration shell.38 Indeed, the C=O and N–H stretching modes of the amide linkage have been widely used to characterize the structure of peptide in solution.57, 58 Here the asymmetric stretching vibration of N-H3 and the stretching vibration of C=O are presented. The values of GA obtained by the saturated first hydration shell model and the CPCM model in the gas phase are listed in Table III compared with the values of un-solvated GA. The first conclusion from the data is that the water effect caused the vibration of C=O and N-H3 shifted to lower wave numbers. The reason is that the bond strength of C=O and N-H3 is greatly reduced in aqueous solution or solvated by water. Second, in the case of the NH3 group, the band of asymmetric stretching vibration in GA · (H2 O)7 is located at 3258 cm−1 by the saturated first hydration shell model; it is red-

Using the MCMM conformational search analysis and the MP2/6-311++G(d,p) approach combined with a full BSSE correction, GlyH+ · (H2 O)n = 1–8 complexes have been investigated step by step to simulate the saturated first hydration shell of the protonated glycine. The first three water molecules strongly bind to the hydrogen-bond donors of the ammonium group and while only one water molecule binds to the hydrogen of the carboxyl group. Five water molecules can fully pack the positive ammonium group and seven water molecules can fully saturate the first hydration shell of GA. We have confirmed that protonation and deprotonation result in explicitly different first hydration shell. Final, we studied the water stabilizing effect and compared our model with the CPCM model in simulating the bulk water effect. Compared to the CPCM model, our first hydration shell model predicts the low hydration energies and the small vibrational redshift by water effect. The main conclusions and perspectives from this work are following:


085103-10

Wei et al.

(1) This work and the previous study on deprotonated glycine illustrate that the competition between the different hydration sites changes with the hydration process and results in an alternative dynamic hydration process. Based on such an insight we can simplify the building process for the micro-hydration and thus save computational costs. (2) The hydration processes of GB and GC illustrate that water plays an excellent role in stabilizing the unstable molecules. Water forms small hydrogen binding network (e.g., the bridging water) to promote association between two groups, which is often involved in molecular recognition. (3) The first hydration shell model is able to model the specific hydration and further reproduce a considerable tendency of bulk water effect, which is comparable to that of the CPCM model. (4) A combination of our first hydration shell model and the polarized continuum model may simulate both the specific hydration effect and bulk water effect.

ACKNOWLEDGMENTS

We are especially thankful for Dr. Mingdong Dong’s instructive discussion and suggestions. This work is supported by the Henan Province Research Office of Education Fund 122102310096, Henan Province Department of Science and Technology Fund 12A180005 and Henan Natural Science Research Office of Education project (2008A140003 and 2010A140001). 1 M.

F. Chaplin, Nat. Rev. Mol. Cell Biol. 7, 861 (2006). Maihom, E. Schuhfried, M. Probst, J. Limtrakul, T. D. Märk, and F. Biasioli, J. Phys. Chem. A 117, 5149 (2013). 3 V. H. Wysocki, G. Tsaprailis, L. L. Smith, and L. A. Breci, J. Mass Spectrom. 35, 1399 (2000). 4 R. I. Cukier and D. G. Nocera, Annu. Rev. Phys. Chem. 49, 337 (1998). 5 G. J. Kearley, F. Fillaux, M. H. Baron, S. Bennington, and J. Tomkinson, Science 264, 1285 (1994). 6 M. K. Green and C. B. Lebrilla, Mass Spectrom. Rev. 16, 53 (1997). 7 M. F. Jarrold, Annu. Rev. Phys. Chem. 51, 179 (2000). 8 R. H. Wu and T. B. McMahon, J. Phys. Chem. B 113, 8767 (2009). 9 R. H. Wu and T. B. McMahon, J. Am. Chem. Soc. 130, 3065 (2008). 10 R. H. Wu and T. B. McMahon, J. Am. Chem. Soc. 129, 569 (2007). 11 J. Murciano-Calles, E. S. Cobos, P. L. Mateo, A. Camara-Artigas, and J. C. Martinez, Biophys. Chem. 158, 104 (2011). 12 F. Garczarek and K. Gerwert, Nature (London) 439, 109 (2006). 13 A. S. Mahadevi and G. N. Sastry, Chem. Rev. 113, 2100 (2013). 14 M. Feldstein, P. Kireeva, T. Kiseleva, B. Gdalin, M. Novikov, Y. Anosova, G. Shandryuk, P. Singh, and G. Cleary, Polym. Sci., Ser. A 51, 799 (2009). 15 K. E. Feldman, M. J. Kade, T. F. A. de Greef, E. W. Meijer, E. J. Kramer, and C. J. Hawker, Macromolecules 41, 4694 (2008). 16 L. B. Sagle, Y. J. Zhang, V. A. Litosh, X. Chen, Y. H. Cho, and P. S. Cremer, J. Am. Chem. Soc. 131, 9304 (2009). 17 J. L. Li and K. S. Zhao, J. Phys. Chem. B 117, 11843 (2013). 18 D. J. Wolstenholme, J. J. Weigand, E. M. Cameron, and T. S. Cameron, Phys. Chem. Chem. Phys. 10, 3569 (2008). 19 P. Cordier, F. Tournilhac, C. Soulié-Ziakovic, and L. Leibler, Nature (London) 451, 977 (2008). 20 A. Katrusiak, Models, Mysteries and Magic of Molecules (Kruger National Park, South Africa, 2008), p. 209. 21 J. D. Chodera, N. Singhal, V. S. Pande, K. A. Dill, and W. C. Swope, J. Chem. Phys. 126, 155101 (2007). 2 T.

J. Chem. Phys. 140, 085103 (2014) 22 D.

R. Carl, B. K. Chatterjee, and P. B. Armentrout, J. Chem. Phys. 132, 044303 (2010). 23 H. Wincel, Chem. Phys. Lett. 488, 219 (2010). 24 H. Wincel, J. Am. Soc. Mass. Spectrom. 23, 1479 (2012). 25 H. Wincel, J. Phys. Chem. A 111, 5784 (2007). 26 H. Wincel, J. Am. Soc. Mass. Spectrom. 19, 1091 (2008). 27 A. Yamakata, E. Soeta, T. Ishiyama, M. Osawa, and A. Morita, J. Am. Chem. Soc. 135, 15033 (2013). 28 S. Toal, D. Meral, D. Verbaro, B. Urbanc, and R. Schweitzer-Stenner, J. Phys. Chem. B 117, 3689 (2013). 29 D. J. Miller and J. M. Lisy, J. Am. Chem. Soc. 130, 15381 (2008). 30 H. Fricke, K. Schwing, A. Gerlach, C. Unterberg, and M. Gerhards, Phys. Chem. Chem. Phys. 12, 3511 (2010). 31 O. H. Kwon, T. H. Yoo, C. M. Othona, J. A. Van Deventer, D. A. Tirrell, and A. H. Zewail, Proc. Natl. Acad. Sci. U.S.A. 107, 17101 (2010). 32 P. Sasisanker and H. Weingartner, ChemPhysChem 9, 2802 (2008). 33 M. F. Bush, J. S. Prell, R. J. Saykally, and E. R. Williams, J. Am. Chem. Soc. 129, 13544 (2007). 34 N. Derbel, B. Hernandez, F. Pfluger, J. Liquier, F. Geinguenaud, N. Jaidane, Z. Ben Lakhdar, and M. Ghomi, J. Phys. Chem. B 111, 1470 (2007). 35 B. Gao, T. Wyttenbach, and M. T. Bowers, J. Am. Chem. Soc. 131, 4695 (2009). 36 A. Kamariotis, O. V. Boyarkin, S. R. Mercier, R. D. Beck, M. F. Bush, E. R. Williams, and T. R. Rizzo, J. Am. Chem. Soc. 128, 905 (2006). 37 K. J. Jalkanen, I. M. Degtyarenko, R. M. Nieminen, X. Cao, L. A. Nafie, F. Zhu, and L. D. Barron, Theor. Chem. Acc. 119, 191 (2008). 38 E. Tajkhorshid, K. J. Jalkanen, and S. Suhai, J. Phys. Chem. B 102, 5899 (1998). 39 W.-G. Han, K. J. Jalkanen, M. Elstner, and S. Suhai, J. Phys. Chem. B 102, 2587 (1998). 40 I. M. Degtyarenko, K. J. Jalkanen, A. A. Gurtovenko, and R. M. Nieminen, J. Phys. Chem. B 111, 4227 (2007). 41 M. Dracinsky, J. Kaminsky, and P. Bour, J. Phys. Chem. B 113, 14698 (2009). 42 J. Thar, S. Zahn, and B. Kirchner, J. Phys. Chem. B 112, 1456 (2008). 43 M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112, 8910 (2000). 44 P. Cabral do Couto, S. G. Estácie, and B. J. Costa Cabral, J. Chem. Phys. 123, 054510 (2005). 45 C. Michaux, J. Wouters, E. A. Perpete, and D. Jacquemin, J. Phys. Chem. B 112, 2430 (2008). 46 D. Jacquemin, C. Michaux, E. A. Perpète, and G. Frison, J. Phys. Chem. B 115, 3604 (2011). 47 C. Michaux, J. Wouters, E. A. Perpete, and D. Jacquemin, J. Phys. Chem. B 112, 7702 (2008). 48 C. Michaux, J. Wouters, E. A. Perpete, and D. Jacquemin, J. Phys. Chem. B 112, 9896 (2008). 49 C. Michaux, J. Wouters, E. A. Perpete, and D. Jacquemin, J. Am. Soc. Mass. Spectrom. 20, 632 (2009). 50 Y. H. Yao, D. Chen, S. Zhang, Y. L. Li, P. H. Tu, B. Liu, and M. D. Dong, J. Phys. Chem. B 115, 6213 (2011). 51 S. Simons, M. Duran, and J. J. J. Dannenberg, Chem. Phys. 105, 11024 (1996). 52 S. Simon, J. Bertran, and M. Sodupe, J. Phys. Chem. A 105, 4359 (2001). 53 C. M. Aikens and M. S. Gordon, J. Am. Chem. Soc. 128, 12835 (2006). 54 H. Wincel, Chem. Phys. Lett. 439, 157 (2007). 55 J. S. Prell, T. M. Cheang, J. T. O’Brien, and E. R. Williams, J. Am. Chem. Soc. 132, 7811 (2010). 56 C. A. Davey, D. F. Sargent, K. Luger, A. W. Maeder, and T. J. Richmond, J. Mol. Biol. 319, 1097 (2002). 57 M. N. Blom, I. Compagnon, N. C. Polfer, G. von Helden, G. Meijer, S. Suhai, B. Paizs, and J. Oomens, J. Phys. Chem. A 111, 7309 (2007). 58 S. Kumar, A. K. Rai, V. B. Singh, and S. B. Rai, Spectrochim. Acta A 61, 2741 (2005). 59 See supplementary material at http://dx.doi.org/10.1063/1.4862985 for the IR spectrums calculated from the structure of un-solvated GA, solvated GA by the CPCM model and GA · (H2 O)7 . The corresponding structural parameters of GA · (H2 O)7 and GACPCM are listed to compare; the hydration Gibbs free energies for each additional water molecule; the hydrogen bonding distance and the corresponding bond angles between the bridging water with the ammonium group (bond S) and the carbonyl oxygen (bond L); the energy of favorable conformations of hydrated complexes of GA, GB, and GC; the MP2 coordinates, total energy, Counterpoise corrected energy, and zero point energy of all favorable conformations.


Ab-initio SCF investigation ofγ-aminobutyric acid.

Ab initio investigation of the thermal decomposition of n-butylcyclohexane.

Elimination of methane from protonated acetaldehyde: The Ab initio transition state.

Elasticity of grossular-andradite solid solution: an ab initio investigation.

Ab initio investigation of CaO-ZnO alloys under high pressure.

A comprehensive study of cold protonated tyramine: UV photodissociation experiments and ab initio calculations.

Infrared signature of micro-hydration in the organophosphate Sarin: an ab initio study.

Hydration of the cyanide ion: an ab initio quantum mechanical charge field molecular dynamics study.

Hydration of porphyrin and Mg-porphyrin: ab initio quantum mechanical charge field molecular dynamics simulations.

The hydration shell of myoglobin.

Ab initio RNA folding.

Investigation of the hydroxylation mechanism of noncoupled copper oxygenases by ab initio molecular dynamics simulations.

Communication: comparing ab initio methods of obtaining effective U parameters for closed-shell materials.

Ab initio investigation of the structural and electronic properties of amorphous HgTe.

Hydration energies of protonated and sodiated thiouracils.

Investigation of the structure and interaction of cellulose triacetate I crystal using ab initio calculations.

Ab initio investigation of electronic and vibrational contributions to linear and nonlinear dielectric properties of ice.

Structural investigation of protonated azidothymidine and protonated dimer.

Ni3Al(111): an ab initio investigation.

Ab Initio Screening of CO2-philic Groups.

Ab initio MP2 and density functional theory computational study of AcAlaNH2 peptide hydration: a bottom-up approach.

Ab initio calculation of the neutron-proton mass difference.

Ab initio molecular dynamics study of the aqueous HOO(-) ion.

AIDA: ab initio domain assembly server.