16436
J. Phys. Chem. B 2010, 114, 16436–16442
Theoretical Studies of Salt-Bridge Formation by Amino Acid Side Chains in Low and Medium Polarity Environments Peter I. Nagy*,†,‡ and Paul W. Erhardt†,‡ Center for Drug Design and DeVelopment, The UniVersity of Toledo, Toledo, Ohio 43606-3390, United States, and Department of Medicinal and Biological Chemistry, The UniVersity of Toledo, Ohio 43606, United States ReceiVed: April 13, 2010; ReVised Manuscript ReceiVed: October 13, 2010
Salt-bridge formation between Asp/Glu · · · Lys and Asp/Glu · · · Arg side chains has been studied by model systems including formic and acetic acids as proton donors and methylamine, guanidine, and methylguanidine as proton acceptors. Calculations have been performed up to the CCSD(T)CBS//MP2/aug-cc-pvtz level with formic acid proton donors. Complexes formed with acetic acid were studied at the CCSD(T)/aug-cc-pvdz// MP2/aug-cc-pvdz level. Protein environments of low and moderate polarity were mimicked by a continuum solvent with dielectric constants (ε) set to 5 and 15, respectively. Free energy differences, ∆Gtot, were calculated for the neutral, hydrogen-bonded form and for the tautomeric ion pair. These values predict that a salt bridge is not favored for the Asp/Glu · · · Lys pair, even in an environment with ε as large as 15. In contrast, the Asp/Glu · · · Arg salt bridge is feasible even in an environment with ε ) 5. Charge transfers for the complexes were calculated on the basis of CHELPG and AIM charges. Introduction The protonation state of ionizable amino-acid side chains is an intricate problem in low polarity environments; namely, when imbedded within a protein. If the side chains of Asp, Glu, Arg, or Lys reside in the surrounding aqueous solution, they are expected to be mostly ionized. Determination of their protonation states becomes difficult, however, when the side chains are buried within the protein, where the effective dielectric constant (ε) of the environment is fairly low (e.g., ε ) 5-15). Calculation of absolute or relative pKa’s of small molecules has been the subject of many theoretical papers in the past decades.1-9 A more complicated problem is the determination of the protonation state of ionizable amino acid side chains in the depth of a protein. Interactions of these side chains with each other and considering the polarization effect of the environment has been studied using different theoretical methods based mainly on continuum solvent approaches.10-17 Recent experimental studies, including 1H, 13C, 15N NMR methods; UV/ vis spectrophotometry; and neutron crystallography, have provided data for comparison with the theoretical results.18-21 The general conclusions from the theoretical studies emphasize that geometric data and conformational issues as well as molecular mechanic force fields and parametrizations22-27 have sensitive effect on the obtained results. In some cases, the Asp or Glu side chain faces an Arg or Lys side chain in a protein such that two kinds of strong interactions are possible. If the residues adopt their neutral forms, then a hydrogen bond becomes favorable at a proper separation of the polar centers. An alternative interaction is possible when this complex adopts a so-called “salt-bridge” arrangement. A salt bridge corresponds to a hydrogen-bonded form, as well, but the partners are formally ionic (Scheme 1). This situation becomes important during molecular-mechanics-based modeling * Corresponding author. E-mail:
[email protected]. † Center for Drug Design and Development. ‡ Department of Medicinal and Biological Chemistry.
SCHEME 1: FA · · · MeNH2 (1a), FA- · · · MeNH3+ (1b), FA · · · Gua (2aa), FA · · · Gua (2ab), FA- · · · GuaH+ (2b), AA · · · MeNH2 (3a), AA- · · · MeNH3+ (3b), AA · · · MeGua (4a), AA- · · · MeGuaH+ (4b), Propanoic acid (5a), and Propanoate (5b)
in which the formal charges for the residues are preset and do not change throughout the simulation. The formal atomic charges are part of the applied force field and must be determined in harmony with other force-field parameters. In a force field by Duan et al.,28 atomic charges
10.1021/jp103313s 2010 American Chemical Society Published on Web 11/19/2010
Salt-Bridge Formation
J. Phys. Chem. B, Vol. 114, No. 49, 2010 16437
were assigned on the basis of a continuum dielectric solvent approach, in which ε was set to 4. This specific value was accepted as being relevant for dry protein powders. Asp, Glu, Arg, and Lys were considered with ionized side chains. Other studies, however, predicted a larger ε of 10-12 when water is embedded in the protein.29,30 Liljefors and Norrby31 performed MP4SDQ/6-311+G(3df,2p)// HF/6-31+G* gas-phase calculations, augmented with an SCIPCM/HF/6-31G* estimation of the solvent effect for the ion pair formation in a formic acid-trimethylamine complex. These authors found that ε > 9 is needed for predominant formation of the ion pair compared with the neutral hydrogen-bonded form. They also determined that the threshold value for ε is 4-6 when a single water molecule is present. The ion pair is not stable in the gas phase. Recently, the present authors studied hydrogen bonds between cyclic ethers, isopropyl alcohol, and phenol as one partner, with amino acid mimics as the other partner.32,33 Consideration of a dielectric continuum environment having ε ) 15 changed remarkably the association energies for the model complexes when compared with the gas phase. The calculations also called attention to the importance of the basis set effect as well as to the method used for the geometry optimization. The present study considers two levels of modeling for the interaction of acidic and basic amino acid side chains in protein environments. The simplest models are the formic acid · · · methylamine and the formic acid · · · guanidine systems. The relative free energies with ε ) 5 and 15 were calculated at the CCSD(T)/ aug-cc-pvdz and the CCSD(T)/CBS (complete basis set limit) levels, considering the MP2/aug-cc-pvdz and MP2/aug-cc-pvtz optimized geometries in solution, respectively. Because these results were similar, the acetic acid · · · methylamine and the acetic acid · · · methylguanidine systems were compared in solution at the CCSD(T)/aug-cc-pvdz//MP2/aug-cc-pvdz level. To our knowledge, these studies represent the highest-level calculations performed on such systems to date.
the converged wave function of the solute obtained from the in-solution calculation, and ∆Gdr is the relative dispersionrepulsion free energy. The cavity formation term was not considered because our calculations aim to model a protein environment wherein the space is already available for the amino acid residues. Nonetheless, the cavity formation energies with detailed solvation free energy contributions are provided in Tables S1 and S2 as Supporting Information. The immediate protein environment was modeled by an acetone solvent accounting for the dispersion and repulsion free energies. The dielectric constant of the solvent was set to 5 and 15 so as to mimic a low and a medium polarity environment, respectively. The relative thermal correction, ∆Gth, was calculated as
∆Gth ) ∆ZPE + ∆(H(T) - ZPE) - T∆S(T)
(3)
where ZPE is the zero-point energy; H(T) and S(T) are the enthalpy and the entropy, respectively, at T ) 298 K; and p ) 1 atm, calculated in the rigid-rotor, harmonic oscillator approximation.48 The CCSD(T) (coupled cluster for singles and doubles with noniterative triples)49 complete basis set in-solution relative internal energy was calculated using the formula by Hobza50 CCSD(T) MP2 ∆ECBS ) ∆ECBS + (∆ECCSD(T) - ∆EMP2)
(4)
The (∆ECCSD(T) - ∆EMP2) post-MP2 correction was calculated utilizing the aug-cc-pvdz basis set at the aug-cc-pvtz and augcc-pvdz optimized geometries for the FA and AA complexes, MP2 , the complete basis set limit MP2 energy respectively. ∆ECBS difference, was extrapolated by means of the formula51-53
E(CBS) ) E(X) - A/X3
(5)
Methods and Calculations The orientations of the partners in the complexes are shown in Scheme 1. All structures are hydrogen-bonded complexes and may be considered as tautomeric systems, depending on the position of the acidic proton. Henceforth, structures will be referred to as neutral or ion pair complexes when the proton resides on the acid or the base partner, respectively. Geometries of the complexes were determined by using the IEF-PCM approach, the integral-equation formalism for the polarizable continuum method.34-41 Applying scaled Bondi radii42 for the cavity formation as described previously,43 the geometry optimizations were performed at the MP244 level utilizing the aug-cc-pvdz and aug-cc-pvtz basis sets.45-47 The relative free energy for a tautomeric pair was calculated as s ∆Gtot ) (∆Eint + ∆Gth) + (∆Eelst + ∆Gdr) ≡ s + ∆Gth) + ∆G(solv) (∆Eint
(1)
s Eint ) < Ψ|H|Ψ>
(2a)
Eelst )
(2b)
wherein H is the solute’s Hamiltonian, V is the solvent reaction field generated by the fully polarized solute in solution, Ψ is
For calculating A and then E(CBS), aug-cc-pvdz single-point calculations were performed at the aug-cc-pvtz optimized geometries for the tautomeric FA complexes, and X was set to 2 and 3 with reference to the aug-cc-pvdz and aug-cc-pvtz energies, respectively. Calculations have been performed by means of the Gaussian 03 package.54 For the FA · · · Gua neutral complex, each of the 2aa and 2ab orientations correspond to a local energy minimum with ε ) 5. Structure 2ab was lower in energy by ∼4 kcal/mol at the IEFPCM/MP2/aug-cc-pvdz level, so 2aa was not further considered in the case when the dielectric constant was set to 5. However, calculations with ε ) 15 found only 2aa as a stable neutral complex, whereas optimization of 2ab led to a proton jump, and the structure moved toward 2b. Thus, the 2aa-type structure, confirmed by frequency analysis, was accepted for the neutral FA · · · Gua complex with ε ) 15. For complexes with acetic acid, no local energy minima were reached except for the neutral acetic acid · · · methylamine structures, even after very long IEF-PCM/MP2/aug-cc-pvdz minimizations. Relatively small imaginary frequencies in the range of 24i-48i were obtained, related to acetyl-methyl torsions. The optimizations stopped when the predicted energy changes were on the order of 10-7 au. Our computer resources did not allow tighter than default convergence criteria with the applied basis set. The main intermolecular bond lengths were stable up to about 0.001 Å when the MP2 energies varied on the 10-6 au scale. The lowest energy structures for the
16438
J. Phys. Chem. B, Vol. 114, No. 49, 2010
Nagy and Erhardt
TABLE 1: IEF-PCM/MP2 Optimized Intermolecular Geometric Parametersa FA · · · MeNH2 neutral
aug-cc-pvdz ε ) 5
H· · ·N O· · ·N 1.577
ion pair
O-H · · · N
O· · ·H O· · ·N
176.4
1.477
2.624 ε ) 15
1.537
176.4
1.563
176.8
O-H · · · N
169.3
1.512
173.8
2.576
1.536
171.2
176.8
2.585
1.526
174.1
2.582
1.426
171.6
2.545
1.522
1.476
174.2
2.550
1.505
171.6
2.593
1.518
173.6
2.572
AA · · · MeNH2
aug-cc-pvdz ε ) 5
ε ) 15
O· · ·H O· · ·N
O · · · H-N
H· · ·N O· · ·N
O-H · · · N
1.611
175.6
1.399
170.9
1.533
173.2
2.535 176.1
O · · · H-N 179.4 179.4 178.9 179.2 179.4 179.4 179.2 179.0
ion pairb
O-H · · · N
1.584 2.627
1.656 1.657 2.709 2.710 1.694 1.695 2.740 2.741 1.649 1.650 2.696 2.697 1.688 1.689 2.729 2.729
neutralb
H· · ·N O· · ·N
2.646
O· · ·H O· · ·N
AA · · · MeGua ion pairb
neutral
ion pair
O · · · H-N
2.619
2.610 ε )15
neutral H· · ·N O· · ·N
2.577
2.599 aug-cc-pvtz ε ) 5
FA · · · Gua
1.491 2.589
2.587
O· · ·H O· · ·N 1.628 1.643 2.686 2.698
O · · · H-N 179.4 179.8
171.2
a Distances in angstroms, angles in degrees. dO · · · H-N hydrogen bonds (see Scheme 1) with ε ) 5 for (1) FA · · · Gua/aug-cc-pvdz: O · · · H ) 1.917 Å, O · · · N ) 2.932 Å, O · · · H-N ) 170.8°; (2) FA · · · Gua/aug-cc-pvtz: O · · · H ) 1.896 Å, O · · · N ) 2.904 Å, O · · · H-N ) 170.1°; (3) AA · · · Gua/aug-cc-pvdz: O · · · H ) 1.892 Å, O · · · N ) 2.911 Å, O · · · H-N ) 172.6°. b Accepted from the structure with energy change on the 10-6 au scale in IEF-PCM/MP2/aug-cc-pvdz optimization.
complexes were accepted for further consideration, and the problem of estimating ∆ZPE and the relative thermal corrections in the absence of the lowest vibrational frequency is discussed in the next section. Charge transfers were calculated by deriving net atomic charges on the basis of the CHELPG fit55 to the respective molecular electrostatic potential in the gas phase and in solution and by applying the AIMAll software56 utilizing the quantum theory of atoms in molecules.57 Results and Discussion Complexes with Formic Acid. Our selected formic acid complexes (1a-2b) were considered as the smallest relevant model systems for studying the proton relocation in hydrogen bonds of the Asp/Glu · · · Lys and Asp/Glu · · · Arg amino acid side chains. In the very-minimum formic acid · · · ammonia complex, the important N-C bond is missing, which may lead to unreliable predictions. The Cδ-Nε bond is also missing in our FA · · · Gua system modeling the Asp/Glu · · · Arg side-chain interaction, but this simplification did not lead to a qualitative difference when compared with the acetic acid · · · methylguanidine (AA · · · MeGua) model (see below). On the other hand, calculations for the formic acid-containing systems allowed for a comparison of the CCSD(T)/aug-cc-pvdz and CCSD(T)CBS relative free energies, derived at the MP2/aug-cc-pvtz optimized geometries for the latter case. The main intermolecular geometric parameters for local energy minima structures are compared in Table 1. All results are from symmetry-unrestricted optimizations. The optimized
H · · · N separation decreases in the neutral FA · · · MeNH2 complex, and the O · · · H distance increases for the ionic FA- · · · MeNH3+ complex when the dielectric constant changes from ε ) 5 to 15. For the FA · · · Gua complex, both the H · · · N distance of the neutral form and the O · · · H distance for the ion pair increase with increasing ε. The O · · · H · · · N bond angles are slightly bent with a maximum deviation of ∼10° from a linear one. Comparing the two FA complexes, the most conspicuous difference has been found for the ion pairs. The O · · · H separations are larger by at least 0.15 Å in FA- · · · GuaH+ compared with that in the FA- · · · MeNH3+ pair. The O · · · H separations were calculated at a surprisingly low value for the FA- · · · MeNH3+ ion pair with the aug-cc-pvtz basis set. To check whether the obtained structures corresponded to local energy minima, frequency analysis were performed for the ionic species with ε ) 5 and 15. All frequency values were positive, and the corresponding ZPE and the Gth values differed only by up to 0.2 kcal/mol when calculated with the two basis sets. Relative energies/free energies are compared in Table 2. The basis set effect on the ∆Eint values is only 0.35 kcal/mol, a very promising result. It suggests that reliable relative energies may be obtained at the CCSD(T)/aug-cc-pvdz level, without the timeconsuming aug-cc-pvtz geometry optimization. ∆G(solv) is more negative by up to 0.6 kcal/mol with the aug-cc-pvdz compared with the aug-cc-pvtz basis set. Formic acid-containing model complexes predict positive and negative ∆Gtot for FA- · · · MeNH3+ and FA- · · · GuaH+, respec-
Salt-Bridge Formation
J. Phys. Chem. B, Vol. 114, No. 49, 2010 16439
TABLE 2: Energy Terms for the Ion Pair Relative to the Hydrogen-Bonded Neutral Complexa MP2/aug-cc-pvdz opt. ∆Eint/MP2 ∆Eint/CCSD(T) ∆G(solv) ∆Gth ∆Gtotb FA- · · · MeNH3+ ε)5 ε ) 15 FA- · · · GuaH+ ε)5 ε ) 15 AA- · · · MeNH3+ ε)5 ε ) 15 AA- · · · MeGuaH+ ε)5
8.59 10.31
9.46 11.04
-8.05 -11.63
3.75 -2.21
4.33 -1.71
-8.58 -7.52
8.77 11.28
9.73 12.08
-6.72 -10.90
3.11
3.57
-7.25
0.67 0.87
2.08 0.28
-0.02 -3.86c -0.14 -8.96c -0.08 1.00
2.93 2.18
-0.04 -3.72
MP2/aug-cc-pvtz opt. ∆Eint/MP2 ∆Eint/CCSD(T)CBS ∆G(solv) FA- · · · MeNH3+ ε)5 ε ) 15 FA- · · · GuaH+ ε)5 ε ) 15
∆Gtotd
8.47 10.28
9.53 11.21
-7.48 -11.05
2.72 1.03
3.52 -2.19
3.98 -1.66
-8.32 -7.42
-3.95c -8.81c
a Energies in kcal/mol. b ∆Gtot ) ∆Eint/CCSD(T) + ∆G(solv) + ∆Gth. c RT ln 2 ) 0.41 kcal/mol was added to ∆Gtot, accounting for the rotational symmetry number of 2 for the ion pair with the C2V symmetry. d ∆Gtot ) ∆Eint/CCSD(T)CBS + ∆G(solv) + ∆Gth.
tively. The projected conclusion is that the Asp/Glu · · · Lys side chains form neutral hydrogen bonds, whereas the ion pair is stable for the Asp/Glu · · · Arg system. Formic acid is, however, a stronger acid than either aspartic or glutamic acid, so it may behave as a too-strong proton donor. Furthermore, an N-C bond was replaced by an N-H in the Arg model, which may make ∆G(solv) unrealistically negative. Altogether, ∆Gtot for the FA- · · · GuaH+ model may be overly negative. Its relative stability of at least of 4 kcal/mol predicts almost exclusive presence of the ion pair in the system. ∆Gtot ) 0.3-1.0 kcal/ mol for FA- · · · MeNH3+ (ε ) 15), however, allows the presence of this form in a non-negligible fraction. For a better estimation of the prevailing tautomer, acetic acid-containing structures were also studied. Complexes with Acetic Acid. The acetic acid + methylamine complexes, AA · · · MeNH2 and AA- · · · MeNH3+ were considered with association modes of 3a and 3b, respectively. The aug-cc-pvdz intermolecular geometric parameters in Table 1 show that the H · · · N distance increases in the neutral form, whereas the O · · · H distance decreases with respect to the FA- · · · MeNH3+ ion pair. The O · · · H · · · N bond angles of the corresponding FA and AA complexes remain approximately unchanged. ∆Eint values (Table 2) increase slightly at both the MP2 and the CCSD(T) levels compared with that for the FA- · · · MeNH3+/ FA · · · MeNH2 model with ε ) 5; however, the increase in ∆Eint reaches 1 kcal/mol with ε ) 15. Furthermore, ∆G(solv) values are definitely less negative in the case of the AA · · · MeNH2 model. As a result, ∆Eint/CCSD(T) + ∆G(solv) increases from -0.59 to 1.18 kcal/mol (ε ) 15) when acetic acid is considered instead of formic acid in the model complex. In total, the ionic tautomer is less likely than estimated by means of the FA · · · MeNH2/FA- · · · MeNH3+ pair of models. Even a conservative estimate of ∆Gtot ) 1 kcal/mol at T ) 298 K allows only 15% for the salt-bridge form when the Asp/Gly · · · Lys side chains interact. Since the FA- · · · GuaH+ relative free energy was significantly negative with ε ) 15, study of the proton relocation seemed to
be critical only for the acetic acid · · · methylguanidine complexes in a model environment with ε ) 5. All relative energy terms for AA- · · · MeGuaH+ (4b) decreased in absolute value compared to the corresponding terms for the FA complex. Whereas decreases of the ∆Eint/MP2 and ∆Eint/CCSD(T) terms are not trivial, the change of ∆G(solv) from -8.58 to -7.25 kcal/mol is not surprising. Consideration of the methyl group on the guanidine moiety, a better model for the amino acid Arg, makes the system less polar; thus, the electrostatic component of ∆G(solv) had to become less negative (Table S1 of the Supporting Information). The ∆Gdr component is about zero, although the individual Gd and Gr components increased in absolute value as a consequence of the larger number of atoms in the system. Summing up all relative energy/free energy contributions, the calculated ∆Gtot is considerably negative, ensuring the prevalence of the ion pair form, and suggests the dominance of the ionic salt-bridge formation in proteins with interacting Asp/Glu · · · Arg side chains. Estimation of ∆ZPE and ∆Gth. As mentioned above, an imaginary frequency was obtained repeatedly for the optimized structure of each complex with acetic acid, being the only exception the neutral acetic acid · · · methylamine system. Visualization of the corresponding vibrations indicated that the problem is related to the rotation/vibration of the methyl group of acetic acid. Previous studies showed43,58-60 that the calculated wave numbers for the torsions of a methyl or a -NH3+ group change only slightly through conformational or tautomeric transitions. Since only the differences of the thermal corrections should be considered in those cases, a remarkable cancellation of the errors calculated for the individual terms were expected. The present tautomeric equilibrium, however, differs from those above. Table S4 of the Supporting Information shows the change of the 〈Ψ|H|Ψ〉 energy for the gas-phase acetate ion and the 〈Ψ|H + 0.5V|Ψ〉 energy for the ion in solution (ε ) 5, 15) as a function of the H-C-C-O torsion angle. The torsion potentials are of 6-fold symmetry; thus, they were considered only in the 0-30° torsion energy range. All geometric parameters were optimized, whereas the HCCO torsion angle was kept at the indicated value. The calculated barriers, V0, for the gas phase and for solutions with ε ) 5, 15 are about 0.010, 0.002, and 0.004 kcal/mol, respectively. At T ) 0 K, where all molecular motions are reduced to just vibration, the requirement of 0.5ω < V0 should be met. The corresponding wavenumber is ω < 3 cm-1 with ε ) 15. Not surprisingly, assignment of such a small value is very difficult with default convergence criteria in solution. Although no rotatable methyl groups exist in Asp or Glu, acetic acid mimic was used to these residues’ side chains to keep our modeling calculations tractable. A better approach is if the H-C-C-O torsion frequencies of acetic acid or the acetate ion are replaced by the corresponding values for the C-C-C-O moiety in propanoic acid (5a) and the propanoate ion (5b). Scheme 1 shows the optimized conformations for these molecules. The torsion frequencies were calculated by optimizing these species in solution. The lowest positive or the imiginary frequency for a given AA complex was replaced by the calculated value for the C-C-C-O moiety: 45.8 and 50.1 cm-1 (ε ) 5) for the ion and the neutral acid, respectively. The corresponding values are 50.5 and 50.6 cm-1 in the case of ε ) 15. Both ∆ZPE (Table S4) and ∆Gth (Table 2) were calculated accordingly. Charge Transfer. It has been clearly established that hydrogen bond formation entails an increase and decrease in
16440
J. Phys. Chem. B, Vol. 114, No. 49, 2010
Nagy and Erhardt
TABLE 3: Net Charges for the Elements of the Model Complexes in the Gas Phase and in Solutiona
TABLE 4: Derived Atomic Charges on the Basis of IEF-PCM/MP2/aug-cc-pvdz Calculations for Hydrogen-Bonded Complexes with Formic Acid
MP2/aug-cc-pvdz FA
MeNH2
FA-
MeNH3+
gas -0.177 0.177 ε ) 5 -0.250 0.250 -0.837 ε ) 15 -0.276 0.276 -0.871 ε ) 5 -0.110 0.110 -0.864 ε ) 15 -0.123 0.123 -0.884 AA
MeNH2
AA-
FA
FA-
Gua
neutral
CHELPG 0.837 0.871 AIM 0.864 0.884
MeNH3+
-0.099 0.099b -0.960 0.960 -0.091 0.091c -0.979 0.979 -0.104 0.104 -0.867 0.867 -0.121 0.121 -0.881 0.881 AA
AA-
MeGua
CHELPG -0.034 0.034 -0.777 d 0.802 -0.051 0.051b -0.897 0.851 AIM ε ) 5 -0.100 0.100 -0.842 0.842 -0.097 0.097 -0.861 ε ) 15 -0.108 0.108 -0.871 0.871
gas -0.159 0.159 ε ) 5 -0.225 0.225 -0.802 ε ) 15 -0.242 0.242 -0.851
MeGuaH+
0.777 0.897
0.861
MP2/aug-cc-pvtz FA
MeNH2
FA-
ε ) 5 -0.255 0.255 -0.821 ε ) 15 -0.283 0.283 -0.863
MeNH3+
0.821 0.863
CHELPG
GuaH+
FA
FA-
Gua
GuaH+
ε)5
ε ) 15
AIM ion pair
ε)5
ε ) 15
neutral ε)5
FA · · · NH2-Me O -0.609 -0.636 -0.740 -0.765 -1.218 C 0.694 0.719 0.744 0.769 1.650 O(H) -0.628 -0.651 -0.784 -0.806 -1.234 a H(O/N) 0.295 0.297 0.318 0.308 0.647 N -0.314 -0.277 -0.153 -0.101 -1.163 FA · · · guanidine O -0.705b -0.653c -0.855 -0.872 -1.226b C 0.741 0.716 0.800 0.818 1.660 O(H) -0.736 -0.738 -0.856 -0.870 -1.240 H(O/N)a 0.588 0.581 0.684 0.674 0.653 N -1.046 -1.027 -1.237 -1.233 -1.264 H(N) 0.539 0.678 0.678 0.535 N(H) -1.042 -1.219 -1.233 -1.308
ε ) 15
ion pair ε)5
ε ) 15
-1.227 -1.266 -1.274 1.651 1.666 1.670 -1.239 -1.250 -1.261 0.644 0.599 0.596 -1.166 -1.213 -1.213 -1.230c -1.265 -1.271 1.655 1.680 1.679 -1.246 -1.265 -1.271 0.654 0.595 0.587 -1.264 -1.353 -1.355 0.595 0.587 -1.353 -1.354
a Transferring proton. b Structure 2ab; see Scheme 1. c Structure 2aa (2aa structure with H · · · O ) separation of 3.41 Å); see Scheme 1.
-0.102 0.102 -0.954 0.954 -0.091 0.091c -0.972 0.972 b
a
Values in atomic charge units, geometries optimized at the indicated levels. b Structure 2ab; see Scheme 1. c Structure 2aa; see Scheme 1. d The structure is not stable in the gas phase.
the electron density on the proton donor and acceptor molecule, respectively.61-63 Umeyama and Morokuma64 relates this process mainly to electron transfer from bonding orbitals of the proton acceptor to antibonding orbitals of the proton donor. Calculated net charges for the species in Table 3 follow this pattern for the gas-phase complexes and for complexes in environments characterized with different dielectric constants by keeping in mind that the acid is the proton donor for neutral complexes, whereas the protonated base is the donor for ion pairs. Charge transfers for all systems were derived at geometries optimized at the indicated level. For formic acid complexes with the CHELPG charges, the corresponding net charges are very close at the MP2/aug-cc-pvdz and MP2/aug-cc-pvtz levels. This finding suggests that the MP2/aug-cc-pvdz basis set is satisfactorily large for estimating the transferred charges. The common feature of the data in this table is that the neutral acid gains, while the neutral base loses electron density. For the elements of the ion pair, the limiting charge value is (1 without electron transfer. Departures from these values indicate the direction and magnitude of the charge transfer. The solvent effect is consistent for the ion pairs: the larger the dielectric constant, the smaller the charge transfer: the ion charges are closer to (1 with ε ) 15 rather than ε ) 5. The transferred CHELPG charges for the FA complexes are smaller for the ion pairs than for the corresponding neutral complexes. The charge transfer is remarkably larger for the neutral MeNH2 compared with the Gua complex. The more polar solvent increases the charge transfer for the neutral FA · · · MeNH2 system but decreases it for the FA · · · Gua complex. Upon replacement of the FA component by AA, the charge transfer decreases and increases for the neutral and ionic complexes, respectively. The most striking feature of this table is the large amount of the transferred CHELPG charges. Grabowski61 characterized the hydrogen bonds of different strengths and found that the general amount of the transferred charge is equal to 0.01-0.03 units for typical complexes,63 although much larger values, up to 0.4 charge units, were calculated for some short and strong hydrogen
bonds.62 The 0.01-0.03 unit charge transfers were found in former studies for gas-phase complexes, in which the transferred amount was not calculated by means of CHELPG charges fitted to the MP2/aug-cc-pvdz or MP2/aug-cc-pvtz molecular electrostatic potentials. Nonetheless, the calculated charge transfers are surprisingly large for the neutral amine complexes. The solvent/protein environment leads to very considerable transfer of 0.15 to 0.24 units upon the modeled Asp/Glu · · · Lys interaction when AA is the acidic partner. A smaller, but non-negligible effect of the environment on the charge transfer, 0.05-0.10 e units, was calculated for the AA · · · MeGua complex as the Asp/ Glu · · · Arg model. These values are in contrast to charges of +1/-1 for Arg, Lys/Asp, Glu, as accepted by Duan et al. in their force field.28 The transferred charges scatter in the range of 0.10-0.16 units calculated upon the AIM atomic charges. The range with the CHELPG sets is 0.02-0.28 charge units (Table 3). This moderate difference may be surprising if the CHELPG and AIM atomic charges are compared for some selected atoms in the molecules. (Tables 4, S5). The AIM charges are almost always larger than the CHELPG values, and the difference may amount to ∼1.1 units (see N charges for the ionic FA · · · NH2Me complex in Table 4). The deviations in the atomic charges are always large for the elements of the carboxylic group. However, the differences in the transferred charges are much smaller. A remarkabe compensation must be in effect for the AIM charges with large absolute values but different signs, leading to only a relatively small net charge for the species. The derived charges reflect the philosophy behind the two methods. The CHEPLG procedure considers fitting points to the molecular electrostatic potential beyond the defined van der Waals distances from the individual atoms, whereas the AIM method determines atomic net charges within the region of the molecule belonging to each atom. The usefulness of the different charge sets in force fields could be evaluated on the basis of their performances in MD simulations, where comparisons with experimental data are possible. Aliste et al.65 also used +1/-1 charges for the above side chains in small peptides when studying salt-bridge formation in a lipid bilayer through MD simulations. The capacity of both Lys and Arg for forming salt bridges was pointed out, but the authors did not consider the possibility for the hydrogen-bonded
Salt-Bridge Formation neutral bonds. In a more recent publication by the Tieleman group, both the neutral and the ionized forms were considered for the above amino acids throughout their partitioning between the aqueous phase and a lipid bilayer.66 Moving toward the center of an ∼37 Å-wide, low-polarity bilayer, the free energy of the nonionized form becomes lower for Lys, Asp, and Glu. The free energies of the neutral and ionized forms of Arg were nearly equal in the bilayer at a depth of ∼20 Å. The explanation for all cases was that the surrounding water molecules, stabilizing the ionized form in aqueous solution, were gradually detached. The above results suggest that the neutral forms of the Arg, Lys, Asp, and Glu side chains not only could be present, but may exhibit the thermodynamically favored moiety in lowpolarity, nonaqueous environment. Presence of a small number of nearly residing water molecules may shift this equilibrium. Whether the corresponding pairs of the above amino acids form a salt bridge or a hydrogen-bonded moiety in the case of a favorable geometric position depends on the overall stabilizing effect of the environment. The present study suggests that different net charges should be allowed for these side chains if they can enter an intermolecular hydrogen bond. Application of altering charges in MD simulations is probably beyond the present-day practice. However, possible use of the AMOEBA force-field is a promising step in this direction.25,27 Conclusions The possible salt-bridge formation upon the interactions of the Asp/Glu · · · Lys and Asp/Glu · · · Arg side chains in proteins has been studied by model systems. Different polarization effects of different protein environments were mimicked by a dielectric continuum solvent with dielectric constants set to 5 and 15. Simpler models consisting of formic acid + methylamine and formic acid + guanidine allowed the estimation of the relative free energies for the neutral and ionic tautomers at both the CCSD(T)/aug-cc-pvdz//MP2/aug-cc-pvdz and the CCSD(T)CBS// MP2/aug-cc-pvtz levels. The similar results validated the capacity of the lower-level calculations and made theoretical considerations feasible for considering more representative models; namely, acetic acid + methylamine and acetic acid + methylguanidine. These latter complexes account for the C-Ccarboxyl and Cδ-Nε bonds for Asp/Glu and Arg, respectively. ∆Gtot for the more elaborate models predict that the salt bridge is thermodynamically not favored for the Asp/Glu · · · Lys pair, even in an environment with ε as large as 15. In contrast, the Asp/Glu · · · Arg salt bridge is preferred, even in an environment with ε as low as 5. These results are related to the inherent features of the complexes: the acetic acid + methylamine ion pair is not stable in the gas phase, in contrast to the acetate · · · methylguanidinium system. On the basis of the calculated remarkable charge transfers, the present study suggests that different net charges should be applied for the Asp, Glu, Lys, and Arg side chains if they can enter an intermolecular hydrogen bond. In closing, it should be noted that all of these results were obtained without considering the mediating effect of one or more water molecules. The possible presence of water molecule(s) around the interacting site of the side chains may need additional examination in future theoretical studies. Acknowledgment. The authors thank Dr. Todd A. Keith for valuable discussion related to the AIM charge derivation and
J. Phys. Chem. B, Vol. 114, No. 49, 2010 16441 the Ohio Supercomputer Center for the granted computer time. A USDA grant support provided to P.E. is also acknowledged. Supporting Information Available: Solvation free energy components, relative internal energy contributions and ZPE values are provided in Tables S1-S3; methyl rotation potentials, in Table S4 and Figure 1S; atomic charges in complexes with acetic acid, in Table S5. This information is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Lim, C.; Bashford, D.; Karplus, M. J. Phys. Chem. 1991, 95, 5610– 5620. (2) Richardson, W. H.; Peng, C.; Bashford, D.; Noodleman, L.; Case, D. A. Int. J. Quantum Chem. 1997, 61, 207–217. (3) Lee, I.; Kim, C. K.; Han, I. S.; Lee, H. W.; Kim, W. K.; Kim, Y. B. J. Phys. Chem. B 1999, 103, 7302–7307. (4) Klicic, J. J.; Friesner, R. A.; Liu, S.-Y.; Guida, W. C. J. Phys. Chem. A 2002, 106, 1327–1335. (5) Liptak, M. D.; Shields, G. C. J. Am. Chem. Soc. 2001, 123, 7314– 7319. (6) Liptak, M. D.; Gross, K. C.; Seybold, P. G.; Feldgus, S.; Shields, G. C. J. Am. Chem. Soc. 2002, 124, 6421–6427. (7) Kaminski, G. A. J. Phys. Chem. B 2005, 109, 5884–5890. (8) Lu, H.; Chen, X.; Zhan, C.-G. J. Phys. Chem. B 2007, 111, 10599– 10605. (9) Liu, S.; Pedersen, L. G. J. Phys. Chem. A 2009, 113, 3468–3655. (10) Bashford, D.; Case, D. A.; Dalvit, C.; Tennant, L.; Wright, P. E. Biochemistry 1993, 32, 8045–8056. (11) You, T. J.; Bashford, D. Biophys. J. 1995, 69, 1721–1733. (12) Antosiewicz, J.; McCammon, J. A.; Gilson, M. K. Biochemistry 1996, 35, 7819–7833. (13) Dillet, V.; Dyson, H. J.; Bashford, D. Biochemistry 1998, 37, 10298–10306. (14) Carlson, H. A.; Briggs, J. M.; McCammon, J. A. J. Med. Chem. 1999, 42, 109–117. (15) Stote, R. H.; Dejaegere, A. P.; Lefevre, J.-F.; Karplus, M. J. Phys. Chem. B 2000, 104, 1624–1636. (16) Dillet, V.; Van Etten, R. L.; Bashford, D. J. Phys. Chem. B 2000, 104, 11321–11333. (17) Khandogin, J.; Brooks, C. L., III Biochemistry 2006, 45, 9363– 9373. (18) Wang, J.; Rabenstein, D. L. Anal. Chem. 2007, 79, 6799–6803. (19) Andre, I.; Linse, S.; Mulder, F. A. A. J. Am. Chem. Soc. 2007, 129, 15805–15813. (20) Clark, A. T.; Smith, K.; Muhandiram, R.; Edmondson, S. P.; Shriver, J. W. J. Mol. Biol. 2007, 372, 992–1008. (21) Ishikawa, T.; Chatake, T.; Ohnishi, Y.; Tanaka, I.; Kurihara, K.; Kuroki, R.; Niimura, N. Chem. Phys. 2008, 345, 152–158. (22) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187–217. (23) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D., Jr; Pastor, R. W. J. Phys. Chem. B 2010, 114, 7830–7843. (24) Case, D. A.; Darden, T. A. ; Cheatham,, T. E.,III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; Wang, B. Hayik, S.; Roitberg, A.; Seabra, G.; Kolossva´ry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M.-J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Seetin, M. G.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. (2010), AMBER 11, University of California, San Francisco. (25) Ren, P. Y.; Ponder, J. W. J. Comput. Chem. 2002, 23, 1497–1506. (26) Xie, W.; Gao, J. J. Chem. Theory Comput. 2007, 3, 1890–1900. (27) For a recent review of force fields, see: Ponder, J. W.; Wu, C.; Ren, P.; Pande, V. S.; Chodera, J. D.; Schnieders, M. J.; Haque, I.; Mobley, D. L.; Lambrecht, D. S.; DiStasio, R. A., Jr.; Head-Gordon, M.; Clark, G. N. I.; Johnson, M. E.; Head-Gordon, T. J. Phys. Chem. B 2010, 114, 2549–2564. (28) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J.; Kollman, P. J. Comput. Chem. 2003, 24, 1999–2012. (29) Garcia-Moreno, B. E.; Dwyer, J. J.; Gittis, A. G.; Lattman, E. E.; Spencer, D. S.; Stites, W. E. Biophys Chem. 1997, 64, 211–224. (30) Dwyer, J. J.; Gittis, A. G.; Karp, D. A.; Lattman, E. E.; Spencer, D. S.; Stites, W. E.; Garcia-Moreno, B. Biophys. J. 2000, 79, 1610–1620. (31) Liljefors, T.; Norrby, P.-O. J. Am. Chem. Soc. 1997, 119, 1052– 1058. (32) Nagy, P. I.; Erhardt, P. W. J. Phys. Chem. A 2006, 110, 13923– 13932.
16442
J. Phys. Chem. B, Vol. 114, No. 49, 2010
(33) Nagy, P. I.; Erhardt, P. W. J. Phys. Chem. A 2008, 112, 4342– 4354. (34) Miertus, S.; Scrocco, E.; Tomasi, J. J. Chem. Phys. 1981, 55, 117– 129. (35) Tomasi, J.; Persico, M. Chem. ReV. 1994, 94, 2027–2094. (36) Cramer, C. J.; Truhlar, D. G. Chem. ReV. 1999, 99, 2161–2200. (37) Orozco, M.; Luque, F. J. Chem. ReV. 2000, 100, 4187–4225. (38) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. ReV. 2005, 105, 2999– 3094. (39) Cance`s, E.; Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 107, 3032–3041. (40) Cance`s, E.; Mennucci, B. J. Chem. Phys. 1998, 109, 249–259. (41) Cance`s, E.; Mennucci, B. J. Chem. Phys. 1998, 109, 260–266. (42) Bondi, A. J. Phys. Chem. 1964, 68, 441–451. (43) Nagy, P. I.; Alagona, G.; Ghio, C. J. Chem. Theory Comput. 2007, 3, 1249–1266. (44) Hehre, W. J.; Radom, L.; Schleyer, P. v. R. ; Pople, J. A. Ab Initio Molecular Orbital Theory; Wiley: New York, 1986. (45) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007–1023. (46) Kendall, R. A.; Dunning, T. H., Jr.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796–6806. (47) Peterson, K. A. Annu. Rep. Comput. Chem. 2007, 3, 195–206. (48) McQuarrie, D. A. Statistical Mechanics; University Science Books: Sausalito, CA, 2000. (49) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. Chem. Phys. Lett. 1989, 157, 479–483. (50) Hobza, P. Annu. Rep. Prog. Chem., Sect. C 2004, 100, 3–27. (51) Halkier, A.; Koch, H.; Jorgensen, P.; Christiansen, O.; Beck-Nielsen, I. M.; Helgaker, T. Theor. Chem. Acc. 1997, 97, 150–157. (52) Wilson, A.; Dunning, T. H., Jr. J. Chem. Phys. 1997, 106, 8718– 8726. (53) Helgaker, T.; Klopper, W.; Koch, H.; Noga, J. J. Chem. Phys. 1997, 106, 9639–9646. (54) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.;
Nagy and Erhardt Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, ReVision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (55) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361– 373. (56) Keith, T. A. AIMAll, Version 10.09.30; 2010; aim.tkgristmill.com. (57) Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Oxford University Press: Oxford, 1990. (58) Nagy, P. I.; Alagona, G.; Ghio, C. J. Am. Chem. Soc. 1999, 121, 4804–4815. (59) Alagona, G.; Ghio, C.; Nagy, P. I.; Durant, G. J. J. Phys. Chem. A 1999, 103, 1857–1867. (60) Alagona, G.; Ghio, C.; Nagy, P. I. J. Chem. Theory Comput. 2005, 1, 801–816. (61) Grabowski, S. J. Annu. Rep. Prog. Chem., Sect. C 2006, 102, 131– 165. (62) Grabowski, S. J.; Sokalski, W. A.; Leszczynski, J. J. Phys. Chem. A 2005, 109, 4331–4341. (63) Scheiner, S. Hydrogen Bonding, A Theoretical PerspectiVe; Oxford University Press: Oxford, 1997. (64) Umeyama, H.; Morokuma, K. J. Am. Chem. Soc. 1977, 99, 1316– 1332. (65) Aliste, M. P.; MacCallum, J. L.; Tieleman, D. P. Biochemistry 2003, 42, 8976–8987. (66) MacCallum, J. L.; Bennett, W. F. D.; Tieleman, D. P. Biophys. J. 2008, 94, 3393–3404.
JP103313S
