- Email: [email protected]
Monte Carlo simulation study on the conformation and interaction of the glycine zwitterion in aqueous solution
THEO CHEM ELSEVIER
Journal of Molecular Structure (Theochem) 397 (1997)113- 119
Monte Carlo simulation study on the conformation and interaction of the glycine zwitterion in aqueous solution Toshio Watanabe, Kazuhiko Hashimoto,
Hideto Takase, Osamu Kikuchi*
Department of Chemistry, University of Tsukuba, Tsukuba 305, Japan
Received 5 June 1996; accepted 24 October 1996
Abstract A Monte Carlo simulation with the statistical perturbation method has been employed to determine the relative stabilities among four representative conformations of the glycine zwitterion (GLYZ) in aqueous solution. The (+,0) = (0,O) and (60,O) conformations, in which the dihedral angle 19between the NCC and CC0 planes is O’, were found to be more stable than the (0,90) and (60,90) conformations, in which 0 =90”. The relative stabilities of these four conformations are well elucidated on the basis of solvent-solute interactions using two-dimensional radial distribution functions. Ab initio MOlGB calculations including the solvent effect via the generalized Born formula gave similar results for the relative stabilities of the four conformations in solution. 0 1997 Elsevier Science B.V. Keywords: Monte Carlo simulation; interactions
Glycine zwitterion;
Conformational
1. Introduction Knowledge of the conformation and interactions of amino acids in water is important for the ultimate understanding of protein hydration and the role of water in biological systems. Glycine is found in the form of a zwitterion in the solid state [I] and in aqueous solution [2]. The conformation of the glycine zwitterion (GLYZ) in aqueous solution has been studied by molecular orbital (MO) calculations in which the solvent effect is included by continuum models [3,4]. Monte Carlo (MC) simulations for the GLYZ in water have been carried out by Roman0 and Clementi [5], Mezei et al. [6], and Alagona et al. [7]. * Corresponding author.
method; Molecular
These simulations concerned the solvation structure for a particular conformation, and revealed a microscopic interaction between the GLYZ and the water molecules. However, little attention has been paid to the conformational energy of the GLYZ in water. Conformational analysis of the GLYZ by MC simulation requires the calculation of the free energy change along the rotation of the C-C bond and/or the C-N bond. There have been no reports of such an analysis. In this study, an MC simulation was applied to four typical conformations of the GLYZ, the (&I!?>= (0, 0), (60,0), (60,90), and (0,90) conformations shown in Fig. 1. Differences in the free energy of solvation were calculated using statistical perturbation theory (WI’) [8,9] and relative stabilities among the four conformations in aqueous solution were predicted.
0166-1280/97/$17.000 1997 Elsevier Science B.V. All rights reserved
PII SO166-1280(96)04942-l
stability; Statistical perturbation
114
T. Watanabe et al./Joumal of Molecular Structure (Theochem) 397 (1997) 113-l 19
Fig. 1. Definition of the two rotation angles, 4 and 8, for the zwitterionic form of glycine: 4 indicates the rotation of the NH3 fragment, while 19 indicates the rotation of the CO* fragment. The conformation shown is (6, 0) = (0,O).
The analysis of the interaction between the GLYZ and solvent molecules introduces special difficulties, since the GLYZ has two hydrogen bonding groups, NHf and CO;, and the calculated distribution functions contain admixtures from various distinctly different modes of solvation [6]. In order to avoid this difficulty, we employed here two-dimensional distribution functions (2D-rdfs) which are represented by two distances between specific atoms. The relative stability of the four conformers in aqueous solution is well elucidated on the basis of the solute-solvent interaction, which is clearly displayed on the 2D-rdf maps.
2. Calculation Two potential functions were used for the water molecules in the MC simulation. One is the threepoint model, which has been determined by ab initio MIDI-4(d) calculations in our previous paper [lo], and the other is the TIP4P model [ll]. The intermolecular potential functions between the GLYZ and Hz0 were newly determined. In this procedure, four conformations were used for the GLYZ, fixing the bond lengths and angles at the standard values [3,7,12]. The MIDI-4(d) calculations were carried out for dimer orientations which were selected in the vicinity of several energy-minimum geometries and for those selected randomly. About 500 configurations were selected for the GLYZ-Hz0 dimer. On the basis of the calculated interaction energies for the dimers, the potential function between the GLYZ and H20, which is constructed by Lennard-Jones (12-6) terms and the Coulomb term,
was determined. The parameters determined are listed in Table 1 and Table 2. The MC simulations were performed in the NVT ensemble according to the standard Metropolis method [13] for one GLYZ solute molecule and 212 water molecules in the cubic cell. The volume of the cell was determined using a density of 1 g cm-3 for pure water. The temperature of the system was set at 298 K. The Owicki-ScheragaJorgensen preferential sampling technique [ 141 was employed. The equilibrium of the system was established in the preceding at least 3000 K steps, and the averaging was carried out in the following 30-6000 K steps. The interaction between the GLYZ and the solvent molecules was analyzed on the basis of the 2D-rdf which was calculated from the twodimensional normalized histogram. The changes in free energy of solvation were calculated along the rotation of 4 and 19in steps of 5” using SPT with the double-width sampling technique [9]. For the ab initio calculations on the GLYZ, the MIDI-4(d) basis set (see Ref. [15]; the exponents of the polarization functions were taken from Ref. [16]) and the 6-31++G* basis set (see Ref. [17]; the exponents of the polarization and diffuse functions were taken from Ref. [ 181) were used. All MO calculations were carried out using our ABINIT program, and the MC simulation was performed using our SIMPLS program on HP-735 workstations and GAIA-300 personal supercomputers.
3. Results and discussion Relative energies of the four conformations of the GLYZ in the gas phase are shown in Table 3. The (0,O) and (60,O) conformations have comparable Table 1 Glycine zwitterion parameters three-point model of water Site
C(a) WC) C N 0 H(N)
A ((kcal A” mol-I)“‘) I946 1123 387 690 _
in the potential
functions
c ((kcal A” moI-‘)‘n) 79.0 48.2 11.1 21.0
for the
t) -0.3500 0.1885 0.6987 -0.7079 -0.7129 0.4693
115
T. Watanabe et al./Journal of Molecular Structure (Theochem) 397 (1997) 113-119 Table 2 Glycine zwitterion parameters TIP4P model of water
in the potential
A ((kcal Al2 mol-‘) In)
c
C(a) II(C) C N 0
2294
166.7
H(N)
-
Site
1294 391.1 794.7
functions
((kcal A6 mol-‘)‘R)
Table 3 Relative energies (kcal mol-‘) of four conformations zwitterion in the gas phase”
for the
k) -0.4793 0.1822 0.5848 -0.6136 -0.5866 0.4390
39.75 4.111 42.19
Basis
(0.0)
(6090)
(60.90)
(0.90)
MIDL4(d) 6-31+ + G*
0.0 0.0
0.9 1.2
10.2 9.7
15.5 13.8
a Standard geometry
energies; the (60,O) conformation has a slightly higher energy (by about 1 kcal mol-‘) than the (0,O) conformation. The 13= 90” conformations are more unstable than the corresponding 8 = 0” conformations by 9- 16 kcal mol-‘. These results were obtained using standard geometry [3,7,12]. The molecular structure of the GLYZ in the gas phase cannot be determined by geometry optimization, since the neutral form is more stable than the zwitterionic form in the gas phase. The MIDI-4(d) geometry optimization of the (0,O) conformation of the GLYZ caused proton transfer, giving neutral glycine. The instable nature of the GLYZ in the gas phase has been pointed out by several groups 1191.
of the glycine
was used.
The variation in the free energy of hydration along the (0,O) - (60,O) - (60,90) - (0,90) - (0,O) pathway of conformational changes and that along the pathway in the opposite direction are given in Fig. 2. They were calculated by SPT with the three-point water model; the TIP4P model gave very similar profiles for the energy variation. The energies of solvation for the four conformations are listed in Table 4. The relative conformational stabilities in aqueous solution were thus estimated from the sum of the solvation energy and the gas-phase energy. They are listed in Table 5. The 8=0” conformations are more stable than the 13= 90” conformations in aqueous solution, although solvation stabilizes the 0 = 90” conformations more than the 8 = 0” conformations. The most stable conformation in aqueous solution was predicted to be the (60,O) conformation, which is 3 kcal mol-’ 0
-
forward
-----backward
-6-
-6
-8-
-8
-10
I,,,,
(O,O)
(60,O)
(60,90)
(0,90)
(Osu
conformation Fig. 2. Variation in the free energy of salvation with respect to conformational changes. The solid line shows the energy change along the forward path, while the broken line show that in the reverse direction. The statistical perturbation method was applied for every 5” of 6 or 8.
116
T. Watanabe et al./Joumal
of Molecular Structure (Theochem) 397 (1997) 113-119
Table 4 Relative free energies of solvation (kcal mol-‘) for four conformations of the glycine zwitterion in water calculated by Monte Carlo simulation with SPTS
Table 6 Averaged solute-solvent interaction energies (kcal mole’) and number of solvent molecules (in parentheses) in the A, B, and C regions in the 2D-rdf maps”
Model
(0,O)
(6070)
(60.90)
(0,90)
Region
(0.0)
(6070)
(6090)
(0.90)
Three-point model TIP4P model
0.0 0.0
-4.2 -4.0
-7.2 -6.6
-7.8 -7.0
A
-12.14 (0.65) -15.23 (2.26) -14.54 (2.91) -7.22 (7.93) -9.19 (10.84)
-15.98 (1.95) -13.05 (1.24) -14.84 (3.19) -7.52 (7.70) -9.67 (10.89)
-18.01 (1.92) -15.49 (1.16) -17.06 (3.08) -7.45 (7.61) -10.22 (10.69)
-18.75 (1.03) -17.00 (I .98) -17.60 (3.01) -7.57 (8.37) -10.22 (11.38)
B a Average of the forward and backward
values. A+B
lower in energy than the (0,O) conformation. The barrier to rotation about the C-C bond is 6-9 kcal mall’, and the 0 rotation is still largely restricted in aqueous solution. Differences in the free energies of solvation among the four conformations and the hydrogen bonding structures can be well understood from the 2D-rfd maps. Fig. 3 shows contour line maps of the distribution of the solvent molecules. They are drawn with respect to two atomic distances, the distance between the N atom of the GLYZ and the 0 atom of H20, RN_O, and between the C atom of CO2 in the GLYZ and the 0 atom of H20, Rc_o. Three characteristic regions are found on each map. Two of these, A and B, show sharp distributions and are located near C, is R N-O=2.4A, while a broader distribution, found in the region 3 A < Rc_0 < 4 A and 3 A < RN_O < 6 A. The number of solvent molecules in the A, B, and C regions and the averaged solute-solvent interaction energies for one solvent molecule are listed in Table 6. On the 2D-rdf maps, the RN_O distances of the fi and B regions are restricted in the range 2.4 + 0.2 A, indicating that Hz0 molecules in the A and B regions are hydrogen bonding with the NH3 hydrogen atoms. The 0 atom of Hz0 in the A region is 3.0-3.5 A distant from the CO* carbon atom and the A region represents Hz0 in the interior region. On the other Table 5 Relative energies (kcal mall’) for four conformations of the glycine zwitterion in water calculated by Monte Carlo simulation with SPT” Conformation
(030)
(60.0)
(60,90)
(090)
Three-point model TIP4P model
3.3 3.1
0.0 0.0
6.2 6.7
10.9 11.6
‘Relative values with respect to the (60,O) conformation. The values are derived from the solvation energies in Table 4 and MIDW(d) energies in the gas phase which are shown in Table 3.
C A+B+C
a The solvent molecules counted are in the ranges Rc_o < 3.9 A and RN_,, < 3.0 A for the A region, Rco > 3.9 A and RN_O < 3.0 A for the B region, and Rc_0 < 3.9 A and RN_O > 3.0 A for the C region.
hand, the B region locates at 4.2 to about 4.7 A from the CO2 carbon atom and represents Hz0 in the exterior region. The number of solvent molecules hydrogen bonding with the NH3 group is 2.9 to about 3.2. The ratio of the populations in the A and B regions reflects the NH3 conformation; one interior N-H hydrogen atom and two exterior N-H hydrogen atoms are involved in the 9 = 0” conformations, while two interior and one exterior hydrogen atoms are involved in the 4 =60” conformations. Mezei et al. [6] have studied the (0,O) conformation in water, and the hydrogen bonding energy around the NH3 hydrogen atoms was estimated to be -16 kcal mall’, while that around the CO2 region is -12 kcal mol-‘. The value of -14.5 kcal mall’ for the A and B region in the (0,O) conformation corresponds to hydrogen bonding for the NH3 hydrogen atoms. The C region locates around Rc4 = 3.5 A and R N_O from 3 to 6 A. The C region includes Hz0 molecules that are hydrogen bonding with the CO2 group and those associated with the CH2 group. The number of solvent molecules involved in the C region is 7.6 to about 8.4. The hydrogen bonding energies in this region (-7.2/-7.6 to about -8.4 kcal mol-‘) are smaller than the reported value for the CO2 region, -12 kcal mall’ [16], since the Hz0 molecules associated with the CH2 group, whose hydrogen bonding energy is small, are included in the C region. It is interesting that the total solvent-solute interaction
T. Watanabe et al.Nournal of Molecular Structure (Theochem) 397 (1997) 113-119
9 --
117
uw)
a 7 6 -
4
-2 !z4
I
3
-
2
-
[
----
i=I
‘t-
*I
0
I
’ 12
’
1’ 3
I 4
’ 5
’ 6
’ 7
11
’
8
9
I1
01234567
0
12
nltn
3
Rcu /A
a
4
5
Rca
IA
I
I
I,
6
7
a
LuLl 9
-
I
6 -a I
‘9” 24
3 -
---
-
~-----~
2 -
01234567
r+ I
0
1
II
I
I
0123456789
I
L
I
I
: 01234567
Rcu IA Fig. 3. Two-dimensional
L
0
12
3
4 ~c.0
5
6
7
a
9
01234567
IA
radial distribution for the four conformations of the glycine zwitterion: (a), (0,O) conformation: (b) (60.0) conforma-
tion; (c) (0,90) conformation; (d) (6090)
conformation. The contour lines show the normalized histogram in steps of
1.The
outermost line
corresponds to the value of 0. At the top and right-hand sides, one-dimensional radial distribution functions are also shown. The broken lines indicate the first peaks in the one- dimensional radial distribution functions.
118
T. Watanabe et al./Journal of Molecular Structure (Theochem) 397 (1997) 113-l 19
energy in the A + B + C regions, -9.19 x 10.84 = -99.6 kcal mall’, agrees with the solute-solvent interaction energy calculated for the first hydration shell, -99.6 kcal mol-’ [6]. In the (60,O) and (60,90) conformations, the location and shape of the A and B regions are similar to each other, indicating that the orientations of the solvent molecules in the c~5 = 60” conformations are favorable ones for hydrogen bond formation with the NH3 group and are not disturbed by the CO:! conformation. The shapes of the A and B regions on the 2D-rdf maps for +=O” (Fig. 3(a) and (c)) are different from those for 4=60” (Fig. 3b and 3d). The shift of the B region from Rc4 = 4.7 (4 = 60”) to 4.2 (4 = 0”) is due to the conformational change of the NH3 group. This conformational change of NH3 deduces the Rc_o distance of the A region in the 4 = 0” conformations, as is observed in the (0,90) conformation (Fig. 3~). However, the Rc_~ distance of the A region on the (0,O) conformation map is 3.7 A, which is muchOlarger than that on the (0,90) conformation map, 2.9 A. A distortion of the hydrogen bonding around the NH3 group is thus expected. In the (0,O) conformation, the Hz0 molecule, which is hydrogen bonded with the internal N-H atom, is pushed out by the CO2 oxygen atom, and the hydrogen bonding system N-H...0 is bent in the (0,O) conformation. This was also seen in the stereo view of the hydration complex for the (0,O) conformation reported in the literature [6]. This distortion may be caused by intramolecular hydrogen bonding between the NH3 and CO2 groups in the (0,O) conformation and reduces the solute-solvent interaction energy, as shown in Table 6. The solutesolvent interaction energies in the A + B + C region are well correlated with the free energies of solvation in Table 4. The 2D-rdf maps show clearly the differences of hydrogen bonding at the NH3 group among the four conformations. Alagona et al. analyzed the solvation structure of the (0,O) conformation of the GLYZ using differential distance distribution analysis, and found that each external hydrogen atom is associated with one water molecule, while no water molecule is associated with the internal hydrogen atom [7]. The present results, on the other hand, suggest that hydrogen bonding exists between the internal N-H bond and H20, although the orientation of Hz0 is disturbed by the CO2 group and hydrogen bonds are not formed in the
Table 7 Relative energies (kcal mol-‘) of four conformations of the glycine zwitterion in water calculated by the ab initio GB method” Basis
(0.0)
(60,O)
(60.90)
(0,90)
MIDI-4(d)/GB 6-3 I+ffi*/GB 6-31G/GBb
I .4 0.3 0.7
0.0 0.0 0.0
7.5 6.6 6.2
11.2
aEach geometry was optimized including h From Ref. [4].
9.8 9.2
the solvent effect.
most favorable structure. The conformational stability of the GLYZ in aqueous solution predicted in this study supports previous results that were calculated by the ab initio GBI 6-3 1G method including the solvent effect using the generalized Born formula (GB) [4]. Table 7 shows the relative energies calculated by the ab initio GB method with the MIDI-4(d) and 6-31++G* basis sets. They are parallel to those predicted by the Monte Carlo simulations with SPT, although the ab initio GB calculations predicted that the energy difference between the (0,O) and (60,O) conformations is very small. The ab initio MO/GB method is expected to give reliable conformational energies of amino acids in aqueous solution. The present study revealed the conformational stability of the GLYZ in aqueous solution by using a Monte Carlo simulation with statistical perturbation theory, which supports previous ab initio GB calculations. It was shown that 2D-rdf maps are very useful for elucidating the solvent-solute interaction for a system in which two or more functional groups are hydrogen bonding with solvent molecules, and a solvent molecule interacts with more than one functional group of the solute molecule, as in the GLYZ.
Acknowledgements This research was supported by a Grant-in-Aid for Scientific Research (C) no. 08640633 from the Ministry of Education, Science, Sports and Culture, and by a Tsukuba University Research Project A.
References [I] G.
Albrecht and R.B. Corey, J. Am. Chem. Sot., 61 (1939) 1087. R.E. Marsh, Acta Crystallogr., 11 (1958) 654. J. Almlof, A. Kvick and J.O. Thomas, J. Chem. Phys., 59 (1973) 3901.
T. Watanabe et al/Journal of Molecular Structure (Theochem) 397 (1997) 113- I19 [2] J.S. Gaffney, R.C. Pierce and L. Friedman, J. Am. Chem. Sot., 99 (1977) 4293. [3] R. Bonaccorsi, P. Palla and .I. Tomasi, J. Am. Chem. Sot., 106 (1984) 1945. 0. Kikuchi, T. Natsui and T. Kozaki, J. Mol. Struct. Theochem., 207 (1990) 103. B. Deprick-Cote, .I. Lnaglet, J. Caillet, J. Berges, E. Kassab and R. Constanciel, Theoret. Chim. Acta, 82 (1992) 435. [4] 0. Kikuchi, T. Matsuoka, H. Sawahataand 0. Takahashi, J. Mol. Struct. Theochem., 305 (1994) 79.0. Takahashi, H. Sawahata, Y. Ogawa and 0. Kikuchi, J. Mol. Struct. Theochem., in press. [5] S. Roman0 and E. Clementi, Int. J. Quantum Chem., 14 (1978) 839. [6] M. Mezei, P.K. Mehrotra and D.L. Beveridge, J. Biomol. Struct. Dynamics, 2 (1984) 1. [7] G. Alagona, C. Ghio and P.A. Kollman, J. Mol. Struct. Theothem., 166 (1988) 385. G. Alagona and C. Ghio, J. Mol. Liq., 47 (1990) 139. [8] R.W. Zwanzig, J. Chem. Phys., 22 (1954) 1420. M. Mezei and D.L. Beveridge, Ann. N.Y. Acad. Sci., 482 (1986) 1. [9] W.L. Jorgensen and C. Ravinohan, J. Chem. Phys., 83 (1985) 3050. [IO] H. Takase and 0. Kikuchi, Chem. Phys., 181 (1994) 57. [l l] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey and M.L. Klein, J. Chem. Phys., 79 (1983) 926. 1121 L. Carozzo, G. Corongiu, C. Petrongolo and E. Clementi, J. Chem. Phys., 68 (1978) 787.
119
[I31 N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller and E. Teller, I. Chem. Phys., 21 (1953) 1087. M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. [14] J.C. Owicki and H.A. Scheraga, Chem. Phys. Lett., 47 (1977) 600. J.C. Owicki, ACS Symp. Ser. 86, American Chemical Society, Washington DC, 1978, p. 159. W.L. Jorgensen, B. Bigot and J. Chandrasekhar, J. Am. Chem. Sot., 104 (1982) 4584. [15] H. Tatewaki and S. Hujinaga, J. Comput. Chem., l(l980) 205. Y. Sakai, H. Tatewaki and S. Hujinaga, J. Comput. Chem., 2 (1981) 100, 108. [16] S. Hujinaga, Gaussian Basis Sets for Molecular Calculations, Elsevier, Amsterdam, 1984. [17] W.J. Hehre, R. Ditchfield and J.A. Pople, J. Chem. Phys., 56, 2257 (1972). [18] W.J. Hehre, L. Radom, P.v.R. Schleyer and J.A. Pople, Ab Initio Molecular Orbital Theory, John Wiley, New York, 1986. [19] Y.-C. Tse, M.D. Newton, S. Vishveshwara and J.A. Pople, J. Am. Chem. Sot., 100 (1978) 4329. Y. Ding and K. KroghJespersen, Chem. Phys. Len, 199 (1992) 261. F.R. Tortonda, J.L. Pascual-Ahuir, E. Silla and I. Tui?bn, Chem. Phys. L&t., 260 (1996) 21.
Journal of Molecular Structure (Theochem) 397 (1997)113- 119
Monte Carlo simulation study on the conformation and interaction of the glycine zwitterion in aqueous solution Toshio Watanabe, Kazuhiko Hashimoto,
Hideto Takase, Osamu Kikuchi*
Department of Chemistry, University of Tsukuba, Tsukuba 305, Japan
Received 5 June 1996; accepted 24 October 1996
Abstract A Monte Carlo simulation with the statistical perturbation method has been employed to determine the relative stabilities among four representative conformations of the glycine zwitterion (GLYZ) in aqueous solution. The (+,0) = (0,O) and (60,O) conformations, in which the dihedral angle 19between the NCC and CC0 planes is O’, were found to be more stable than the (0,90) and (60,90) conformations, in which 0 =90”. The relative stabilities of these four conformations are well elucidated on the basis of solvent-solute interactions using two-dimensional radial distribution functions. Ab initio MOlGB calculations including the solvent effect via the generalized Born formula gave similar results for the relative stabilities of the four conformations in solution. 0 1997 Elsevier Science B.V. Keywords: Monte Carlo simulation; interactions
Glycine zwitterion;
Conformational
1. Introduction Knowledge of the conformation and interactions of amino acids in water is important for the ultimate understanding of protein hydration and the role of water in biological systems. Glycine is found in the form of a zwitterion in the solid state [I] and in aqueous solution [2]. The conformation of the glycine zwitterion (GLYZ) in aqueous solution has been studied by molecular orbital (MO) calculations in which the solvent effect is included by continuum models [3,4]. Monte Carlo (MC) simulations for the GLYZ in water have been carried out by Roman0 and Clementi [5], Mezei et al. [6], and Alagona et al. [7]. * Corresponding author.
method; Molecular
These simulations concerned the solvation structure for a particular conformation, and revealed a microscopic interaction between the GLYZ and the water molecules. However, little attention has been paid to the conformational energy of the GLYZ in water. Conformational analysis of the GLYZ by MC simulation requires the calculation of the free energy change along the rotation of the C-C bond and/or the C-N bond. There have been no reports of such an analysis. In this study, an MC simulation was applied to four typical conformations of the GLYZ, the (&I!?>= (0, 0), (60,0), (60,90), and (0,90) conformations shown in Fig. 1. Differences in the free energy of solvation were calculated using statistical perturbation theory (WI’) [8,9] and relative stabilities among the four conformations in aqueous solution were predicted.
0166-1280/97/$17.000 1997 Elsevier Science B.V. All rights reserved
PII SO166-1280(96)04942-l
stability; Statistical perturbation
114
T. Watanabe et al./Joumal of Molecular Structure (Theochem) 397 (1997) 113-l 19
Fig. 1. Definition of the two rotation angles, 4 and 8, for the zwitterionic form of glycine: 4 indicates the rotation of the NH3 fragment, while 19 indicates the rotation of the CO* fragment. The conformation shown is (6, 0) = (0,O).
The analysis of the interaction between the GLYZ and solvent molecules introduces special difficulties, since the GLYZ has two hydrogen bonding groups, NHf and CO;, and the calculated distribution functions contain admixtures from various distinctly different modes of solvation [6]. In order to avoid this difficulty, we employed here two-dimensional distribution functions (2D-rdfs) which are represented by two distances between specific atoms. The relative stability of the four conformers in aqueous solution is well elucidated on the basis of the solute-solvent interaction, which is clearly displayed on the 2D-rdf maps.
2. Calculation Two potential functions were used for the water molecules in the MC simulation. One is the threepoint model, which has been determined by ab initio MIDI-4(d) calculations in our previous paper [lo], and the other is the TIP4P model [ll]. The intermolecular potential functions between the GLYZ and Hz0 were newly determined. In this procedure, four conformations were used for the GLYZ, fixing the bond lengths and angles at the standard values [3,7,12]. The MIDI-4(d) calculations were carried out for dimer orientations which were selected in the vicinity of several energy-minimum geometries and for those selected randomly. About 500 configurations were selected for the GLYZ-Hz0 dimer. On the basis of the calculated interaction energies for the dimers, the potential function between the GLYZ and H20, which is constructed by Lennard-Jones (12-6) terms and the Coulomb term,
was determined. The parameters determined are listed in Table 1 and Table 2. The MC simulations were performed in the NVT ensemble according to the standard Metropolis method [13] for one GLYZ solute molecule and 212 water molecules in the cubic cell. The volume of the cell was determined using a density of 1 g cm-3 for pure water. The temperature of the system was set at 298 K. The Owicki-ScheragaJorgensen preferential sampling technique [ 141 was employed. The equilibrium of the system was established in the preceding at least 3000 K steps, and the averaging was carried out in the following 30-6000 K steps. The interaction between the GLYZ and the solvent molecules was analyzed on the basis of the 2D-rdf which was calculated from the twodimensional normalized histogram. The changes in free energy of solvation were calculated along the rotation of 4 and 19in steps of 5” using SPT with the double-width sampling technique [9]. For the ab initio calculations on the GLYZ, the MIDI-4(d) basis set (see Ref. [15]; the exponents of the polarization functions were taken from Ref. [16]) and the 6-31++G* basis set (see Ref. [17]; the exponents of the polarization and diffuse functions were taken from Ref. [ 181) were used. All MO calculations were carried out using our ABINIT program, and the MC simulation was performed using our SIMPLS program on HP-735 workstations and GAIA-300 personal supercomputers.
3. Results and discussion Relative energies of the four conformations of the GLYZ in the gas phase are shown in Table 3. The (0,O) and (60,O) conformations have comparable Table 1 Glycine zwitterion parameters three-point model of water Site
C(a) WC) C N 0 H(N)
A ((kcal A” mol-I)“‘) I946 1123 387 690 _
in the potential
functions
c ((kcal A” moI-‘)‘n) 79.0 48.2 11.1 21.0
for the
t) -0.3500 0.1885 0.6987 -0.7079 -0.7129 0.4693
115
T. Watanabe et al./Journal of Molecular Structure (Theochem) 397 (1997) 113-119 Table 2 Glycine zwitterion parameters TIP4P model of water
in the potential
A ((kcal Al2 mol-‘) In)
c
C(a) II(C) C N 0
2294
166.7
H(N)
-
Site
1294 391.1 794.7
functions
((kcal A6 mol-‘)‘R)
Table 3 Relative energies (kcal mol-‘) of four conformations zwitterion in the gas phase”
for the
k) -0.4793 0.1822 0.5848 -0.6136 -0.5866 0.4390
39.75 4.111 42.19
Basis
(0.0)
(6090)
(60.90)
(0.90)
MIDL4(d) 6-31+ + G*
0.0 0.0
0.9 1.2
10.2 9.7
15.5 13.8
a Standard geometry
energies; the (60,O) conformation has a slightly higher energy (by about 1 kcal mol-‘) than the (0,O) conformation. The 13= 90” conformations are more unstable than the corresponding 8 = 0” conformations by 9- 16 kcal mol-‘. These results were obtained using standard geometry [3,7,12]. The molecular structure of the GLYZ in the gas phase cannot be determined by geometry optimization, since the neutral form is more stable than the zwitterionic form in the gas phase. The MIDI-4(d) geometry optimization of the (0,O) conformation of the GLYZ caused proton transfer, giving neutral glycine. The instable nature of the GLYZ in the gas phase has been pointed out by several groups 1191.
of the glycine
was used.
The variation in the free energy of hydration along the (0,O) - (60,O) - (60,90) - (0,90) - (0,O) pathway of conformational changes and that along the pathway in the opposite direction are given in Fig. 2. They were calculated by SPT with the three-point water model; the TIP4P model gave very similar profiles for the energy variation. The energies of solvation for the four conformations are listed in Table 4. The relative conformational stabilities in aqueous solution were thus estimated from the sum of the solvation energy and the gas-phase energy. They are listed in Table 5. The 8=0” conformations are more stable than the 13= 90” conformations in aqueous solution, although solvation stabilizes the 0 = 90” conformations more than the 8 = 0” conformations. The most stable conformation in aqueous solution was predicted to be the (60,O) conformation, which is 3 kcal mol-’ 0
-
forward
-----backward
-6-
-6
-8-
-8
-10
I,,,,
(O,O)
(60,O)
(60,90)
(0,90)
(Osu
conformation Fig. 2. Variation in the free energy of salvation with respect to conformational changes. The solid line shows the energy change along the forward path, while the broken line show that in the reverse direction. The statistical perturbation method was applied for every 5” of 6 or 8.
116
T. Watanabe et al./Joumal
of Molecular Structure (Theochem) 397 (1997) 113-119
Table 4 Relative free energies of solvation (kcal mol-‘) for four conformations of the glycine zwitterion in water calculated by Monte Carlo simulation with SPTS
Table 6 Averaged solute-solvent interaction energies (kcal mole’) and number of solvent molecules (in parentheses) in the A, B, and C regions in the 2D-rdf maps”
Model
(0,O)
(6070)
(60.90)
(0,90)
Region
(0.0)
(6070)
(6090)
(0.90)
Three-point model TIP4P model
0.0 0.0
-4.2 -4.0
-7.2 -6.6
-7.8 -7.0
A
-12.14 (0.65) -15.23 (2.26) -14.54 (2.91) -7.22 (7.93) -9.19 (10.84)
-15.98 (1.95) -13.05 (1.24) -14.84 (3.19) -7.52 (7.70) -9.67 (10.89)
-18.01 (1.92) -15.49 (1.16) -17.06 (3.08) -7.45 (7.61) -10.22 (10.69)
-18.75 (1.03) -17.00 (I .98) -17.60 (3.01) -7.57 (8.37) -10.22 (11.38)
B a Average of the forward and backward
values. A+B
lower in energy than the (0,O) conformation. The barrier to rotation about the C-C bond is 6-9 kcal mall’, and the 0 rotation is still largely restricted in aqueous solution. Differences in the free energies of solvation among the four conformations and the hydrogen bonding structures can be well understood from the 2D-rfd maps. Fig. 3 shows contour line maps of the distribution of the solvent molecules. They are drawn with respect to two atomic distances, the distance between the N atom of the GLYZ and the 0 atom of H20, RN_O, and between the C atom of CO2 in the GLYZ and the 0 atom of H20, Rc_o. Three characteristic regions are found on each map. Two of these, A and B, show sharp distributions and are located near C, is R N-O=2.4A, while a broader distribution, found in the region 3 A < Rc_0 < 4 A and 3 A < RN_O < 6 A. The number of solvent molecules in the A, B, and C regions and the averaged solute-solvent interaction energies for one solvent molecule are listed in Table 6. On the 2D-rdf maps, the RN_O distances of the fi and B regions are restricted in the range 2.4 + 0.2 A, indicating that Hz0 molecules in the A and B regions are hydrogen bonding with the NH3 hydrogen atoms. The 0 atom of Hz0 in the A region is 3.0-3.5 A distant from the CO* carbon atom and the A region represents Hz0 in the interior region. On the other Table 5 Relative energies (kcal mall’) for four conformations of the glycine zwitterion in water calculated by Monte Carlo simulation with SPT” Conformation
(030)
(60.0)
(60,90)
(090)
Three-point model TIP4P model
3.3 3.1
0.0 0.0
6.2 6.7
10.9 11.6
‘Relative values with respect to the (60,O) conformation. The values are derived from the solvation energies in Table 4 and MIDW(d) energies in the gas phase which are shown in Table 3.
C A+B+C
a The solvent molecules counted are in the ranges Rc_o < 3.9 A and RN_,, < 3.0 A for the A region, Rco > 3.9 A and RN_O < 3.0 A for the B region, and Rc_0 < 3.9 A and RN_O > 3.0 A for the C region.
hand, the B region locates at 4.2 to about 4.7 A from the CO2 carbon atom and represents Hz0 in the exterior region. The number of solvent molecules hydrogen bonding with the NH3 group is 2.9 to about 3.2. The ratio of the populations in the A and B regions reflects the NH3 conformation; one interior N-H hydrogen atom and two exterior N-H hydrogen atoms are involved in the 9 = 0” conformations, while two interior and one exterior hydrogen atoms are involved in the 4 =60” conformations. Mezei et al. [6] have studied the (0,O) conformation in water, and the hydrogen bonding energy around the NH3 hydrogen atoms was estimated to be -16 kcal mall’, while that around the CO2 region is -12 kcal mol-‘. The value of -14.5 kcal mall’ for the A and B region in the (0,O) conformation corresponds to hydrogen bonding for the NH3 hydrogen atoms. The C region locates around Rc4 = 3.5 A and R N_O from 3 to 6 A. The C region includes Hz0 molecules that are hydrogen bonding with the CO2 group and those associated with the CH2 group. The number of solvent molecules involved in the C region is 7.6 to about 8.4. The hydrogen bonding energies in this region (-7.2/-7.6 to about -8.4 kcal mol-‘) are smaller than the reported value for the CO2 region, -12 kcal mall’ [16], since the Hz0 molecules associated with the CH2 group, whose hydrogen bonding energy is small, are included in the C region. It is interesting that the total solvent-solute interaction
T. Watanabe et al.Nournal of Molecular Structure (Theochem) 397 (1997) 113-119
9 --
117
uw)
a 7 6 -
4
-2 !z4
I
3
-
2
-
[
----
i=I
‘t-
*I
0
I
’ 12
’
1’ 3
I 4
’ 5
’ 6
’ 7
11
’
8
9
I1
01234567
0
12
nltn
3
Rcu /A
a
4
5
Rca
IA
I
I
I,
6
7
a
LuLl 9
-
I
6 -a I
‘9” 24
3 -
---
-
~-----~
2 -
01234567
r+ I
0
1
II
I
I
0123456789
I
L
I
I
: 01234567
Rcu IA Fig. 3. Two-dimensional
L
0
12
3
4 ~c.0
5
6
7
a
9
01234567
IA
radial distribution for the four conformations of the glycine zwitterion: (a), (0,O) conformation: (b) (60.0) conforma-
tion; (c) (0,90) conformation; (d) (6090)
conformation. The contour lines show the normalized histogram in steps of
1.The
outermost line
corresponds to the value of 0. At the top and right-hand sides, one-dimensional radial distribution functions are also shown. The broken lines indicate the first peaks in the one- dimensional radial distribution functions.
118
T. Watanabe et al./Journal of Molecular Structure (Theochem) 397 (1997) 113-l 19
energy in the A + B + C regions, -9.19 x 10.84 = -99.6 kcal mall’, agrees with the solute-solvent interaction energy calculated for the first hydration shell, -99.6 kcal mol-’ [6]. In the (60,O) and (60,90) conformations, the location and shape of the A and B regions are similar to each other, indicating that the orientations of the solvent molecules in the c~5 = 60” conformations are favorable ones for hydrogen bond formation with the NH3 group and are not disturbed by the CO:! conformation. The shapes of the A and B regions on the 2D-rdf maps for +=O” (Fig. 3(a) and (c)) are different from those for 4=60” (Fig. 3b and 3d). The shift of the B region from Rc4 = 4.7 (4 = 60”) to 4.2 (4 = 0”) is due to the conformational change of the NH3 group. This conformational change of NH3 deduces the Rc_o distance of the A region in the 4 = 0” conformations, as is observed in the (0,90) conformation (Fig. 3~). However, the Rc_~ distance of the A region on the (0,O) conformation map is 3.7 A, which is muchOlarger than that on the (0,90) conformation map, 2.9 A. A distortion of the hydrogen bonding around the NH3 group is thus expected. In the (0,O) conformation, the Hz0 molecule, which is hydrogen bonded with the internal N-H atom, is pushed out by the CO2 oxygen atom, and the hydrogen bonding system N-H...0 is bent in the (0,O) conformation. This was also seen in the stereo view of the hydration complex for the (0,O) conformation reported in the literature [6]. This distortion may be caused by intramolecular hydrogen bonding between the NH3 and CO2 groups in the (0,O) conformation and reduces the solute-solvent interaction energy, as shown in Table 6. The solutesolvent interaction energies in the A + B + C region are well correlated with the free energies of solvation in Table 4. The 2D-rdf maps show clearly the differences of hydrogen bonding at the NH3 group among the four conformations. Alagona et al. analyzed the solvation structure of the (0,O) conformation of the GLYZ using differential distance distribution analysis, and found that each external hydrogen atom is associated with one water molecule, while no water molecule is associated with the internal hydrogen atom [7]. The present results, on the other hand, suggest that hydrogen bonding exists between the internal N-H bond and H20, although the orientation of Hz0 is disturbed by the CO2 group and hydrogen bonds are not formed in the
Table 7 Relative energies (kcal mol-‘) of four conformations of the glycine zwitterion in water calculated by the ab initio GB method” Basis
(0.0)
(60,O)
(60.90)
(0,90)
MIDI-4(d)/GB 6-3 I+ffi*/GB 6-31G/GBb
I .4 0.3 0.7
0.0 0.0 0.0
7.5 6.6 6.2
11.2
aEach geometry was optimized including h From Ref. [4].
9.8 9.2
the solvent effect.
most favorable structure. The conformational stability of the GLYZ in aqueous solution predicted in this study supports previous results that were calculated by the ab initio GBI 6-3 1G method including the solvent effect using the generalized Born formula (GB) [4]. Table 7 shows the relative energies calculated by the ab initio GB method with the MIDI-4(d) and 6-31++G* basis sets. They are parallel to those predicted by the Monte Carlo simulations with SPT, although the ab initio GB calculations predicted that the energy difference between the (0,O) and (60,O) conformations is very small. The ab initio MO/GB method is expected to give reliable conformational energies of amino acids in aqueous solution. The present study revealed the conformational stability of the GLYZ in aqueous solution by using a Monte Carlo simulation with statistical perturbation theory, which supports previous ab initio GB calculations. It was shown that 2D-rdf maps are very useful for elucidating the solvent-solute interaction for a system in which two or more functional groups are hydrogen bonding with solvent molecules, and a solvent molecule interacts with more than one functional group of the solute molecule, as in the GLYZ.
Acknowledgements This research was supported by a Grant-in-Aid for Scientific Research (C) no. 08640633 from the Ministry of Education, Science, Sports and Culture, and by a Tsukuba University Research Project A.
References [I] G.
Albrecht and R.B. Corey, J. Am. Chem. Sot., 61 (1939) 1087. R.E. Marsh, Acta Crystallogr., 11 (1958) 654. J. Almlof, A. Kvick and J.O. Thomas, J. Chem. Phys., 59 (1973) 3901.
T. Watanabe et al/Journal of Molecular Structure (Theochem) 397 (1997) 113- I19 [2] J.S. Gaffney, R.C. Pierce and L. Friedman, J. Am. Chem. Sot., 99 (1977) 4293. [3] R. Bonaccorsi, P. Palla and .I. Tomasi, J. Am. Chem. Sot., 106 (1984) 1945. 0. Kikuchi, T. Natsui and T. Kozaki, J. Mol. Struct. Theochem., 207 (1990) 103. B. Deprick-Cote, .I. Lnaglet, J. Caillet, J. Berges, E. Kassab and R. Constanciel, Theoret. Chim. Acta, 82 (1992) 435. [4] 0. Kikuchi, T. Matsuoka, H. Sawahataand 0. Takahashi, J. Mol. Struct. Theochem., 305 (1994) 79.0. Takahashi, H. Sawahata, Y. Ogawa and 0. Kikuchi, J. Mol. Struct. Theochem., in press. [5] S. Roman0 and E. Clementi, Int. J. Quantum Chem., 14 (1978) 839. [6] M. Mezei, P.K. Mehrotra and D.L. Beveridge, J. Biomol. Struct. Dynamics, 2 (1984) 1. [7] G. Alagona, C. Ghio and P.A. Kollman, J. Mol. Struct. Theothem., 166 (1988) 385. G. Alagona and C. Ghio, J. Mol. Liq., 47 (1990) 139. [8] R.W. Zwanzig, J. Chem. Phys., 22 (1954) 1420. M. Mezei and D.L. Beveridge, Ann. N.Y. Acad. Sci., 482 (1986) 1. [9] W.L. Jorgensen and C. Ravinohan, J. Chem. Phys., 83 (1985) 3050. [IO] H. Takase and 0. Kikuchi, Chem. Phys., 181 (1994) 57. [l l] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey and M.L. Klein, J. Chem. Phys., 79 (1983) 926. 1121 L. Carozzo, G. Corongiu, C. Petrongolo and E. Clementi, J. Chem. Phys., 68 (1978) 787.
119
[I31 N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller and E. Teller, I. Chem. Phys., 21 (1953) 1087. M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. [14] J.C. Owicki and H.A. Scheraga, Chem. Phys. Lett., 47 (1977) 600. J.C. Owicki, ACS Symp. Ser. 86, American Chemical Society, Washington DC, 1978, p. 159. W.L. Jorgensen, B. Bigot and J. Chandrasekhar, J. Am. Chem. Sot., 104 (1982) 4584. [15] H. Tatewaki and S. Hujinaga, J. Comput. Chem., l(l980) 205. Y. Sakai, H. Tatewaki and S. Hujinaga, J. Comput. Chem., 2 (1981) 100, 108. [16] S. Hujinaga, Gaussian Basis Sets for Molecular Calculations, Elsevier, Amsterdam, 1984. [17] W.J. Hehre, R. Ditchfield and J.A. Pople, J. Chem. Phys., 56, 2257 (1972). [18] W.J. Hehre, L. Radom, P.v.R. Schleyer and J.A. Pople, Ab Initio Molecular Orbital Theory, John Wiley, New York, 1986. [19] Y.-C. Tse, M.D. Newton, S. Vishveshwara and J.A. Pople, J. Am. Chem. Sot., 100 (1978) 4329. Y. Ding and K. KroghJespersen, Chem. Phys. Len, 199 (1992) 261. F.R. Tortonda, J.L. Pascual-Ahuir, E. Silla and I. Tui?bn, Chem. Phys. L&t., 260 (1996) 21.