Like-Charge Ion Pairing in Water: An Ab Initio Molecular Dynamics

Letter pubs.acs.org/JPCL

Like-Charge Ion Pairing in Water: An Ab Initio Molecular Dynamics Study of Aqueous Guanidinium Cations Mario Vazdar,*,†,‡ Frank Uhlig,† and Pavel Jungwirth† †

Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic, Flemingovo nám. 2, CZ-16610 Prague 6, Czech Republic ‡ Division of Organic Chemistry and Biochemistry, Rudjer Bošković Institute, POB 180, HR-10002 Zagreb, Croatia ABSTRACT: The existence of like-charge guanidinium−guanidinium contact ion pairs in water is established by ab initio molecular dynamics simulations. Despite direct electrostatic repulsion, a contact ion pair is observed between two guanidinium cations, stabilized primarily by their amphiphilic behavior and van der Waals interactions. In a control simulation performed for two aqueous ammonium cations, no such contact ion pair is formed. This is the strongest computational evidence so far for the existence of specific contact ion pairing between guanidinium cations in water, with important implications for biological processes involving arginine-rich proteins.

SECTION: Biophysical Chemistry and Biomolecules

I

aqueous bulk continuum, this ion pair was calculated to be stabilized by about 0.09 eV.17 In order to obtain more conclusive computational evidence and characterize the electronic structural origin of the stability of the Gdm+− Gdm+ like-charge ion pair, we present here the results from extensive ab initio molecular dynamics (AIMD) simulations of a pair of guanidinium cations in ambient liquid water, contrasting its behavior to a control system consisting of two aqueous ammonium cations. AIMD simulations were performed for 2 guanidinium cations and 2 chloride anions dissolved in 62 water molecules in the unit cell. In parallel, simulations where the guanidiniums were replaced by two ammonium cations were also performed. Two different starting configurations were prepared; in the first one, Gdm+ cations were put in a stacked configuration at a distance of 3.4 Å and then hydrated. In the second one, guanidinium cations were placed further apart, at the distance of 6 Å, and in contrast to the first starting configuration, water molecules were placed also in between guanidinium cations. A similar starting condition was also prepared for the ammonium cations. Each of the trajectories was propagated for 50 ps with a time step of 0.5 fs. We employed the B-LYP density functional19 with an empirical dispersion correction (B-LYP-D),20 with the electronic density converged in each step (Born−Oppenheimer dynamics). As shown previously, dispersion is important for correct energetics of the interaction between guanidinium

on-specific effects play a vital role in living systems. They are important for essential biological functions, such as maintaining the electric potential between the cell exterior and interior1 and enzymatic activity in different ionic solutions.2,3 A very counterintuitive example of ion-specific interactions can be found in pairing between like-charged ions. It has been suggested previously that pairing between positively charged arginine residues is abundant in proteins4 and may be one of the driving forces of formation of protein dimers.5 Arginine−arginine pairing can also be among the factors responsible for binding between species B adenovirus and its receptor6 and for the enzymatic activity of cytochrome c oxidase.7 Direct experimental detection of like-charge aqueous ion pairs formed by two guanidinium cations (Gdm+), with Gdm+ being the constitutive component of the arginine side chains in proteins, remains elusive.8,9 Nevertheless, indirect experimental evidence, such as the slow down of electrophoretic mobility of oligoarginines by guanidinium cations,10 implies that likecharge Gdm+ ion pairs may indeed be stable in water. Likecharge ion pairing has also been studied by various computational methods, including molecular dynamics,11−14 Monte Carlo simulations,13,15 and quantum chemical calculations,16,17 which indicate that Gdm+ forms contact ion pairs in water. However, because the molecular dynamics simulations with different force fields show varying ion-pairing strength between the two guanidinium cations,18 it is important to address this issue also using ab initio approaches. A very recent quantum mechanical cluster study16 showed that the Gdm+−Gdm+ contact ion pair with 12 water molecules represents a weakly stable (by about 0.04 eV) minimal structure, while in an © 2012 American Chemical Society

Received: June 13, 2012 Accepted: July 18, 2012 Published: July 18, 2012 2021

dx.doi.org/10.1021/jz3007657 | J. Phys. Chem. Lett. 2012, 3, 2021−2024

The Journal of Physical Chemistry Letters

