17670
J. Phys. Chem. 1996, 100, 17670-17677
Determination of Potential Parameters for Amino Acid Zwitterions Oh Young Kwon, Su Yeon Kim, and Kyoung Tai No*,† Department of Chemistry, Soong Sil UniVersity, Sang Do 5 Dong 1-1, Dong Jak Gu, Seoul 156-743, Korea
Young Kee Kang Department of Chemistry, Chungbuk National UniVersity, Cheongju, Chungbuk 361-763, Korea
Mu Shik Jhon† Department of Chemistry, Korea AdVanced Institute of Science and Technology, 373-1 Kusung-dong Yusung-gu, Taejon 305-701, Korea
Harold A. Scheraga* Baker Laboratory of Chemistry, Cornell UniVersity, Ithaca, New York 14853-1301 ReceiVed: April 24, 1996; In Final Form: August 16, 1996X
An empirical potential energy function for description of the interactions between amino acid zwitterions has been developed. Potential energy surfaces of molecular ion pairs, methylammonium ion-acetate ion, glycine zwitterion-methyl ammonium ion, and glycine zwitterion-acetate ion, were obtained with 6-31G** ab initio molecular orbital calculations. The point charges of the ions were obtained with the modified partial equalization of orbital electronegativity (M-PEOE) method, the coefficients of the attractive part of the nonbonded potential were calculated from charge-dependent effective atomic polarizabilities, and the coefficients of the repulsive part of the nonbonded potential were obtained from the equilibrium conditions of molecules in molecular crystals and by reproducing the lattice energies of amino acid molecular crystals. The hydrogenbond model proposed by No et al. was introduced, and the potential parameters of the model were determined in order that the potential energy function could reproduce both the ab initio potential energy surfaces and the experimental structures of some amino acid molecular crystals. The minimum-energy crystal structures obtained through crystal packing computations with the empirical potential energy function agreed well with X-ray crystal structures.
I. Introduction Amino acids in aqueous solution are ionized and can act as acids or bases. Several charged species exist at equilibrium depending on the pH of the solution. R-Amino acids having one amino group and one carboxylic acid group crystallize from neutral aqueous solutions as fully ionized species known as zwitterions or dipolar ions, and each amino acid has both a positive and a negative charge. Since the N-H+‚‚‚-O-C hydrogen bonds are very strong in amino acid crystals, they stabilize the zwitterion form and play a major role in determining the structures of amino acid crystals. Therefore, this type of hydrogen bond is characteristic of amino acid crystals. Investigation of the interactions -COO-‚‚‚+H3N-, -COO-‚‚‚H2O, and -NH3+‚‚‚OH2 is very important for illumination of the behavior of peptides in solution. Therefore, reliable potential energy functions for these hydrogen-bonded pairs are necessary for theoretical investigations of the behavior of peptides in solution or in crystals. Such zwitterions, for example, γ-aminobutyric acid (GABA), act as essential neurotransmitters in the central nervous system.1,2 Therefore, investigation of the behavior of zwitterions in aqueous solution and of the interaction of zwitterions with their receptors is very important. The ammonium ion is an essential functional group in almost all neurotransmitters and in many † Member of the Center for Molecular Science, Korea. * Authors to whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, October 15, 1996.
S0022-3654(96)01180-X CCC: $12.00
other important biological molecules.3 Therefore, the precise expression of the interaction between -NH3+ and -CO2- is very important for the investigation of receptor-ligand interactions that occur in biological systems.4 Pair potential energy functions describing the interactions between water and biomolecules have been investigated by several workers.5-14 Analytical potential energy functions have been determined for both water-neutral amino acid5a-c,7d,8a,c,10b and water-amino acid zwitterion5d-g,6,7d-g,8b,c,10a,c,11-14 pairs. There are numerous studies of the interactions between -NH3+ and water and between -CO2- and water by molecular orbital calculations4,5a,d-g,6,7d-g,8b,9,11,12 and by computer simulations.5c,8b,10,13,14 Analytical potential energy functions for ion-molecule pairs have also been proposed.5d-g,6,7d-f,8b However, only a few molecular orbital calculations4,6,11b have been carried out for the interaction between -CO2- and -NH3+, and the potential energy surface for this interaction6,7b is not yet described with a proper analytical function. Experimentally, information about the structures and binding energies between molecules that contain these ionic groups can be obtained from X-ray structures and heats of sublimation of crystals of amino acids. Since the hydrogen bond in a zwitterion pair is very strong, the information inherent in crystal structures is confined to a limited number of configurations, i.e., to only energetically stable configurations. For this reason, it is necessary to carry out molecular orbital calculations to obtain information about the whole potential energy surface for the interaction between ionic pairs. © 1996 American Chemical Society
Potential Parameters for Amino Acid Zwitterions
J. Phys. Chem., Vol. 100, No. 44, 1996 17671
Since the contribution of the electrostatic interaction between the two dipoles, µN-H and µO-C, to the N-H+‚‚‚-O-C interaction energy is dominant, a physically realistic representation of the electrostatic interaction is very important for describing the angular dependence and the energy of the hydrogen bond. Momany et al.7c calculated the minimumenergy packing structures and the lattice energies of crystals of various amino acids with their empirical potential energy function, ECEPP. In the calculations, all potential parameters of the atoms in amino acid zwitterions, except the net atomic charges, were taken from the potential energy function that had been developed for neutral peptides.7b Although they concluded that “the potential energy functions are applicable to conformational studies of polypeptides and proteins”, the resulting binding energies that correspond to the lattice energies were too small. This is due mainly to low estimates of the electrostatic interactions between amino acids. Voogd et al.15 and No et al.16 calculated the contribution of the electrostatic energy to the lattice energy of crystals of amino acids. Such contributions, estimated from heat of sublimation data,17 were predominant (over 60%) compared with the other interaction terms. The purpose of this work is to develop a set of potential energy functions for the -NH3+‚‚‚-O2C- pair. For simplicity, the point charges are placed on the centers of atoms, and the point charges are calculated by the modified partial equalization of orbital electronegativity (M-PEOE) method.18 The hydrogenbond model, which was proposed by No et al.,19 is utilized because it reproduces the angular dependence of hydrogen bonds formed between amides and between carboxylic acids. Moreover, the merit of the hydrogen-bond model is that it does not contain an angle variable in its functional form. The nonbonded potential parameters are determined by using the equilibrium conditions of molecular crystals. Finally, to test the reliability of the newly developed set of potential energy functions, crystal packing studies of some amino acid zwitterions are carried out and compared with experimental data. II. Computations a. Model Compounds and Calculation of Potential Energy Surface. The methylammonium cation, acetate anion, and glycine zwitterion were used as rigid-geometry model compounds (Figure 1). The potential energy surfaces (PES) for the intermolecular interactions of the ion pairs of the model compounds, methylammonium ion-acetate ion (Figure 2a), glycine zwitterion-methyl ammonium ion (Figure 2b), and glycine zwitterion-acetate ion (Figure 2c), were calculated with the 6-31G** ab initio MO method at the Hartree-Fock level. The minimum-energy conformations of the model compounds and their ion pairs were also optimized with 6-31G** ab initio MO calculations. In the calculation of the PES of the molecular ion pairs, the geometry and conformation of each monomer were fixed at their minimum-energy conformations. All the ab initio MO calculations were carried out with the Gaussian92 program.20 To describe the relative orientation of the two molecules that participate in the hydrogen bond in the binary complex, six degrees of freedom are necessary if the two molecules are assumed rigid. As shown in Figure 2a-c, the geometrical parameters are the following: the N5‚‚‚O5 distance (rNO), the Cam 1 -N5‚‚‚O5 bond angle (θCNO), the H5‚‚‚O5-C5 bond angle (θHOC), the Cam 1 -H1‚‚‚N5‚‚‚O5 dihedral angle (φCHNO), the H5-N5‚‚‚O5‚‚‚C5 dihedral angle (φHNOC), and the N5‚‚‚O5-C5ac am Cac 1 dihedral angle (φNOCC). The C1 atoms, C1 and C1 , in Figure 2a-c are bonded to the ammonium ion and acetate ion,
Figure 1. Optimized geometry, bond lengths (in angstroms), and bond angles (in degrees), the M-PEOE net atomic charges, δ’s (in esu), and the dispersion coefficients, Cii’s (in kcal Å6/mol), of the model compounds treated here, namely, (a, top) methyl ammonium ion, (b, middle) acetate ion, and (c, bottom) glycine zwitterion.
respectively. Among these six geometrical parameters, three are varied for the calculation of potential energy surface of the methylammonium ion-acetate ion and glycine zwitterionacetate ion complexes. Four geometrical parameters are varied for the glycine zwitterion-methylammonium ion complex. The other geometrical parameters are fixed at the values of the equilibrium geometries (Figure 1). For each complex, about 120 configurations were sampled, and the ab initio stabilization energy (Vst), which corresponds to the dimerization energy at 0 K, was calculated at each sampled point. The values of Vst obtained from the ab initio calculations correspond to the difference between the energy of the ionic complex and the sum of two monomer energies, Vst ) Vcomplex - (Vmonomer1 + Vmonomer2). In Table 1, the sampling regions and intervals of the variables are summarized for each pair. b. Empirical Potential Energy Function. The stabilization
17672 J. Phys. Chem., Vol. 100, No. 44, 1996
Kwon et al.
TABLE 1: Sampled Configurations of the Model Binary Complexes for Calculating the ab Initio Potential Energy Surfaces rNO, ∆ra (Å)
θCNO,b ∆θa (deg)
methylammonium-acetate ion
1.5-4.5, 0.2
0-170, 10
glycine zwitterion-methylammonium ion
1.5-4.5, 0.2
0-170, 10
glycine zwitterion-acetate ion
1.5-4.5, 0.2
0-170, 10
model complex
φCHNO,c ∆φa (deg)
φHNOC,d ∆φa(deg) -87 to -177, 30 -102 to -162, 30 69 to 159, 30 84 to 174, 30 -48 to -78, 30 -138 to -153, 15
e 99-159, 30 114-144, 30 e
a ∆ represents the grid intervals for r, θ, and φ, in the PES calculations. b θCNO is the bond angle formed by the C1am, N5, and O5 atoms. c φCHNO is the dihedral angle formed by the C1am, H1, N5, and O5 atoms. d φHNOC is the dihedral angle formed by the H5, N5, O5, and C5 atoms. e φCHNO was not changed; it was fixed at the value of the equilibrium conformation (see Figure 2).
TABLE 2: Classification of Atom Types and Their Descriptions for Zwitterions atom typea (this work) H1 H5 C1 C5 N5 O5
description
designation for M-PEOE methodb
designation for CDEAP methodc
aliphatic bonded to N+ aliphatic in carboxylate ion (CO2-) in ammonium ion (sp3) in carboxylate ion (CO2-)
H H(N+) C(sp3) C(O-) N+ O-
H1 H1 C4 C3 N4 O1
a
The atom types H1 and C1 were defined in our previous paper, ref 22. b The atomic species classified for the M-PEOE method are summarized in Table 3 of the Appendix of ref 18b. c Reference 21.
Figure 2. Minimum-energy conformations of model binary complexes: (a, top) CH3-NH3+‚‚‚-O2C-CH3, (b, middle) CH3NH3+‚‚‚glycine zwitterion, (c, bottom) CH3-CO2-‚‚‚glycine zwitterion.
energy due to intermolecular interaction was described by the following empirical effective pair potential energy function.
Vst ) ∑ i>j
qiqj rij
+ VHB + ∑ i>j
(
Aij
rij12
-
Cij
)
rij6
(1)
where qi and qj are the charges on atoms i and j, rij is their separation distance, VHB is a hydrogen-bond potential defined below, and the coefficients Aij and Cij are defined below. The atom types, classified according to the chemical environments
of the atoms in the model compounds, are listed in Table 2. In the third and fourth columns are given the atom types used in the M-PEOE18 and the charge-dependent effective atomic polarizability (CDEAP) calculation method,21 respectively. Since the atom types of the CDEAP method are classified mainly according to the valence states of the atoms in molecules and those of the M-PEOE method are classified according to the chemical environments of the atoms in molecules, the two methods use different classifications of the atom types. It was assumed that the net atomic charges are located on the atomic centers, and they were calculated with the M-PEOE18 method because there are no experimental gas phase electrical moments for the model compounds introduced in this work. M-PEOE is an empirical point charge calculation method in which the empirical parameters are adjusted in order that the calculated point charges can reproduce experimental electrical multipole moments and ab initio electrostatic potentials. The monomer charges (in both the monomeric units and the complex) were fixed throughout the calculations. The relative contribution of the electrostatic interaction energy between amino acid zwitterions to the lattice energy of amino acid molecular crystal was investigated by No et al.16a with several point charge sets including the M-PEOE charges. The difference between the electrostatic interaction energy and the lattice energy at 0 K corresponds to the sum of the hydrogen-bond, nonbonded, polarization, etc., energies. It is very important to balance (by repeated optimization trials) the electrostatic interaction energy component with the other energy components, for a correct prediction of the physical properties of molecules, i.e., conformation, interaction energy, etc. From this point of view, it was concluded that the point charges taken as the M-PEOE charges are physically realistic. The point charges are summarized in Table 3. Lennard-Jones 6-12 type functions were used for the nonbonded interactions. The attractive coefficients, Cii’s, were calculated with the Slater-Kirkwood formula23
[ ]
Ri 4 ep Cii ) 3 m 1/2 (R /N )1/2 e i i 2
(2)
Potential Parameters for Amino Acid Zwitterions
J. Phys. Chem., Vol. 100, No. 44, 1996 17673
TABLE 3: Net Atomic Point Charges Calculated with the M-PEOE Method and the Attractive Nonbonded Potential Parameters, Cii, Calculated with the CDEAP Method and the Slater-Kirkwood Formula molecule
atom typea
point charge (esu)
Cii (kcal Å6/mol)
H1 C1 H5 N5 H1 C1 C5 O5 H1 C1 H5 N5 C5 O5
0.1161 0.0854 0.3410 -0.4566 -0.0055 -0.0984 0.6808 -0.7816 0.0123 0.1945 0.3368 -0.5314 0.8685 -0.7845
40.68 427.43 30.67 577.16 40.68 427.43 223.01 450.12 40.68 427.43 30.67 577.16 223.01 450.12
methyl ammonium ion
acetate ion
glycine zwitterion
a
Figure 3. Coordinate system, rNH (in angstroms), θNHO (in degrees), θHOC (in degrees), and φNHOC (in degrees), which is usually used to describe a hydrogen-bonded system.
From column 1 of Table 2.
where e, p, and me have their usual meanings. Ni and Ri are the number of effective electrons and effective atomic polarizability, respectively.
VNB ) ∑ i>j
( ) Aij
rij12
-
Cij rij6
(3)
Standard combination rules were used.
ij ) 0.25Cij2/Aij ij ) (iijj)1/2
σij ) (Aij/Cij)1/6
(4)
σij ) (σii + σjj)/2
(5)
The effective atomic polarizabilities, Ri’s, were calculated with the CDEAP method.21 The polarizability of the ith atomic species in a molecule, Ri*, is described as a function of the net atomic charges,
R*i ) R* i,0 - aidqi
(6)
where R* i,o, ai, and dqi are the inherent effective atomic polarizability, charge coefficient, and net atomic charge of the ith atomic species in a molecule. The Cii’s for the atomic species appearing in this work are listed in the last column of Table 3. The repulsive coefficients, Aii’s, were adjusted to satisfy the equilibrium conditions of the molecules in molecular crystals and to reproduce the lattice energies of the amino acid molecular crystals. In the calculations, the molecules in the crystals were assumed to be rigid bodies. The details of the constraint method and the fitting procedures are described in our previous paper.22 Since zwitterions in amino acid crystals change to a neutral form upon sublimation through intramolecular proton transfer, the lattice energy cannot be estimated17 simply from the experimental enthalpy of sublimation. For the estimation of the lattice energies of some amino acid crystals, Voogd et al.15 and No et al.16 carried out extensive ab initio calculations of proton transfer energies of several amino acids. The lattice energies of several amino acids were taken from the work of No et al.16a and were used as constraints for determining the repulsive coefficients, Aii’s. For the description of the hydrogen bonds formed between -CO2- and -NH3+, the hydrogen-bond model proposed by No et al.19 was introduced. The hydrogen-bond model describes the ab initio PES of hydrogen-bonded molecular pairs, amide-
Figure 4. Coordinate system, rHO, rNO, rCH, and rNC, which are introduced in this work for describing the hydrogen bond. Distances are in angstroms.
amide, amide-carboxylic acid, and carboxylic acid-carboxylic acid very well. With that model, the angular dependence of the hydrogen-bond energy could be described without introducing an explicit angular-dependent potential energy term. Figure 3 shows the geometrical parameters that are generally adopted for describing the hydrogen-bond system, N-H+‚‚‚-O-C. If rNH and rOC are fixed, then four degrees of freedom, Viz., rHO, θNHO, θHOC, and φNHOC, are necessary for describing the relative orientations in the hydrogen-bonded system. Then, the hydrogenbond potential energy can be described as a function of the four geometrical parameters.
V′HB ) f (rHO,θNHO,θHOC,φNHOC)
(7)
If the coordinate system described in Figure 4 is introduced instead of the coordinate system of Figure 3, the hydrogenbond potential can also be described with the four distances, rHO, rNO, rHC, and rNC.
V′HB(rHO,θNHO,θHOC,φNHOC) = VHB(rHO,rNO,rHC,rNC) (8) For simplicity, in order to use effective two-body potentials, VHB was approximated as a sum of two-body potential functions.
VHB = hHO(rHO) + hNO(rNO) + hHC(rHC) + hNC(rNC) (9) For h, a 1-6-12 type function was used.
hij(rij) )
qiqj Bij Dij - 6 + 12 rij r r ij
(10)
ij
where the subscript ij represents one of the four atomic pairs, HO, NO, HC, and NC, respectively. Since the 1-4 interaction potential, hNC(rNC), does not contribute appreciably to VHB (as observed during the energy calculations), hNC was not included in the optimization of the parameters of VHB. Therefore, the optimum values of the B and D of the 1-2 pair, (H‚‚‚O), and of the 1-3 pairs (N‚‚‚O and H‚‚‚C), were determined to reproduce the ab initio PES of the hydrogen-bonded molecular pairs. Since B and D for the 1-4 pair were not included in the optimization, ANC and CNC of eq 3 were used for BNC and DNC, respectively, of eq 10.
17674 J. Phys. Chem., Vol. 100, No. 44, 1996
Kwon et al.
HO,NO,HC,NC
∑k
VHB )
hk(rk)
(11)
For VHO, the 1-2 pair in the hydrogen bond, Hagler et al.24 introduced both 1-6-9 and 1-6-12 type functions, and McGuire et al.7a used a 1-10-12 type function. c. Procedure for Optimization of Parameters. In this work, the potential parameters in VNB, i.e., the Aii’s, and in VHO, i.e., the Bk’s and Dk’s, were determined iteratively with two optimization procedures because these two parameter sets influence each other. The stabilization energy, Vst, at a configuration Rl, can be written as follows.
Vst ) Vst(Rl,{Aii},{Bk,Dk}) ) VNB(Rl,{Aii}) + VHB(Rl,{Bk,Dk}) + Vel(Rl) (12) For the optimization of the potential parameters, the two sets of parameters {Aii} and {Bk,Dk} were determined iteratively. In the first iterative step, {Aii} was fixed at {Aii(0)}, and {Bk,Dk} which satisfies the following conditions were obtained as (1) {B(1) k ,Dk }.
∂ ∂Rm
(
∑l
)
(n) |Vab(Rl) - Vst(Rl,{Aii(n)},{B(n) k ,Dk })|
|Vab(Rl)|
)0 for all Rm (13)
where the superscript (0), (1), and (n) represent initial, first, and nth iterative steps. Rm represents one of the B(n) k ’s and ab D(n) k ’s, and V (Rl) is the ab initio stabilization energy at Rl. (1) Next, in the second iterative step, {B(1) k ,Dk } were fixed to (2) obtain {Aii} as a {Aii }, which satisfies the following equilibrium conditions of the crystal. Nm
o,ka (n) (n) (n) FM({Aii(n)},{B(n) k ,Dk }) ) ∑[Wt|Ft ({Aii },{Bk ,Dk })| + (n)
ka
(n) o,cal - ELexpt| (14) Wτ|τo,ka({Aii(n)},{B(n) k ,Dk })|] + WL|EL
( ) ∂FM
∂Aii(n)
)0
for all ii
(15)
jj*ii
where Fto,ka and τo,ka represent the magnitude of the net translational force and net torque of the kath molecule in the crystal at its equilibrium position. The superscripts o and ka represent the equilibrium position and the kath molecule, respectively, in the asymmetric unit of the unit cell of the crystal. ELo,cal and ELexpt represent the calculated and experimental lattice energy of the crystal at 0 K. Wt, Wτ, and WL are weighting factors. Details of this procedure are described in our previous (3) paper.22 In the third iterative step, {B(3) k ,Dk } which satisfy eq (2) 13, were determined at fixed {Aii }. This iterative procedure continued until the equilibrium conditions of the crystal were satisfied within a certain error bound, and the ab initio PES are described well with the optimized parameters. The lattice energies and the information about the amino acid crystals which were used as constraints are summarized in Table 4. d. Test of the Reliability of the Potential Energy Function through Crystal Packing Studies. Crystal packing computations were carried out with the crystal packing computer program Lattice Minimization (LMIN).26 This program minimizes the energy of a crystal subject to no constraint other than the existence of a fully variable lattice. The molecules in the
TABLE 4: Experimental Crystal Information for Some Amino Acid Crystals and Their Lattice Energies Estimated from the Experimental Enthalpies of Sublimation19
crystal
experimental space group and number of molecules lattice energy (kcal/mol) ref in asymmetric unit
R-glycine zwitterion
P21/n, Z ) 4
68.0
L-alanine
P212121, Z ) 4, 23 K
65.5
Pbca, Z ) 8 P21/c, Z ) 4
70.2
zwitterion
β-alanine zwitterion DL-valine zwitterion
Tda L-cysteine zwitterion P212121, Z ) 4 a
17 25a 17 25b 25c 19 25d 25e
Td means tetragonal.
unit cell are treated as independent rigid bodies; energy minimization is carried out with the translational and rotational (Euler angle) coordinates of all molecules in the unit cell, and the six lattice parameters, as independent variables. Symmetry relations between the molecules in the unit cell are computed before and after energy minimization; comparison of these computed symmetries shows whether the initial space group symmetry was retained after energy minimization. In the determination of the parameters of the potential (section IIa), the molecular geometry and conformation were optimized. However, for the crystal packing studies, the molecular geometries of the molecules in the crystals were taken from the experimental crystal structures. The initial positions and attitudes of the molecules in the unit cell were calculated by applying the transformations of the appropriate space group to the coordinates of the asymmetric unit, as described in ref 26. The parameters A and B, which determine the cutoff distance and feather for nonbonded interactions,26 were set to 4.0 and 4.5, respectively. III. Results and Discussion The optimized hydrogen-bonding potential parameters are listed in Table 5. The depth, k, and the equilibrium distance, r°k, of the potential well of hHO, the 1-2 pair of the VHB, were obtained as 4.211 kcal/mol and 1.850 Å, respectively. The probability distribution of the O‚‚‚H bond distance, usually called the hydrogen-bond distance, in amino acid and peptide crystals was investigated by several workers.27,28 The resulting maximum probability was between 1.9 and 2.0 Å for r°max. The value of r°k ) 1.850 Å is a little shorter than r°max because both O and H atoms are located inside of the repulsive region of the 1-3 interaction potential functions. The maximum distribution range of rNO in the crystals was found to lie between 2.8 and 2.9 Å. However, r°NO of hNO was obtained as 3.902 Å. This means that the N+ atom is located inside the repulsive core of the oxygen atom, O-, and vice versa. The same situation arises for the C and H atomic pair. This is characteristic of the hydrogen-bond model introduced in this work to describe the angular dependence of the hydrogen-bonded system.19 The hydrogen bond is formed inside the repulsive cores of the 1-3 potential functions, hNO and hHC. The equilibrium distance of the nonbonded potential functions for the N‚‚‚O and H‚‚‚C atomic pairs are 3.145 and 3.120 Å, respectively. The repulsion between the 1-3 atomic pairs in the hydrogen-bonded system plays an important role in the description of the angular dependence of the hydrogen-bonded potential energy surface. The role of the 1-3 interactions in describing the angular dependence of the hydrogen-bond energy was discussed in our previous paper.19 The optimized nonbonded potential parameters are summarized in Table 6. The Cii’s were obtained from the effective
Potential Parameters for Amino Acid Zwitterions
J. Phys. Chem., Vol. 100, No. 44, 1996 17675
TABLE 5: Optimized Potential Parameters for 1-2 and 1-3 Atomic Pairs k of the Hydrogen-Bonded Systema atom pairb
Dkc (kcal Å12/mol)
Bkc (kcal Å6/mol)
kHB c,d (kcal/mol)
σkHB c,d (Å)
rko,HB c,e (Å)
rijo f (Å)
H5‚‚‚O5 N5‚‚‚O5 H5‚‚‚C5
6 760 896 086 58 517
337 508 82.4
4.211 0.072 0.029
1.648 3.476 2.987
1.850 3.902 3.353
3.145 3.120
a All parameters listed here pertain to all three complexes. b The atom types are described in column 1 of Table 2. c The potential parameters for the hydrogen-bonded potential, eq 10. d Calculated with the combination rules for Dk and Bk, similar to those of eq 4. e rko,HB ) 21/6σkHB. f The equilibrium distance of the nonbonded potential function for the atom pair i-j. σii and σjj were taken from Table 6, and the combination rule of eq 5 was used.
distances show some deviation from the ab initio results. Especially, the calculated rHO’s, 1.89-1.93 Å, are much longer than the ab initio values, 1.66-1.79 Å. The maximum probable value of rHO was observed between 1.9 and 2.0 Å in amino acid molecular crystals.28 Since both the ab initio PES and the equilibrium conditions of the molecular crystal were used as constraints, the calculated rHO’s are longer than the ab initio results and a little shorter than the average rHO distance in amino acid crystals; i.e., the optimized values are a compromise between these two constraints. The ab initio energies of the complexes are larger than those calculated with the empirical potential functions. The electrostatic energy makes the dominant contribution to the energy of the complex. If the potential parameters are optimized to reproduce the ab initio complex structures and energies, then the calculated lattice energies are overestimated, and the calculated cell volumes of the crystals become too small compared with the experimental results. To overcome this difficulty, this is one of the reasons why the crystal equilibrium conditions were also introduced as constraints together with the ab initio PES.
TABLE 6: Optimized Nonbonded Potential Parameters for the Atoms in the Model Compounds atom typea
Aiib (kcalÅ12/mol)
Ciic (kcalÅ6/mol)
iid (kcal/mol)
σiid (Å)
H1 H5 C1 C5 N5 O5
13 460 2 384 1 095 000 141 500 215 000 280 500
40.68 30.67 427.43 223.01 577.16 450.12
0.031 0.098 0.042 0.088 0.388 0.181
2.628 2.067 3.697 2.931 2.682 2.922
a From column 1 of Table 2. b Optimized in this work. c Calculated with eq 3 and 6. d Calculated from Aii and Cii using eq 4.
atomic polarizability by using eq 2. The effective atomic polarizabilities were calculated with the CDEAP method, eq 6. The Aii’s were optimized to satisfy the conditions described by eq 15. The Aii’s of H1 and C1 were taken from the crystal packing results of our previous work.22 Since the atoms in the ionic functional groups mostly participate in hydrogen bonds in the amino acid crystals, there is a little chance for those atoms to contact each other through van der Waals interactions. For this reason, one cannot expect to obtain accurate values of the repulsive nonbonded potential parameters for those atomic pairs from the equilibrium conditions of the crystal. In Table 7, the calculated minimum-energy conformations of the model complexes are compared with those of the ab initio MO results. The relative orientations of the complexes calculated with the empirical potential functions developed here agree well with the ab initio results. However, the hydrogen-bonded
In Table 8, the lattice constants, Euler angles, and the lattice energies, calculated with LMIN,26 for some amino acid crystals, are listed and compared with experimental data. Since the potential parameters for the hydrogen bond between the ionic groups in the zwitterion and the neutral functional groups, which are parts of the amino acid side chains, i.e., -COOH, -NH2, and -OH, were not yet determined, the crystal packing study
TABLE 7: Minimum-Energy Configurations of the Complexes Obtained Both with ab Initio (6-31G**) and with the Set of Empirical Potential Functions Developed Herea dimers methylammonium ion-acetate ion glycine zwitterion-methylammonium ion glycine zwitterion-acetate ion
ab initio empirical potential ab initio empirical potential ab initio empirical potential
rHO
rN‚‚‚O
θCNO
θNOC
1.66 1.89 1.79 1.93 1.78 1.92
2.57 2.72 2.72 2.77 2.73 2.78
133.81 140.63 131.64 138.53 47.19 47.55
91.32 90.77 93.57 90.42 94.13 90.01
a
Eelec
ENB
E1-3
EHB
-99.5
-0.1
9.7
-10.2
-37.6
-0.1
7.0
-9.3
-41.9
-0.1
8.0
-9.3
Eminb -116.32 -100.13 -47.38 -39.94 -45.63 -43.32
b
Distances are in angstroms, angles are in degrees, and energies are in kcal/mol. Emin is the binding energy at the minimum-energy configuration. Emin ) Edimer - (Emonomer1 + Emonomer2).
TABLE 8: Crystal Packing Results for Amino Acid Zwitterionsa lattice parameters molecular crystals R-glycine L-alanine
β-alanine DL-valine
Td L-cysteine
Elattice expb this work expb this work expb this work expb this work expb this work
-68.0 -64.8 -65.5 -63.1
Eelec -43.1 -43.0
-83.1 -70.2 -66.2
-63.5
-73.1
-50.6
-41.5
a
b
c
R
β
γ
volume
5.11 5.19 5.93 5.96 9.88 9.90 5.21 5.28 8.12 7.96
11.97 11.69 12.26 12.20 13.81 13.24 22.10 21.07 12.19 12.22
5.47 5.57 5.79 5.79 6.09 6.10 5.41 5.49 5.43 5.53
90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0
111.7 114.0 90.0 90.0 90.0 90.0 109.2 110.3 90.0 90.0
90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0 90.0
310 309 421 420 830 799 588 572 537 538
Euler angle 77.3 77.4 10.2 7.6 -98.4 -101.2 153.5 152.9 -105.4 -102.1
128.6 127.6 84.5 83.2 83.9 84.8 79.0 75.3 110.9 111.1
16.1 16.4 152.6 154.7 123.2 122.3 -132.7 -132.8 -165.3 -165.7
are in kcal/mol, the cell parameters a, b, and c are in Å, R, β, and γ and the Euler angles are in deg, and the volumes are in Å3. b The references to the experimental crystal structures are given in Table 4. aEnergies
17676 J. Phys. Chem., Vol. 100, No. 44, 1996 was limited to a few amino acid crystals in which only N-H+‚‚‚-O-C hydrogen bonds are formed. R-Glycine. The cell constants of R-glycine agree well with the experimental values and the relative orientations of the molecules in the crystal are also well reproduced. The error in the cell volume is about 0.3%. The lattice energy is estimated a little less well, 4.7%, compared with the estimated experimental lattice energy. L-Alanine. Since the errors in the cell lengths are less than 0.5%, and the lattice angles are the same as the experimental values, the cell volume agrees well with the experimental value. The alanines in the crystal are a little rotated, about 1-3 °, from their experimental positions. The lattice energies are also reproduced well. The contributions of the electrostatic energy to the lattice energy of R-glycine and L-alanine crystals are almost the same, ∼43 kcal/mol. β-Alanine. Since the cell parameter b was estimated less well compared with experiment, with an error of about 4%, the cell volume is also estimated less well by about 4%. The other cell parameters also agree well with the experimental values. The calculated lattice energy is much larger, 20 kcal/mol, than that of L-alanine. This extraordinary stabilization is due mainly to the contribution from the electrostatic interactions. DL-Valine. The cell parameters a and c are overestimated, about 0.15%, and b is estimated less well, about 4.5%. The cell angles R, β, and γ agree well with the experimental values. The cell volume is about 3% less than the experimental cell volume. The lattice energy is a little greater than the experimental value obtained with a Born-Haber cycle.21 The contribution of the electrostatic interaction is a little smaller than in the R-glycine and L-alanine crystals. Td L-Cysteine. The cell parameters are relatively well reproduced, and the cell volume is also reproduced well. The electrostatic energy contributes more than in R-glycine, L-alanine, and DL-valine. In general, the energy-minimized lattice constants, Euler angles, and lattice energies agree reasonably well with the experimental values. The calculated lattice energies are estimated less well in all the crystals. The sums of the contribution of the nonelectrostatic terms, hOH, hNO, hHC, hNC, and VNB, to the lattice energy were found to lie between 20 and 25 kcal/ mol. The set of empirical potential functions developed in this work reproduce the crystal structures well even though the potential functions from the ab initio minimum energy complex structures do not reproduce the dimer conformations very well, as shown in Table 7. The crystal structure itself does not contain any information about regions of the potential surface other than near the minimum-energy corresponding to the observed crystal structure. In particular, it provides no information about the angular dependence of the hydrogen bond. Gas-phase potentials from ab initio calculations, on the other hand, do provide this information about the angular dependence and can reproduce the structure of the complex, but not the crystal structure. Therefore, we have derived our empirical potential functions with a double constraint, Viz., to match the crystal structure and reproduce the ab initio structure of the complex. While our empirical potential does not reproduce complex structures very well, it does reproduce crystal structures and lattice energies very well. We regard the empirical potential function as a good one for polypeptides because the environment of a polypeptide in solution is much more like that in a crystal than the vacuum environment of the ab initio calculations. If the empirical potential parameters are adjusted to reproduce the ab initio minimum-energy complex structures well, the
Kwon et al. energy-minimized cell volumes of the crystals are reduced to about 80% of the experimental volumes. The cell volumes cannot be increased simply by increasing the repulsive core parameters of the nonbonded potential. If the σ of the LennardJones potential function is increased, then the molecules in the crystal rotate and become more compact with respect to their surroundings, the cell volume decreases, and the cell constants change in an unpredictable manner. Since the crystals do not contain much information about the angular dependence of the optimized hydrogen bond, during the optimization procedure, the angular dependence of the empirical potential energy function was mainly adjusted to the ab initio potential energy surface, and the bond distance was mainly adjusted to the information that comes from the crystals. IV. Conclusions A potential energy function for amino acid zwitterions has been developed. The electrostatic interaction energy was calculated with the point charges obtained from the M-PEOE method.18b The attractive coefficient of the r-6 term in the nonbonded potential was calculated with the Slater-Kirkwood formula.23 The atomic polarizabilities appearing in the SlaterKirkwood formula were calculated with the CDEAP method.21 The hydrogen-bond potential function for the zwitterions was also developed, based on our proposed hydrogen-bond model.19 The repulsive coefficient of the r-12 term in the nonbonded potential was obtained through crystal packing calculations on amino acid crystals. Both the crystal equilibrium conditions and the lattice energies were used as constraints for determining the repulsive parameters of the nonbonded potential. The set of potential energy functions does not reproduce the structure of ab initio complexes of the model compounds well. However, it does reproduce the crystal structures of some amino acid zwitterions. Since the purpose of this work is to develop potential energy functions for describing the interactions between amino acid zwitterions and the role of the amino acids in protein folding, we place stress on describing the crystal structures rather than the structures of the complexes in the hydrogen-bonding models. Acknowledgment. This work was supported by the Ministry of Science and Technology, Korea (N81540), by the Korea Science and Engineering Foundation (Korea-U.S. Cooperative Science Program), and by research grants from the National Science Foundation (MCB95-13167 and INT93-06345). The computations were carried out by using the Cornell National Supercomputer Facility, a resource of the Center for Theory and Simulation in Science and Engineering at Cornell University, which is funded in part by the National Science Foundation, New York State, the IBM Corporation, and members of its Corporate Research Institute, and also by using the computers of the Chemistry Department of Soong Sil University. References and Notes (1) Turner, A. J.; Whittle, S. R. Biochem. J. 1983, 209, 29. (2) Robert, E. In Bowery, N. G. Ed. Actions and Interactions of GABA and Benzodiazepines; Raven Press: New York, 1984; p 1. (3) Wolff, M., Ed. The Basis of Medicinal Chemistry. In Burger’s Medicinal Chemistry, 4th ed.; Wiley: New York, 1979. (4) Remko, M. J. Mol. Struct. (THEOCHEM) 1989, 201, 287. (5) (a) Clementi, E.; Cavallone, F.; Scordamaglia, R. J. Am. Chem. Soc. 1977, 99, 5531. (b) Scordamaglia, R.; Cavallone, F.; Clementi, E. J. Am. Chem. Soc. 1977, 99, 5545. (c) Bolis, G.; Clementi, E. J. Am. Chem. Soc. 1977, 99, 5550. (d) Carazzo, L.; Corongiu, G.; Petrongolo, C.; Clementi, E. J. Chem. Phys. 1978, 68, 787. (e) Corongiu, G.; Clementi, E. Gazz. Chim. Ital. 1978, 108, 273. (f) Clementi, E.; Corongiu, G.; Ranghino, G. J. Chem. Phys. 1981, 74, 578. (g) Ranghino, G.; Clementi, E.; Romano, S. Biopolymers 1983, 22, 1449.
Potential Parameters for Amino Acid Zwitterions (6) Ni, X.; Shi, X.; Ling, L. Int. J. Quantum Chem. 1988, 34, 527. (7) (a) McGuire, R. F.; Momany, F. A.; Scheraga, H. A. J. Phys. Chem. 1972, 76, 375. (b) Momany, F. A.; Carruthers, L. M.; McGuire, R. F.; Scheraga, H. A. J. Phys. Chem. 1974, 78, 1595. (c) Momany, F. A.; Carruthers, L. M.; Scheraga, H. A. J. Phys. Chem. 1974, 78, 1621. (d) Nemenoff, R. A.; Snir, J.; Scheraga, H. A. J. Phys. Chem. 1978, 82, 2504. (e) Kim, S.; Jhon, M. S.; Scheraga, H. A. J. Phys. Chem. 1988, 92, 7216. (f) Wee, S. S.; Kim, S.; Jhon, M. S.; Scheraga, H. A. J. Phys. Chem. 1990, 94, 1656. (g) Chipot, C.; Angyan, J. G.; Maigret, B.; Scheraga, H. A. J. Phys. Chem. 1993, 97, 9797. (8) (a) Jorgensen, W. L.; Swenson, C. J. J. Am. Chem. Soc. 1985, 107, 1489. (b) Jorgensen, W. L.; Gao, J. J. Phys. Chem. 1986, 90, 2174. (c) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657. (9) Bonaccorsi, R.; Palla, P.; Tomasi, J. J. Am. Chem. Soc. 1984, 106, 1945. (10) (a) Alagona, G.; Ghio, C.; Kollman, P. J. Am. Chem. Soc. 1986, 108, 185. (b) Singh, U. C.; Brown, F. K.; Bash, P. A.; Kollman, P. A. J. Am. Chem. Soc. 1987, 109, 1607. (c) Alagona, G.; Ghio, C.; Kollman, P. A. J. Mol. Struct. (THEOCHEM) 1988, 166, 385. (11) (a) Kokpol, S. U.; Doungdee, P. B.; Hannongbua, S. V.; Rode, B. M.; Limtrakul, J. P. J. Chem. Soc., Faraday Trans. 2 1988, 84, 1789. (b) Remko, M.; Liedl, K. R.; Rode, B. M. J. Mol. Struct. (THEOCHEM) 1995, 336, 7. (12) Rzepa, H. S.; Yi, M.-Y. J. Chem. Soc., Perkin Trans. 2 1991, 531. (13) Stamato, F. M. L. G.; Goodfellow, J. M. Int. J. Quantum Chem.: Quantum Biol. Symp. 1986, 13, 277. (14) Smith, P. E.; Dang, L. X.; Pettitt, B. M. J. Am. Chem. Soc. 1991, 113, 67. (15) Voogd, J.; Derissen, J. L.; van Duijneveldt, F. B. J. Am. Chem. Soc. 1981, 103, 7701. (16) (a) No, K. T.; Cho, K. H.; Kwon, O. Y.; Jhon, M. S.; Scheraga, H. A. J. Phys. Chem. 1994, 98, 10742. (b) Kwon, O. Y.; Kim, S. Y.; No, K. T. Bull. Korean Chem. Soc. 1995, 16, 410. (17) Power, L. F.; Turner, K. E.; Moore, F. H. Acta Crystallogr. 1976, B32, 11. (18) (a) No, K. T.; Grant, J. A.; Scheraga, H. A. J. Phys. Chem. 1990, 94, 4732. (b) No, K. T.; Grant, J. A.; Jhon, M. S.; Scheraga, H. A. J. Phys. Chem. 1990, 94, 4740.
J. Phys. Chem., Vol. 100, No. 44, 1996 17677 (19) No, K. T.; Kwon, O. Y.; Kim, S. Y.; Jhon, M. S.; Scheraga, H. A. J. Phys. Chem. 1995, 99, 3478. (20) Gaussian program, Revision D.2: Frisch, M. J.; Trucks, G. W.; Head-Gordon, M.; Gill, P. M. W.; Wong, M. W.; Foresman, J. B.; Johnson, B. G.; Schlegel, H. B.; Robb, M. A.; Replogle, E. S.; Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binkley, J. S.; Gonzalez, C.; Martin, R. L.;Fox, D. J.; Defrees, D. J.; Baker, J.; Stewart, J. J. P.; Pople, J. A. Gaussian, Inc., Pittsburgh, PA, 1992. (21) No, K. T.; Cho, K. H.; Jhon, M. S.; Scheraga, H. A. J. Am. Chem. Soc. 1993, 115, 2005. (22) No, K. T.; Kwon, O. Y.; Kim, S. Y.; Cho, K. H.; Yoon, C. N.; Kang, Y. K.; Gibson, K. D.; Jhon, M. S.; Scheraga, H. A. J. Phys. Chem. 1995, 99, 13019. (23) Scott, R. A.; Scheraga, H. A. J. Chem. Phys. 1966, 45, 2091. (24) (a) Hagler, A. T.; Lifson, S. J. Am. Chem. Soc. 1974, 96, 5327. (b) Hagler, A. T.; Leiserowitz, L.; Tuval, M. J. Am. Chem. Soc. 1976, 98, 4600. (c) Bernstein, J.; Hagler, A. T. J. Am. Chem. Soc. 1978, 100, 673. (d) Lifson, S.; Hagler, A. T.; Dauber, P. J. Am. Chem. Soc. 1979, 101, 5111. (e) Hagler, A. T.; Lifson, S.; Dauber, P. J. Am. Chem. Soc. 1979, 101, 5122. (f) Hagler, A. T.; Dauber, P.; Lifson, S. J. Am. Chem. Soc. 1979, 101, 5131. (25) (a) Jo¨nsson, P.-G.; Kvick, A° . Acta Crystallogr. 1972, B28, 1827. (b) Destro, R.; Marsh, R. E.; Bianchi, R. J. Phys. Chem. 1988, 92, 966. (c) Papavinasam, E.; Natarajan, S.; Shivaprakash, N. C. Int. J. Peptide Protein Res. 1986, 28, 525. (d) Mallikarjunan, M.; Rao, S. T. Acta Crystallogr. 1969, B25, 296. (e) Kerr, K. A.; Ashmore, J. P.; Koetzle, T. F. Acta Crystallogr. 1975, B31, 2022. (26) (a) Gibson, K. D.; Scheraga, H. A. J. Phys. Chem. 1995, 99, 3752. (b) Gibson, K. D.; Scheraga, H. A. J. Phys. Chem. 1995, 99, 3765. (27) (a) Vinogradov, S. N. Biopolymers 1979, 18, 1559. (b) Vinogradov, S. N. Int. J. Peptide Protein Res. 1979, 14, 281. (28) Go¨rbitz, C. H. Acta Crystallogr. 1989, B45, 390.
JP961180V