Letter

cations.16 Pseudopotentials21 replaced the core electrons. The Kohn−Sham orbitals were expanded into an atom-centered double-ζ quality DZVP-MOLOPT-GTH basis set 22 in conjunction with an auxiliary plane wave basis set with a cutoff of 280 Å for the electronic density.23 Periodic boundary conditions were used with a fixed unit cell of dimensions of 12.65 Å × 12.65 Å × 12.65 Å. Simulations were performed at 300 K24 using the CP2K program package.25 For the purpose of testing the reliability of the B-LYP-D method used in AIMD simulations, we first benchmarked it against more accurate perturbative MP2 and coupled cluster CCSD(T) calculations on small Gdm+−water clusters. These calculations were performed using the Gaussian09 program.26 Results for the complexation energy ΔEc between Gdm+ and one or three water molecules using different basis sets are presented in Table 1. The comparison in Table 1 shows that

forming a contact ion pair with no water molecules between them. Water molecules (and occasionally also chloride counterions) form hydrogen bonds with Gdm+ exclusively in the molecular plane of the ion, being repelled from the hydrophobic faces of Gdm+. In the control case of two ammonium cations (right panel), no such contact ion pair is formed, with the other NH4+ being separated by water molecule(s). The relatively high occurrence of the ammonium solvent-separated ion pair is, to a large extent, an artifact due to the periodicity of a relatively small system, which does not allow for separations of the two ions larger than half of the unit cell size in each dimension. The present system is, nevertheless, large enough to safely distinguish between contact and solventseparated ion pairs, with guanidinium (but not ammonium) cations clearly preferring the former arrangement. The difference in co-ion pairing and the hydration structure of Gdm+ versus NH4+ ions can be further inferred from radial distribution functions (RDFs) with respect to the carbon atom of Gdm+ and nitrogen atoms of NH4+, presented in Figure 2. For a pair of Gdm+, the peak in the RDF is located at 3.7 Å, which is the same distance as that at which also the peak in the RDF for water oxygen atoms is located (left panel). This implies that there is no water in between paired guanidinium ions. These RDFs thus represent a clear signature of formation of a like-charge Gdm+ contact ion pair. The situation is completely different for the ammonium ion pair where the first peak in the RDF of central nitrogen atoms only appears at a distance of 4.7 Å. Because the first peak in the RDF between NH4+ and water oxygens occurs at 2.9 Å, only solvent-separated NH4+−NH4+ ion pairs exist, with no formation of contact pairs. The RDFs in Figure 2 are statistically significant and meaningful; nevertheless, AIMD simulations may be computationally too demanding to allow for a fully converged sampling of the distribution of distances between the two guanidinium cations. In order to check the potential effect of the choice of initial geometry, two very different starting configurations of the two guanidinium ions were prepared. Starting from the first initial geometry, where the two guanidinium ions were in a stacked contact ion pair, this close contact was maintained throughout the whole 50 ps of the simulation. In case of the other initial configuration, where the guanidinium ions were separated by water molecules, a T-shape Gdm+ ion pair was formed first. This T-shaped ion pair with an average central carbon distance of 4.7 Å persisted for another 20 ps (Figure 3), after which a stacked contact ion pair was formed and persisted. This was accompanied with shortening of the carbon−carbon distance and energy stabilization (Figure 3). Accordingly, the angle between the two Gdm+ planes steadily decreased from the value of about 80°, which corresponds to the T-shape structure, to values much closer to 0, which in turn correlate with the stacked structure. Note that the existence of the Tshape structure as a local minimum on the potential energy surface was also predicted by model calculations of a guanidinium ion pair in a dielectric continuum representing water, with the stacked contact ion pair being a global minimum.17 Recently, a very useful analysis of noncovalent interactions (NCIs) has been developed based on the electronic density and its reduced gradient s = [1/(2(3π2)1/3)(|∇ρ|/ρ4/3)].27 This approach is based on the fact that the reduced gradient describes deviations from the homogeneous electronic distribution.28 Low-density and low-gradient regions in a system, together with the sign of a second-derivative

Table 1. Complexation Energies ΔEc for the Guanidinium Cation and One Water Molecule (Gdm + 1W) and the Guanidinium Ion Pair with Three Water Molecules (given in kcal mol−1), Obtained by B-LYP-D and MP2 Geometry Optimizations with Different Basis Sets, Respectivelya Gdm-Gdm + 3Wb

Gdm + 1W ΔEc

B-LYP-D

MP2

CCSD(T)

B-LYP-D

MP2

cc-pVDZ aug-cc-pVDZ cc-pVTZ aug-cc-pVTZ

−25.3 −17.4 −20.2 −16.9

−22.2 −17.9 −19.4 −17.7

−21.8 −17.8 −19.3 −17.6

51.5 56.6 54.7 57.1

55.5 54.8 56.7 57.2

a

CCSD(T) single-point energy calculations were performed at the MP2 geometry with a given basis set. bTaken from ref 16.

augmentation of the basis set by diffuse functions is essential for proper description of the interaction between the Gdm+ cation and surrounding water molecules. Thus, we chose to use in the production runs the DZVP-MOLOPT-GTH basis set, which is comparable in size to the aug-cc-pVDZ basis set and also contains diffuse functions.22 The comparison of complexation energies ΔEc shown in Table 1 indicates that the BLYP-D/ DZVP-MOLOPT-GTH level of theory is adequate for the purpose of the present study. Spatial distribution functions of co-ions and water around the Gdm+ or NH4+ ion from aqueous bulk AIMD simulations are depicted in Figure 1. In case of Gdm+ (left panel), another Gdm+ ion is found with high occurrence in close proximity,

Figure 1. Spatial distribution functions around the Gdm+ ion (left; Gdm+, blue; water oxygen, red; and water hydrogen, white) and NH4+ (right; NH4+, blue; water oxygen, red; and water hydrogen, white). The black arrow points to a water molecule(s) occupying space between two NH4+. 2022

dx.doi.org/10.1021/jz3007657 | J. Phys. Chem. Lett. 2012, 3, 2021−2024

The Journal of Physical Chemistry Letters

Letter

Figure 2. (Left) Radial distribution functions of the central carbon atom of the Gdm+ ion versus another central carbon atom of the Gdm+ ion (green), the water oxygen atom (red), and the water hydrogen atom (black) around the central carbon atom of the Gdm+ ion. (Right) Radial distribution functions of the central nitrogen atom of the NH4+ ion versus another central nitrogen atom of the NH4+ ion (blue), the water oxygen atom (red), and the water hydrogen atom (black) around the central nitrogen atom of the NH4+ ion. Simulations for the Gdm+ ion pair were started from the stacked configuration.

Figure 3. Time evolution of the angle (green) between the two guanidinium ions and the total potential energy ΔEpot (blue). Plotted values are running averages averaged over 2000 fs. The AIMD simulation was started from a solvent-separated Gdm+−Gdm+ ion pair, and the first unequilibrated 10 ps were discarded from the plot. Representative snapshots of a T-shaped Gdm+ and stacked Gdm+ contact ion pair with characteristic distances between central carbon atoms are presented. The transition between these two structures is indicated by a dashed vertical line.

Figure 4. Gradient isosurfaces (s = 0.6 au) for Gdm+ (upper panel) and NH4+ ion pairs (lower panel) are displayed for the selected AIMD snapshots. The surfaces are colored on a color scale according to values of sign(λ2) × ρ (shown on the right), ranging from −0.05 to +0.05 au.

simulation. The key behind the unusual stability of this “electrostatics-defying” structure5 lies in the amphiphilic nature of the guanidinium (Gdm+) cation. Due to its planarity, aromaticity, and nonuniform charge distribution, the faces of the Gdm+ cation are hydrophobic and weakly hydrated, with guanidinium forming hydrogen bonds with water molecules exclusively in the plane.29 This also leads to further stabilization by van der Waals interactions, as evidenced by the present noncovalent interaction analysis. It is important to note that in the case of the water-separated initial configuration, a T-shape structure was formed first. This structure represents a local minimum on the potential energy surface, with an energetically more favorable stacked structure being formed subsequently and being then stable for the rest of the simulation. Formation of like-charge contact ion pairs requires specific geometric and electronic structure properties (planarity and a delocalized inhomogeneous electronic distribution), such as that present for guanidinium cations. “Normal” ions, such as the control ammonium cations, clearly do not form a like-charged contact ion pair but at best solvent-separated pairs, the occurrence of which was enhanced by the finite size of the simulated cell. The

component of the electronic density λ2, are then identified with various types of noncovalent interactions, that is, hydrogen bonds, van der Waals (vdW) interactions, and sterical clashes. These interactions differ in the value of sign(λ2) × ρ, which is positive for steric clashes, negative for hydrogen bonds, and slightly negative for vdW interactions. We performed the NCI analysis for a set of snapshots from the AIMD for aqueous Gdm+ and NH4+ ion pairs (Figure 4). In case of the Gdm+ ion pair, the green color of the isosurfaces, which indicates vdW interactions, is present in all of the displayed snapshots. These interactions get stronger upon closer contacts of the stacked guanidinium cations (compare snapshots at 30 vs 45 ps). In contrast to the guanidinium ion pair, green color (or any other color) isosurfaces between NH+4 cations practically do not occur, which further points to the fact that like-charge contacts are not stabilized for the control system of the pair of aqueous ammonium cations. The present extensive AIMD simulations demonstrate that a stacked like-charge ion pair forms between two guanidinium cations and, once formed, is stable for the duration of the 2023

dx.doi.org/10.1021/jz3007657 | J. Phys. Chem. Lett. 2012, 3, 2021−2024

The Journal of Physical Chemistry Letters

Letter

present findings, based on an electronic structure description of a bulk aqueous guanidinium solution at ambient conditions, provide the most direct and strongest computational evidence so far for the remarkable and counterintuitive formation of likecharge Gdm+−Gdm+ contact ion pairs in water. Establishing this phenomenon with a solid computational as well as experimental backing will allow us to better understand and rationalize biological processes involving interactions within and between arginine-rich proteins.

■

(13) Bader, J. S.; Chandler, D. Computer Simulation Study of the Mean Forces between Ferrous and Ferric Ions in Water. J. Phys. Chem. 1992, 96, 6423−6427. (14) Zangi, R. Attraction between Like-Charged Monovalent Ions. J. Chem. Phys. 2012, 136, 184501. (15) Boudon, S.; Wipff, G.; Maigret, B. Monte Carlo Simulations on the Like-Charged Guanidinium−Guanidinium Ion Pair in Water. J. Phys. Chem. 1990, 94, 6056−6061. (16) Vazdar, M.; Vymetal, J.; Heyda, J.; Vondrasek, J.; Jungwirth, P. Like-Charge Guanidinium Pairing from Molecular Dynamics and Ab Initio Calculations. J. Phys. Chem. A 2011, 115, 11193−11201. (17) Vondrasek, J.; Mason, P. E.; Heyda, J.; Collins, K. D.; Jungwirth, P. The Molecular Origin of Like-Charge Arginine−Arginine Pairing in Water. J. Phys. Chem. B 2009, 113, 9041−9045. (18) Wernersson, E.; Heyda, J.; Vazdar, M.; Lund, M.; Mason, P. E.; Jungwirth, P. Orientational Dependence of the Affinity of Guanidinium Ions to the Water Surface. J. Phys. Chem. B 2011, 115, 12521−12526. (19) Becke, A. D. Density Functional Exchange Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098−3100. (20) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H−Pu. J. Chem. Phys. 2010, 132, 154104. (21) Goedecker, S.; Teter, M.; Hutter, J. Separable Dual-Space Gaussian Pseudopotentials. Phys. Rev. B 1996, 54, 1703−1710. (22) VandeVondele, J.; Hutter, J. Gaussian Basis Sets for Accurate Calculations on Molecular Systems in Gas and Condensed Phases. J. Chem. Phys. 2007, 127, 114105. (23) Lippert, G.; Hutter, J.; Parrinello, M. A Hybrid Gaussian and Plane Wave Density Functional Scheme. Mol. Phys. 1997, 92, 477− 487. (24) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (25) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. QUICKSTEP: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach. Comput. Phys. Commun. 2005, 167, 103−128. (26) Frisch, M. J. et al. Gaussian 09, revision A.1.; Gaussian Inc.: Wallingford, CT, 2009. (27) Johnson, E. R.; Keinan, S.; Mori-Sanchez, P.; Contreras-Garcia, J.; Cohen, A. J.; Yang, W. Revealing Noncovalent Interactions. J. Am. Chem. Soc. 2010, 132, 6498−6506. (28) Cohen, A. J.; Mori-Sanchez, P.; Yang, W. Insights into Current Limitations of Density Functional Theory. Science 2008, 321, 792− 794. (29) Mason, P. E.; Neilson, G. W.; Dempsey, C. E.; Barnes, A. C.; Cruickshank, J. M. The Hydration Structure of Guanidinium and Thiocyanate Ions: Implications for Protein Stability in Aqueous Solution. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 4557−4561.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Support from the Czech Science Foundation (Grant 203/08/ 0114) and the Academy of Sciences is gratefully acknowledged. F.U. thanks the International Max Planck Research School for Dynamical Processess in Atoms, Molecules and Solids for support.

■

REFERENCES

(1) Berg, J. M.; Tymoczko, J. L.; Stryer, L. Biochemistry; W. H. Freeman: New York, 2002. (2) Heyda, J.; Pokorna, J.; Vrbka, L.; Vacha, R.; Jagoda-Cwiklik, B.; Konvalinka, J.; Jungwirth, P.; Vondrasek, J. Ion Specific Effects of Sodium and Potassium on the Catalytic Activity of HIV-1 Protease. Phys. Chem. Chem. Phys. 2009, 11, 7599−7604. (3) Wondrak, E. M.; Louis, J. M.; Oroszlan, S. The Effect of Salt on the Michaelis Menten Constant of the HIV-1 Protease Correlates with the Hofmeister Series. FEBS Lett. 1991, 280, 344−346. (4) Magalhaes, A.; Maigret, B.; Hoflack, J.; Gomes, J. N. F.; Scheraga, H. A. Contribution of Unusual Arginine-Arginine Short-Range Interactions to Stabilization and Recognition in Proteins. J. Protein Chem. 1994, 13, 195−215. (5) Pednekar, D.; Tendulkar, A.; Durani, S. Electrostatics-Defying Interaction between Arginine Termini As a Thermodynamic Driving Force in Protein-Protein Interaction. Proteins: Struct., Funct., Bioinf. 2009, 74, 155−163. (6) Sumowski, C. V.; Schmitt, B. B. T.; Schweizer, S.; Ochsenfeld, C. Quantum-Chemical and Combined Quantum-Chemical/MolecularMechanical Studies on the Stabilization of a Twin Arginine Pair in Adenovirus Ad11. Angew. Chem., Int. Ed. 2010, 49, 9951−9955. (7) Qian, J.; Mills, D. A.; Geren, L.; Wang, K.; Hoganson, C. W.; Schmidt, B.; Hiser, C.; Babcock, G. T.; Durham, B.; Millett, F.; Ferguson-Miller, S. Role of the Conserved Arginine Pair in Proton and Electron Transfer in Cytochrome c Oxidase. Biochemistry 2004, 43, 5748−5756. (8) Mason, P. E.; Neilson, G. W.; Enderby, J. E.; Saboungi, M. L.; Dempsey, C. E.; MacKerell, A. D., Jr; Brady, J. W. The Structure of Aqueous Guanidinium Chloride Solutions. J. Am. Chem. Soc. 2004, 126, 11462−11470. (9) Hunger, J.; Niedermayer, S.; Buchner, R.; Hefter, G. Are Nanoscale Ion Aggregates Present in Aqueous Solutions of Guanidinium Salts? J. Phys. Chem. B 2010, 114, 13617−13627. (10) Kubickova, A.; Krizek, T.; Coufal, P.; Wernersson, E.; Heyda, J.; Jungwirth, P. Guanidinium Cations Pair with Positively Charged Arginine Side Chains in Water. J. Phys. Chem. Lett. 2011, 2, 1387− 1389. (11) Pettitt, B. M.; Rossky, P. J. Alkali Halides in Water: Ion−Solvent Correlations and Ion−Ion Potentials of Mean Force at Infinite Dilution. J. Chem. Phys. 1986, 84, 5836−5844. (12) Dang, L. X.; Montgomery Pettitt, B. A Theoretical Study of Like Ion Pairs in Solution. J. Phys. Chem. 1990, 94, 4303−4308. 2024

dx.doi.org/10.1021/jz3007657 | J. Phys. Chem. Lett. 2012, 3, 2021−2024