CationâÏ Interactions between Quaternary Ammonium Ions and

Article Cite This: J. Phys. Chem. B 2018, 122, 2251−2260

pubs.acs.org/JPCB

Cation−π Interactions between Quaternary Ammonium Ions and Amino Acid Aromatic Groups in Aqueous Solution Esam A. Orabi† and Guillaume Lamoureux* Department of Chemistry and Biochemistry and Centre for Research in Molecular Modeling (CERMM), Concordia University, 7141 Sherbrooke Street West, Montréal, Québec H4B 1R6, Canada S Supporting Information *

ABSTRACT: Cation−π interactions play important roles in the stabilization of protein structures and protein−ligand complexes. They contribute to the binding of quaternary ammonium ligands (mainly RNH3+ and RN(CH3)3+) to various protein receptors and are likely involved in the blockage of potassium channels by tetramethylammonium (TMA+) and tetraethylammonium (TEA+). Polarizable molecular models are calibrated for NH4+, TMA+, and TEA+ interacting with benzene, toluene, 4-methylphenol, and 3methylindole (representing aromatic amino acid side chains) based on the ab initio MP2(full)/6-311++G(d,p) properties of the complexes. Whereas the gas-phase affinity of the ions with a given aromatic follows the trend NH4+ > TMA+ > TEA+, molecular dynamics simulations using the polarizable models show a reverse trend in water, likely due to a contribution from the hydrophobic effect. This reversed trend follows the solubility of aromatic hydrocarbons in quaternary ammonium salt solutions, which suggests a role for cation−π interactions in the salting-in of aromatic compounds in solution. Simulations in water show that the complexes possess binding free energies ranging from −1.3 to −3.3 kcal/mol (compared to gas-phase binding energies between −8.5 and −25.0 kcal/mol). Interestingly, whereas the most stable complexes involve TEA+ (the largest ion), the most stable solvent-separated complexes involve TMA+ (the intermediate-size ion).

1. INTRODUCTION

Cation−π interactions involving quaternary ammonium ions have not received as much attention as those involving metal ions.19−31 Ammonium (NH4+) is in itself a physiologically relevant ion known to bind proteins32,33 but can also be viewed as a minimal model compound for the Lys side chain (RNH3+). Tetramethylammonium (TMA+) is a model for ligands with an RN(CH3)3+ cationic moiety, such as trimethyllysine (one of the methylation products of methyltransferases)34,35 or acetylcholine (known to bind, through its quaternary ammonium head, to the π electrons of a Tyr residue in acetylcholinesterase).16 TMA+ inhibits the activity of most potassium channels, as does the larger tetraethylammonium (TEA+).36,37 Ab initio calculations on the complexation of quaternary ammonium ions and model compounds to the side chains of Phe, Tyr, and Trp are thus important for understanding the strength and directionality of these interactions and for calibrating potential models of cation−π interactions in proteins and in protein−ligand complexes. Although cation−π interactions are strong in the gas phase,19,20,38 their strength is reduced by the successive addition of water molecules to the complex.39 The strength of these interactions is significantly reduced in aqueous solutions.19,20,40−45 Interestingly, although cation−π interac-

Cation−π interactions refer to the noncovalent association between an inorganic cation or the cationic moiety of an organic molecule and the π electrons of alkenes, alkynes, or aromatics.1 In proteins, cation−π interactions occur between the ammonium group of lysine (Lys) or the guanidinium group of arginine (Arg) and the aromatic side chain of phenylalanine (Phe), tyrosine (Tyr), or tryptophan (Trp).2,3 In protein−ligand complexes, they occur either between aromatic amino acids and cations or cationic moieties (such as K+, NH4+, or RNH3+), or between aromatic ligands and the cationic side chains of Lys and Arg. Cation−π interactions in proteins are found more commonly than expected from chance alone.3−7 They are also often found at the protein− protein8 and protein−DNA interfaces.9,10 They contribute to protein stability,11−14 protein−ligand interactions,2,15 and molecular recognition in general.16 The stability of cation−π complexes in gas phase is predominantly due to electrostatic and polarization forces.2,17,18 Contributions from other forces such as dispersion and charge transfer are much smaller. Electronic polarization is important because the cations usually found in these complexes are small and produce strong electric fields, and because the aromatic ligands are electron-rich and highly polarizable.17,19,20 For these interactions, polarizable force fields have been demonstrated to be more accurate than pairwise-additive ones.17−20 © 2018 American Chemical Society

Received: December 5, 2017 Revised: February 1, 2018 Published: February 4, 2018 2251

DOI: 10.1021/acs.jpcb.7b11983 J. Phys. Chem. B 2018, 122, 2251−2260

Article

The Journal of Physical Chemistry B

2.2.1. Parametrization Strategy. Optimization of potential models for the various cation−π and cation−water pairs follows our previously reported approach.19,20,32,56,57 In particular, the optimization consists of modifying the van der Waals potential energy between interacting fragments by adjusting pair-specific Lennard-Jones (LJ) parameters between atom i in the ion and atom j in the aromatic compound, Emin,ij and Rmin,ij, defined by the Lorentz−Berthelot combination rules

tions are typically not as strong as cation−anion interactions in the gas phase, they are less destabilized by the successive addition of water molecules and retain their attractive character even in aqueous solutions. Calculation of the binding affinity between cations and aromatic systems in aqueous solutions is an important step toward understanding their strength and directionality in proteins. Empirical force fields are often used for such purpose,19,20,40−44 given the prohibitive computational cost of ab initio simulations.45 Force fields that correctly describe the various interactions in the system (ion−water, ion−aromatic, and water−aromatic) are however required. The aim of this work is to computationally investigate the binding affinity and binding directionality between quaternary ammonium ions (NH4+, TMA+, and TEA+) and aromatic amino acid side chains (modeled by benzene, toluene, 4methylphenol, and 3-methylindole) in aqueous solution. This effort includes the development of empirical force field parameters for cation−π interactions in the context of the Drude polarizable force field.46,47 The results of our simulations are useful to understand the nature of quaternary ammonium binding in proteins and to understand the phenomenon of salting-in of aromatic hydrocarbons by quaternary ammonium salts.

Emin, ij =

Emin, i × Emin, j and R min, ij =

R min, i + R min, j 2 (1)

Optimization initially targets the ab initio PECs of the interacting pair. Each point of the energy surfaces contributes to the error function χ2 by a Boltzmann-weighted error term.19 Because the PECs are computed using rigid monomer geometries extracted from the optimized dimers, parameters obtained from the minimization of χ2 are further refined to reproduce the complexation energies of the fully relaxed ab initio pairs.19 Pair-specific LJ parameters for NH4+ interactions have been reported previously.19,32 In the present work, no pair-specific LJ parameters are found necessary for TMA+− water interaction. For the interaction of TMA+ with the aromatic ligands (benzene, toluene, 4-methylphenol, and 3methylindole), pair-specific LJ parameters between the nitrogen atom of the ion and the carbon atoms of the sixmembered ring are optimized. For TEA+ interactions, pairspecific LJ parameters are optimized between the methylene carbon atoms of the ion and the oxygen atom of water or the carbon atoms of the six-membered rings of the aromatic ligands. Models for TMA+−benzene and TEA+−benzene are optimized first. The transferability of the optimized parameters to the other ion−aromatic complexes is then tested. For complexes where the ion−benzene parameters are considered nontransferable (if the complexation energies of the complexes are deviating by more than 5% from ab initio data), the parameters are readjusted. For TEA+ complexes, the value of the Emin parameter between the methylene carbon atoms of the ion and the carbon atoms of the six-membered ring is found not to be critical and is set to the ion−benzene optimized value. All other mixed LJ parameters are derived from the Lorentz−Berthelot combination rules. This strategy, achieved using the NBFIX facility of the CHARMM force field,46,47 avoids modifying the original electrostatics of the compounds (partial charges and polarizabilities) and preserves their solvation properties. We will refer to the new polarizable model as “LJ-optimized”. 2.2.2. Molecular Dynamics (MD). All MD simulations are performed with cubic periodic boundary conditions in the isothermal−isobaric ensemble (NpT). Water molecules are described with the polarizable SWM4-NDP water model.52 In all simulations, the water molecules are totally rigid (OH bonds and HOH angles) and all covalent bonds to hydrogen atoms are made rigid using the SHAKE/Roll and RATTLE/ Roll algorithms.58 All simulations are carried out using one ion, one aromatic molecule, and 500 water molecules at T = 298.15 K and p = 1 atm. The electrostatic interactions are computed using the particle-mesh Ewald method,59 with κ = 0.34 for charge screening and a 1.0 Å grid spacing with sixthorder splines for mesh interpolation. Real-space interactions (Lennard-Jones and electrostatic) are cut off at 15 Å and the long-range contribution from the Lennard-Jones term is

2. METHODS 2.1. Ab Initio Calculations. The geometries of NH4+, TMA+, and TEA+ in complex with water, benzene, toluene, 4methylphenol, and 3-methylindole are optimized without the symmetry constraints at the MP2(full)/6-311++G(d,p) level using Gaussian 09.48 The interaction energies are corrected for basis set superposition error (BSSE) by the counterpoise method,49 and referred to as ECP throughout the text. Frequency calculations are performed on all structures to confirm that they are energy minima. The optimization of each pair involves various initial structures, with different binding modes of the ion (uni, bi, and tridentate) and different orientations of the ion with respect to the plane of the aromatic moiety, and only the global minimum structures are considered. Potential energy curves (PECs) of all interacting pairs are calculated at the same level and all potential energies are corrected for BSSE. The curves are calculated by scanning the distance between the nitrogen atom of the ion and its point of projection in the aromatic ring plane (or the oxygen atom of water, in the case of ion− water complexes) between 2 and 10 Å. During the scan, the geometry of each molecule in the dimer and the relative orientation of the two fragments is kept as in the gas-phase MP2(full)/6-311++G(d,p) global minimum optimized geometry of the dimer. 2.2. Molecular Mechanics Calculations. Molecular mechanics calculations are performed with the program CHARMM.50 Polarizable potential models based on classical Drude oscillators51 have been reported for H2O,52 NH4+,19 and aromatic compounds.53,54 Following our previous work,19,20,32,56,57 a polarizable model for interaction of NH4+, TMA+, and TEA+ with water, benzene, toluene, 4methylphenol, and 3-methylindole is optimized based on the ab initio properties of their cation−π and cation−water pairs. The polarizable force field parameters for the isolated compounds55 are reported in the Supporting Information (SI). 2252


Article

The Journal of Physical Chemistry B introduced as an average density-dependent term.60 A dual Nosé−Hoover thermostat61 is applied; the first is coupled to the atomic sites to keep them at 298.15 K and the second is coupled to the Drude particles to keep them at low temperature (1 K) and ensure self-consistent dipole induction.51 The relaxation time of the atomic thermostat is 0.1 ps, whereas that of the Drude thermostat is 0.005 ps. An Andersen−Hoover barostat62 with a relaxation time of 0.2 ps is used to regulate the pressure. The equations of motion are integrated using a 1 fs time step. 2.2.3. Potential of Mean Force Calculations. Potentials of mean force (PMFs) between each cation (NH4+, TMA, or TEA+) and each aromatic compound (benzene, toluene, 4methylphenol, or 3-methylindole) are calculated using umbrella sampling. The distance r between the nitrogen atom of the ion and the center of the six-membered ring of the aromatic compound is used as a reaction coordinate and a harmonic potential of force constant of 10 kcal/(mol Å2) is applied to bias the sampling. The reaction coordinate is sampled between 2.0 and 12.0 Å using 0.5 Å separated windows, and each window is simulated for 2.5 ns, the first 0.5 ns of which is considered equilibration and discarded. Simulations are performed either by letting both solutes free to rotate (unconstrained PMFs), or by constraining them in the so called “en face” and “edge-on” orientations. The en face geometry is enforced by constraining all C−X···N angles to 90° (where C is a carbon atom of the aromatic sixmembered ring, X is the geometric center of the ring, and N is the nitrogen atom of the ion) using a force constant of 100 kcal/(mol rad2). The edge-on geometry is enforced by constraining the center of the ion to be in the plane of the aromatic moiety of the ligand, using a force constant of 100 kcal/(mol Å2). In all simulations, the aromatic molecule is constrained to the center of the simulation box using a 2 kcal/(mol Å2) force constant. Although the aromatic molecules are let free to rotate in the unconstrained and en face constrained simulations, they are kept in a fixed plane through all of the edge-on constrained simulations. The unbiased PMFs are reconstructed using the weighted histogram analysis method,63,64 and the radial variation in the entropy of the solute pairs is taken into account by adding a 2RT ln(r) correction term to the PMFs from the unconstrained and en face simulations and by adding an RT ln(r) correction term to the PMF from the edge-on simulations.

Figure 1. Optimized geometries of TEA+ at the MP2(full)/6-311+ +G(d,p) level of theory. The methyl carbon atoms are in a quasiplanar configuration in conformer 1a and form an asymmetric pyramid in conformer 1b. Conformer 1a has a C2v point symmetry, whereas conformer 1b has a C1 symmetry. Atomic coordinates are provided in the SI.

optimized and default Drude models (EDrude), are reported in Table 1. Table 1 also reports the equilibrium distance between the nitrogen atom of the ion and the geometric center of the six-membered ring of the aromatic compound, Drude for both ab initio (rQM ) eq ) and Drude force field (req calculations. For the three water complexes and for the NH4+−4-methylphenol complex, the distance is between N and O. The MP2(full)/6-311++G(d,p) results for NH4+containing pairs are in close agreement with our previous results obtained at the MP2(FC)/6-311++G(d,p) level.19,32 For each of the five ligands, the binding energies across the ion series follow the trend NH4+ ≫ TMA+ > TEA+. NH4+ binds water in a unidentate fashion.19 Except for 4methylphenol complex (2d), in which NH4+ interacts with oxygen, complexes of NH4+ (2a, 2b, 2c, and 2e) have cation−π geometries, with the ion located on top of the sixmembered ring (en face geometry).19,32 Complex 2d (ECP = −23.80 kcal/mol) possesses an alternative stable structure (not shown) in which the ammonium ion interacts with the π electrons, with ECP = −19.74 kcal/mol, E = −22.55 kcal/mol, and rQM eq = 2.904 Å. In the minimum-energy conformation, TMA+ binds water in a tridentate fashion. Complexes of TMA+ (2f−j) display cation−π geometries in which the ion is slightly displaced from the center of the six-membered ring. Complexes of TEA+ (2l−o) possess binding energies close to those of the corresponding TMA+ complexes (see Table 1). The pyramidal conformer (1b) of TEA+ forms a complex with benzene (structure not shown) that is slightly more stable than complex 2l. This alternative complex is characterized by ECP = −8.98 kcal/mol, E = −14.79 kcal/mol, and rQM eq = 4.313 Å. The BSSE-corrected binding energy between NH4+ and benzene (ECP = −17.55 kcal/mol) and between TMA+ and benzene (ECP = −8.97 kcal/mol) are in agreement with the experimental binding enthalpy in gas phase (−17.122,38 and −9.4 kcal/mol,66 respectively). To the best of our knowledge, no experimental binding enthalpy is available for the TEA+− benzene complex. Ab initio potential energy curves for NH4+, TMA+, and TEA+ in complex with water, benzene, toluene, 4-methylphenol, and 3-methylindole are reported in Figure 3 (solid lines), along with the corresponding curves obtained from the default (dashed lines) and LJ-optimized (dotted lines) Drude models (see Section 3.2). 3.2. Optimized Force Field. Pair-specific LJ parameters for the interaction of NH4+ with water, benzene, toluene, 4methylphenol, and 3-methylindole have been optimized previously.19,32 These are reported in Table 2, along with the newly optimized pair-specific LJ parameters for the

3. RESULTS AND DISCUSSION 3.1. Ab Initio Interaction Energies and Potential Energy Curves. TEA+ is reported to possess two main conformers (Figure 1; see the Supporting Information for coordinates).65 The nitrogen and methyl carbon atoms form a quasi-planar configuration in conformer 1a (termed planar), whereas the methyl carbon atoms form an asymmetric pyramid in conformer 1b (termed pyramidal).65 Optimization of the two structures at the MP2(full)/6-311++G(d,p) level shows that conformer 1a is 1.17 kcal/mol more stable than conformer 1b and thus only conformer 1a is considered for complexes of TEA+. The optimized geometries of the global minimum conformers of all of the studied complexes are presented in Figure 2 (see the SI for coordinates). The BSSE-corrected and uncorrected complexation energies (E CP and E, respectively), as well as the binding energies from the LJ2253


Article


Figure 2. Optimized geometries at the MP2(full)/6-311++G(d,p) level of theory for the complexes of NH4+ with water (2a), benzene (2b), toluene (2c), 4-methylphenol (2d), and 3-methylindole (2e), and for the complexes of TMA+ (2f−2j) and TEA+ (2k−2o) with the same ligands. See the SI for atomic coordinates.

Table 1. Ab Initio Complexation Energies (E: Uncorrected; ECP: BSSE-Corrected) and Equilibrium Distances (rQM eq ) Calculated at the MP2(full)/6-311++G(d,p) Level of Theory, and Corresponding Energies (EDrude) and Distances (rDrude ) eq Calculated Using the Default and LJ-Optimized Polarizable Modelsa,b ab initio ligand

ion

conformer

H2O

NH4+ +

2a 2f 2k 2b 2g 2l 2c 2h 2m 2d 2i 2n 2e 2j 2o

benzene

toluene

4-methylphenol

3-methylindole

TMA TEA+ NH4+ TMA+ TEA+ NH4+ TMA+ TEA+ NH4+ TMA+ TEA+ NH4+ TMA+ TEA+

ECP

E −22.25 −11.65 −10.34 −19.89 −12.90 −13.90 −21.59 −14.15 −15.44 −26.38 −16.95 −18.17 −28.07 −21.28 −22.98

rQM eq

−20.29 −9.78 −7.97 −17.55 −8.97 −8.49 −18.90 −9.99 −9.53 −23.80 −12.62 −11.70 −25.04 −15.38 −14.88

2.703 3.632 3.882 2.917 4.204 4.325 2.920 4.166 4.273 2.641 4.177 4.276 2.822 4.088 4.199

Drude, LJ-optimized

Drude, default

EDrude

EDrude

−20.27 −9.86c −8.25 −17.55 −8.93 −8.77 −19.06 −10.02 −9.76 −21.00 −13.10 −11.76 −24.60 −15.48 −14.98

rDrude eq 2.805 3.565c 3.801 2.892 4.301 4.373 2.901 4.283 4.431 2.702 4.138 4.451 3.016 4.152 4.366

−19.07 −9.86 −10.40 −19.45 −8.14 −9.20 −20.68 −9.31 −10.83 −17.90 −10.77 −13.09 −33.54 −15.16 −16.70

rDrude eq 2.934 3.565 3.776 2.864 4.134 4.292 2.858 4.135 4.272 2.905 4.138 4.290 2.690 4.041 4.137

a Energies are in kcal/mol and distances are in Å. bFor water−ion and 4-methylphenol−NH4+ complexes, req refers to the distance between the N atom of the ion and the O atom of the ligand. For all of the other ligands, it refers to the distance between the N atom and the center of the sixmembered ring. cNo optimization was required for the H2O−TMA+ complex.

interaction of TMA+ and TEA+ with the five ligands. The parameters are first optimized based on the ab initio PECs and then refined to reproduce the ab initio geometry and interaction energy in the global minimum complex (see columns “ECP” and “EDrude” for the optimized model of Table 1). The LJ parameters obtained from the Lorentz−Berthelot combination rules underestimate the stability of NH4+ complexes with water and 4-methylphenol but overestimate

that of the other ammonium complexes (Table 1 and Figure 3). Besides reproducing the ab initio interaction energy, the optimized parameters for H2O−NH4+ interaction have been previously shown to reproduce the experimental hydration free energy of the ion.19 It was also shown that the interaction of NH4+ with the four aromatics is better modeled by including a nonatomic site (X) at the center of the sixmembered ring.19 This site is found unnecessary for the TMA+ and TEA+ complexes. In line with MP2 results, the 2254


Article

The Journal of Physical Chemistry B Figure 3. continued

atom of the ion and the O atom of water (rNO) is scanned between 2 and 10 Å. For complexes with aromatic ligands, the distance between the nitrogen atom of the ion and its point of projection on the aromatic ring plane is scanned between 2 and 10 Å. For consistency, the curves for the aromatic complexes are shown by plotting ECP versus the distance between N and the center of the sixmembered ring (rNX). In all curves, the monomers are kept in the gas-phase internal geometries found in the dimers and the distance is scanned in 0.1 Å increments, keeping the relative orientations of the monomers as in the optimized complex.

optimized model for the interaction between NH4+ and 4methylphenol predicts a second stable structure possessing a cation−π geometry, with EDrude = −21.52 Å and rDrude = 3.237 eq Å. In comparison, the default LJ parameters severely overestimate the ab initio binding energy of this cation−π complex (EDrude = −24.24 Å versus ECP = −19.74). The properties of the TMA+−water complex are adequately described by the LJ parameters obtained from the combination rule, and no pair-specific parameters were introduced. The default LJ parameters slightly underestimate the binding energies for the TMA+−aromatic complexes, especially with 4-methylphenol (Table 1 and Figure 3). The optimized pair-specific LJ parameters for the interaction of TMA+ with benzene are transferable to the interaction of the ion with toluene and 3-methylindole. The LJ parameters obtained from the Lorentz−Berthelot combination rules overestimates the ab initio binding energies for all TEA+ complexes (Table 1 and Figure 3), and thus the model is calibrated. For complexes of TEA+, the methylene carbons rather than the N atom of the ion are used to adjust the pair-specific LJ parameters (see Table 2). The optimized parameters for the TEA+−toluene interaction are transferred to the interaction of the ion with 4-methylphenol. In line with MP2 results, the complexation of benzene with the pyramidal conformer (1b) of TEA+ gives a more stable structure than with the planar conformer (1a). The model predicts a complexation energy of −9.93 kcal/mol and rDrude = 4.411 Å. eq As shown in Figure 3, the models reproduce the ab initio PECs as a whole. 3.3. Cation−π Interactions in Aqueous Solution. The binding affinity of NH4+, TMA+, and TEA+ with the four aromatic ligands in water is estimated from the PMF calculations. The PMFs from the unconstrained and the constrained en face and edge-on orientations are shown in Figure 4. (The positions and depths of the PMF minima are reported in Table S1 of the SI.) The results show that, for a given aromatic ligand, the binding free energy from the unconstrained and en face constrained simulations follows the trend TEA+ > TMA+ > NH4+, which is the reverse of that of the gas-phase binding energies (see Table 1). For a given ligand, the binding free energies from the edge-on constrained simulations are comparable across the ion series, and no systematic trend is observed. These binding free energies are much smaller (less negative) than those from the en face constrained simulations, showing that the en face binding conformation is favored in water as in gas phase (Figure 2). The equilibrium separations between the ion and the center of the six-membered ring are similar for both the unconstrained and en face simulations. These distances are slightly larger than for the gaseous pairs (by about 0.4 Å; see Table

Figure 3. Potential energy curves for NH4+ (black), TMA+ (red), and TEA+ (blue) in complex with (a) water, (b) benzene, (c) toluene, (d) 4-methylphenol, and (e) 3-methylindole calculated at the MP2(full)/6-311++G(d,p) level (solid curves), as well as with the default (dashed curves) and LJ-optimized (dotted curves) Drude models. For complexes with water, the distance between the nitrogen 2255


Article


Table 2. Optimized Pair-Specific Lennard-Jones Parameters for the Interactions of NH4+, TMA+, and TEA+ with Water, Benzene, Toluene, 4-Methylphenol, and 3-Methylindolea pair-specific LJ parameters N(NH4+) molecule

i

H2O benzene

O Xb Cc Xb Cc Xb O Cc Xb Cc

toluene 4-methylphenol

3-methylindole

C(TEA+)d

N(TMA+)

H(NH4+)

Emini,N (kcal/mol)

Rmini,N (Å)

Emini,H (kcal/mol)

Rmini,H (Å)

0.10185 0.14706 − 0.14706 − 0.01872 0.3975 − 0.6632 0.3105

3.7592 3.50006 − 3.50006 − 4.69984 3.3358 − 2.9908 3.73231

0.00924 0.00602 − 0.00602 − 0.0 0.01844 − 0.02714 0.00836

2.88481 3.28084 − 3.28084 − 0.0 2.56374 − 2.77158 3.49493

Emini,N (kcal/mol)

Rmini,N (Å)

− 0.0 0.30971 0.0 0.30971 0.0 − 0.36967 0.0 0.30971

− 0.0 4.80186 0.0 4.80186 0.0 − 4.18721 0.0 4.80186

Emini,C (kcal/mol)

Rmini,C (Å)

0.01349 0.0 0.06216 0.0 0.06216 0.0 − 0.06216 0.0 0.06216

5.01795 0.0 4.30961 0.0 4.42961 0.0 − 4.42961 0.0 4.62561

a Parameters for NH4+ are reproduced from refs 19 and 32. bRefers to a nonatomic site at the center of the six-membered ring (see refs 19, 20, and 32 for details). cRefers to the carbon atoms of the six-membered ring. dRefers to the methylene carbon atoms in TEA+. The dashes in the table indicate that no pair-specific LJ parameters are optimized, and that the corresponding LJ parameters are obtained using the Lorentz−Berthelot combination rules.

Å) in some of the edge-on PMFs of TMA+ and TEA+ correspond to solvent-separated pairs. For the NH4+−benzene complex, the values of −3.3 and −0.9 kcal/mol for en face and edge-on binding affinities are consistent with the values of −3.2 and −0.8 kcal/mol obtained from a two-dimensional PMF reported previously.20 In comparison, Sa et al.,45 based on a fully quantum mechanical simulation of the NH4+−benzene pair in water, have found an en face binding free energy of −5.75 kcal/mol (at 3.25 Å) and an edge-on binding free energy of −0.33 kcal/mol. For the NH4+−toluene complex, Chipot et al.,42 using an additive force field corrected for cation−π interactions, have reported an en face binding free energy of −5.47 kcal/mol (at 3.05 Å) and an average PMF minimum of −2.99 kcal/mol (at 3.16 Å). The en face free energy wells found by Sa et al. and by Chipot et al. are located at almost exactly the same distances as those from Figure 4e (3.25 vs 3.2 Å for benzene and 3.05 vs 3.1 Å for toluene) but are significantly deeper (−5.75 vs −3.3 kcal/mol for benzene and −5.47 vs −3.6 kcal/mol for toluene). For the TMA+−benzene complex, Gao et al.,40 using a hybrid quantum mechanics/molecular mechanics representation, have calculated average binding free energies of −0.6 kcal/mol for the contact pair (at 4.7 Å) and of −1.8 kcal/mol for the solvent-separated pair (at 7.5 Å). For the TMA+− benzene complex as well, Duffy et al.,41 using the OPLS (Optimized Potential for Liquid Simulations) force field, have reported a −3.3 kcal/mol average free energy minimum at contact distance (4.75 Å) and no free energy minimum at solvent-separated distance. For the TMA+−phenol complex, Gaberšcě k and Mavri43 have used the solvent reaction field method to estimate a free energy of association of −11.06 kcal/mol. This value is most likely an overestimation, given that an optimization of the gaseous TMA+−phenol complex at the MP2(full)/6-311++G(d,p) level yields ECP = −11.36 kcal/ mol (this work). The solubility of benzene in aqueous quaternary ammonium salt solutions increases with the size of the ion (NH4+ < TMA+ < TEA+).67 This trend across the three ions is the same we observe for the binding free energy of the cation−π complex involving any given aromatic ligand (including benzene). For example, the free energy minimum of the

1), yet are consistent with contact complexes. The complex of NH4+ with 3-methylindole displays a broad free energy well between ∼3 and ∼6 Å. Figure 4a indicates that the two conformers of TEA+ (see Figure 1) possess similar binding affinities for benzene, especially at the minimum of the PMF. Larger ions create weaker electrostatic fields and smaller induced dipoles.20 Yet, as shown in Figure 4e−h, they have stronger affinity for aromatic compounds, despite their charged center being further away. The increase in binding affinity going from NH4+ to TMA+ and TEA+ can thus be attributed to the hydrophobic effect. This is in contrast to the gas-phase binding where electrostatic and polarization forces dominate.2,17,18 The weak shoulders in Figure 4a-d (at about 7−10 Å for TMA+ and TEA+ complexes) and the local free energy minima in Figure 4e-h (at about 6 Å for NH4+ complexes, 7.5 Å for TMA+ complexes, and 9 Å for TEA+ complexes) correspond to the solvent-separated cation−π complexes. Although both TMA+ and TEA+ ions form direct cation−π complexes with comparable separations (∼4.5 Å), solventseparated cation−π complexes have markedly shorter range for TMA+ than for TEA+ (∼7.5 Å for TMA+ and ∼9 Å for TEA+). This is likely due to the fact that TMA+ (the smaller ion) has a more structured hydration shell than TEA+ and creates a more polarized layer of solvent. This is evidenced from the radial distribution function gNO(r) (between the nitrogen of the ion and the oxygen of the water molecules) calculated from 5 ns MD simulations of each ion in 500 H2O molecules at 1 atm and 298.15 K (Figure 5). The function displays two peaks at 4.4 and 7.3 Å for TMA+, and at 5.9 and 8.6 Å for TEA+, corresponding to the first and second solvation shells of the ion. The minimum of the first peak is located at 6.1 for TMA+ and 7.1 Å for TEA+. This difference in thickness of the first solvation shell (about 1.0 Å thicker for TEA+) is consistent with the difference in position of the free energy minimum for the solvent-separated complexes (about 1.5 Å further apart for TEA+). In comparison to the unconstrained and en face constrained simulations, equilibrium separations from the edge-on constrained simulations (Figure 4i-l) are larger due to steric clashes and electrostatic repulsion between the cation and the aromatic hydrogens. The shallow free energy wells (at about 9 2256


Article


Figure 4. Potentials of mean force (PMFs) between the center of the ions (NH4+, TMA+, and TEA+) and the center of the six-membered ring of (a, e, i) benzene, (b, f, j) toluene, (c, g, k) 4-methylphenol, and (d, h, l) 3-methylindole; (a)−(d) are obtained from unconstrained simulations, (e)−(h) from en face constrained simulations, and (i)−(l) from edge-on constrained simulations. The dashed line in (a) represents the PMF obtained with conformer 1b of TEA+. All of the other PMFs involving TEA+ are calculated with conformer 1a only.

ment of their LJ parameters. These models, together with previously reported models for NH4+ interactions,19,32 were used in MD simulations to investigate the binding affinity of quaternary ammonium ions for aromatic compounds in water. Simulations show that the en face binding geometry is favored over the edge-on geometry for all of the complexes in solutionas it is in gas phase. The calculations show that the trend in binding affinity of the three cations with a given ligand is reversed compared to the gas phase, and follows the trend of benzene solubility in aqueous quaternary ammonium salts,67 which suggests a role for cation−π interactions in the salting-in of aromatic compounds. The present parameters extend the application of the Drude force field to the interactions of cations with various aromatic systems. Work is

unconstrained PMF for the benzene complex gets deeper by 0.9 kcal/mol going from NH4+ to TMA+ and by 0.2 kcal/mol going from TMA+ to TEA+. By comparison, the salting-in constant for benzene increases by 0.47 L/mol on going from NH4Br to TMABr and by 0.21 L/mol on going from the TMABr to TEABr.67 This finding suggests that cation−π interactions play a major role in the observed solubility trend: the higher the affinity of benzene with the cation, the higher its solubility in water.

4. SUMMARY Polarizable models for the interactions of TMA+ and TEA+ with water, benzene, toluene, 4-methylphenol, and 3methylindole were optimized through a pair-specific readjust2257


Article


(4) Minoux, H.; Chipot, C. Cation-π Interactions in Proteins: Can Simple Models Provide an Accurate Description? J. Am. Chem. Soc. 1999, 121, 10366−10372. (5) Burley, S. K.; Petsko, G. A. Amino-Aromatic Interactions in Proteins. FEBS Lett. 1986, 203, 139−143. (6) Singh, J.; Thornton, J. M. SIRIUS. An Automated Method for the Analysis of the Preferred Packing Arrangements between Protein Groups. J. Mol. Biol. 1990, 211, 595−615. (7) Mitchell, J. B.; Nandi, C. L.; McDonald, I. K.; Thornton, J. M.; Price, S. L. Amino/aromatic Interactions in Proteins: Is the Evidence Stacked Against Hydrogen Bonding? J. Mol. Biol. 1994, 239, 315− 331. (8) Crowley, P. B.; Golovin, A. Cation-π Interactions in ProteinProtein Interfaces. Proteins: Struct., Funct., Bioinf. 2005, 59, 231−239. (9) Wintjens, R.; Liévin, J.; Rooman, M.; Buisine, E. Contribution of Cation-π Interactions to the Stability of Protein-DNA Complexes. J. Mol. Biol. 2000, 302, 393−408. (10) Gromiha, M. M.; Santhosh, C.; Ahmad, S. Structural Analysis of Cation-π Interactions in DNA Binding Proteins. Int. J. Biol. Macromol. 2004, 34, 203−211. (11) Shi, Z.; Olson, C. A.; Kallenbach, N. R. Cation-π Interaction in Model α-Helical Peptides. J. Am. Chem. Soc. 2002, 124, 3284−3291. (12) Andrew, C. D.; Bhattacharjee, S.; Kokkoni, N.; Hirst, J. D.; Jones, G. R.; Doig, A. J. Stabilizing Interactions between Aromatic and Basic Side Chains in α-Helical Peptides and Proteins. Tyrosine Effects on Helix Circular Dichroism. J. Am. Chem. Soc. 2002, 124, 12706−12714. (13) Tsou, L. K.; Tatko, C. D.; Waters, M. L. Simple Cation-π Interaction between a Phenyl Ring and a Protonated Amine Stabilizes an α-Helix in Water. J. Am. Chem. Soc. 2002, 124, 14917−14921. (14) Prajapati, R. S.; Sirajuddin, M.; Durani, V.; Sreeramulu, S.; Varadarajan, R. Contribution of Cation-π Interactions to Protein Stability. Biochemistry 2006, 45, 15000−15010. (15) Lummis, S. C. R.; Beene, D. L.; Harrison, N. J.; Lester, H. A.; Dougherty, D. A. A Cation-π Binding Interaction with a Tyrosine in the Binding Site of the GABAC Receptor. Chem. Biol. 2005, 12, 993− 997. (16) Lehn, J.-M.; Meric, R.; Vigneron, J.-P.; Cesario, M.; Guilhem, J.; Pascard, C.; Asfari, Z.; Vicens, J. Binding of Acetylcholine and Other Quaternary Ammonium Cations by Sulfonated Calixarenes. Crystal Structure of a [choline-tetrasulfonated calix[4]arene] Complex. Supramol. Chem. 1995, 5, 97−103. (17) Caldwell, J. W.; Kollman, P. A. Cation-π Interactions: Nonadditive Effects Are Critical in Their Accurate Representation. J. Am. Chem. Soc. 1995, 117, 4177−4178. (18) Archambault, F.; Chipot, C.; Soteras, I.; Luque, F. J.; Schulten, K.; Dehez, F. Polarizable Intermolecular Potentials for Water and Benzene Interacting with Halide and Metal Ions. J. Chem. Theory Comput. 2009, 5, 3022−3031. (19) Orabi, E. A.; Lamoureux, G. Cation-π and π-π Interactions in Aqueous Solution Studied Using Polarizable Potential Models. J. Chem. Theory Comput. 2012, 8, 182−193. (20) Lamoureux, G.; Orabi, E. A. Molecular Modelling of Cation-π Interactions. Mol. Simul. 2012, 38, 704−722. (21) Mavri, J.; Koller, J.; Hadz̆i, D. Ab Initio and AM1 Calculations on Model Systems of Acetylcholine Binding: Complexes of Tetramethylammonium with Aromatics, Neutral and Ionic Formic Acid. J. Mol. Struct.: THEOCHEM 1993, 283, 305−312. (22) Kim, K. S.; Lee, J. Y.; Lee, S. J.; Ha, T.-K.; Kim, D. H. On Binding Forces between Aromatic Ring and Quaternary Ammonium Compound. J. Am. Chem. Soc. 1994, 116, 7399−7400. (23) Lee, J. Y.; Lee, S. J.; Choi, H. S.; Cho, S. J.; Kim, K. S.; Ha, T.-K. Ab Initio Study of the Complexation of Benzene with Ammonium Cations. Chem. Phys. Lett. 1995, 232, 67−71. (24) Pullman, A.; Berthier, G.; Savinelli, R. Theoretical Study of Binding of Tetramethylammonium Ion with Aromatics. J. Comput. Chem. 1997, 18, 2012−2022.

Figure 5. N−O radial distribution functions for TMA+ (black) and TEA+ (red) in water at 298.15 K, obtained from MD simulations of the ions in aqueous solution.

in progress to generalize the parameters and allow for their use in proteins, protein−lipid complexes, and protein−ligand complexes.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.7b11983. Coordinates of the optimized complexes, positions, and depths of the PMF minima from Figure 4, and force field parameters (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: (514) 8482424 ext. 5314. ORCID

Guillaume Lamoureux: 0000-0001-7656-0176 Notes

The authors declare no competing financial interest. † On leave from Department of Chemistry, Faculty of Science, Assiut University, Assiut 71516, Egypt (E.A.O.).

■

ACKNOWLEDGMENTS The authors thank Alexander D. MacKerell, Jr., for reading the manuscript and for his valuable comments. This work was supported in part by an FQRNT Nouveaux chercheurs grant and an NSERC Discovery grant to G.L., and a Garnet Strong Scholarship and a Power Corporation of Canada Graduate Fellowship to E.A.O. Computational resources were provided by Calcul Québec and Compute Canada.

■

REFERENCES

(1) Stauffer, D. A.; Barrans, R. E.; Dougherty, D. A. Concerning the Thermodynamics of Molecular Recognition in Aqueous and Organic Media. Evidence for Significant Heat Capacity Effects. J. Org. Chem. 1990, 55, 2762−2767. (2) Dougherty, D. A. Cation-π Interactions in Chemistry and Biology: A New View of Benzene, Phe, Tyr, and Trp. Science 1996, 271, 163−168. (3) Gallivan, J. P.; Dougherty, D. A. Cation-π Interactions in Structural Biology. Proc. Natl. Acad. Sci. U.S.A. 1999, 96, 9459−9464. 2258


Article

The Journal of Physical Chemistry B (25) Zhu, W.-L.; Tan, X.-J.; Puah, C. M.; Gu, J.-D.; Jiang, H.-L.; Chen, K.; Felder, C. E.; Silman, I.; Sussman, J. L. How Does Ammonium Interact with Aromatic Groups? A Density Functional Theory (DFT/B3LYP) Investigation. J. Phys. Chem. A 2000, 104, 9573−9580. (26) Felder, C. E.; Jiang, H.-L.; Zhu, W.-L.; Chen, K.; Silman, I.; Botti, S. A.; Sussman, J. L. Quantum/Classical Mechanical Comparison of Cation-π Interactions between Tetramethylammonium and Benzene. J. Phys. Chem. A 2001, 105, 1326−1333. (27) Pullman, A.; Berthier, G.; Savinelli, R. Components of the Interaction Energy of Benzene with Na+ and Methylammonium Cations. J. Mol. Struct.: THEOCHEM 2001, 537, 163−172. (28) Kim, D.; Hu, S.; Tarakeshwar, P.; Kim, K. S.; Lisy, J. M. Cation-π Interactions: A Theoretical Investigation of the Interaction of Metallic and Organic Cations with Alkenes, Arenes, and Heteroarenes. J. Phys. Chem. A 2003, 107, 1228−1238. (29) Reddy, A. S.; Sastry, G. N. Cation [M = H+, Li+, Na+, K+, Ca2+, Mg2+, NH4+, and NMe4+] Interactions with the Aromatic Motifs of Naturally Occurring Amino Acids: A Theoretical Study. J. Phys. Chem. A 2005, 109, 8893−8903. (30) Marshall, M. S.; Steele, R. P.; Thanthiriwatte, K. S.; Sherrill, C. D. Potential Energy Curves for Cation-π Interactions: Off-Axis Configurations Are Also Attractive. J. Phys. Chem. A 2009, 113, 13628−13632. (31) Singh, N. J.; Min, S. K.; Kim, D. Y.; Kim, K. S. Comprehensive Energy Analysis for Various Types of π-Interaction. J. Chem. Theory Comput. 2009, 5, 515−529. (32) Wang, S.; Orabi, E. A.; Baday, S.; Bernèche, S.; Lamoureux, G. Ammonium Transporters Achieve Charge Transfer by Fragmenting Their Substrate. J. Am. Chem. Soc. 2012, 134, 10419−10427. (33) Baday, S.; Orabi, E. A.; Wang, S.; Lamoureux, G.; Bernèche, S. Mechanism of NH4+ Recruitment and NH3 Transport in Rh Proteins. Structure 2015, 23, 1550−1557. (34) Couture, J.-F.; Hauk, G.; Thompson, M. J.; Blackburn, G. M.; Trievel, R. C. Catalytic Roles for Carbon-Oxygen Hydrogen Bonding in SET Domain Lysine Methyltransferases. J. Biol. Chem. 2006, 281, 19280−19287. (35) Daze, K. D.; Hof, F. The Cation-π Interaction at ProteinProtein Interaction Interfaces: Developing and Learning from Synthetic Mimics of Proteins That Bind Methylated Lysines. Acc. Chem. Res. 2013, 46, 937−945. (36) Armstrong, C. M. Interaction of Tetraethylammonium Ion Derivatives with the Potassium Channels of Giant Axons. J. Gen. Physiol. 1971, 58, 413−437. (37) French, R. J.; Shoukimas, J. J. Blockage of Squid Axon Potassium Conductance by Internal tetra-N-alkylammonium Ions of Various Sizes. Biophys. J. 1981, 34, 271−291. (38) Deakyne, C. A.; Meot-Ner, M. Unconventional Ionic Hydrogen Bonds. 2. NH+···π Complexes of Onium Ions with Olefins and Benzene Derivatives. J. Am. Chem. Soc. 1985, 107, 474− 479. (39) Xu, Y.; Shen, J.; Zhu, W.; Luo, X.; Chen, K.; Jiang, H. Influence of the Water Molecule on Cation-π Interaction: Ab Initio Second Order Møller-Plesset Perturbation Theory (MP2) Calculations. J. Phys. Chem. B 2005, 109, 5945−5949. (40) Gao, J.; Chou, L. W.; Auerbach, A. The Nature of Cation-π Binding: Interactions between Tetramethylammonium Ion and Benzene in Aqueous Solution. Biophys. J. 1993, 65, 43−47. (41) Duffy, E. M.; Kowalczyk, P. J.; Jorgensen, W. L. Do Denaturants Interact with Aromatic Hydrocarbons in Water? J. Am. Chem. Soc. 1993, 115, 9271−9275. (42) Chipot, C.; Maigret, B.; Pearlman, D. A.; Kollman, P. A. Molecular Dynamics Potential of Mean Force Calculations: A Study of the Toluene-Ammonium π-Cation Interactions. J. Am. Chem. Soc. 1996, 118, 2998−3005. (43) Gaberšcě k, M.; Mavri, J. Phenol Forms Complexes with Tetramethylammonium Ions in Aqueous Solution? Chem. Phys. Lett. 1999, 308, 421−427.

(44) Gallivan, J. P.; Dougherty, D. A. A Computational Study of Cation-π Interactions vs Salt Bridges in Aqueous Media: Implications for Protein Engineering. J. Am. Chem. Soc. 2000, 122, 870−874. (45) Sa, R.; Zhu, W.; Shen, J.; Gong, Z.; Cheng, J.; Chen, K.; Jiang, H. How Does Ammonium Dynamically Interact with Benzene in Aqueous Media? A First Principle Study Using the Car-Parrinello Molecular Dynamics Method. J. Phys. Chem. B 2006, 110, 5094− 5098. (46) Lopes, P. E. M.; Huang, J.; Shim, J.; Luo, Y.; Li, H.; Roux, B.; MacKerell, A. D., Jr. Polarizable Force Field for Peptides and Proteins Based on the Classical Drude Oscillator. J. Chem. Theory Comput. 2013, 9, 5430−5449. (47) Lemkul, J. A.; Huang, J.; Roux, B.; MacKerell, A. D., Jr. An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications. Chem. Rev. 2016, 116, 4983−5013. (48) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian 09, revision B.01; Gaussian, Inc.: Wallingford, CT, 2009. (49) Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553−566. (50) Brooks, B. R.; Brooks, C. L.; Mackerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al. CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30, 1545−1614. (51) Lamoureux, G.; Roux, B. Modeling Induced Polarization with Classical Drude Oscillators: Theory and Molecular Dynamics Simulation Algorithm. J. Chem. Phys. 2003, 119, 3025−3039. (52) Lamoureux, G.; Harder, E.; Vorobyov, I. V.; Roux, B.; MacKerell, A. D., Jr. A Polarizable Model of Water for Molecular Dynamics Simulations of Biomolecules. Chem. Phys. Lett. 2006, 418, 245−249. (53) Lopes, P. E. M.; Lamoureux, G.; Roux, B.; MacKerell, A. D., Jr. Polarizable Empirical Force Field for Aromatic Compounds Based on the Classical Drude Oscillator. J. Phys. Chem. B 2007, 111, 2873− 2885. (54) Lopes, P. E. M.; Lamoureux, G.; MacKerell, A. D., Jr. Polarizable Empirical Force Field for Nitrogen-Containing Heteroaromatic Compounds Based on the Classical Drude Oscillator. J. Comput. Chem. 2009, 30, 1821−1838. (55) P. E. M. Lopes and A. D. MacKerell, Jr., personal communication. (56) Orabi, E. A.; Lamoureux, G. Polarizable Interaction Model for Liquid, Supercritical, and Aqueous Ammonia. J. Chem. Theory Comput. 2013, 9, 2035−2051. (57) Orabi, E. A.; Lamoureux, G. Molecular Dynamics Investigation of Alkali Metal Ions in Liquid and Aqueous Ammonia. J. Chem. Theory Comput. 2013, 9, 2324−2338. (58) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Explicit Reversible Integrators for Extended Systems Dynamics. Mol. Phys. 1996, 87, 1117−1157. (59) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (60) Lagüe, P.; Pastor, R. W.; Brooks, B. R. Pressure-Based LongRange Correction for Lennard-Jones Interactions in Molecular Dynamics Simulations: Application to Alkanes and Interfaces. J. Phys. Chem. B 2004, 108, 363−368. (61) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695−1697. (62) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant Pressure Molecular Dynamics Algorithms. J. Chem. Phys. 1994, 101, 4177− 4189. (63) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011−1021. 2259


Article

The Journal of Physical Chemistry B (64) Souaille, M.; Roux, B. Extension to the Weighted Histogram Analysis Method: Combining Umbrella Sampling with Free Energy Calculations. Comput. Phys. Commun. 2001, 135, 40−57. (65) Luzhkov, V. B.; Åqvist, J. Mechanisms of Tetraethylammonium Ion Block in the KcsA Potassium Channel. FEBS Lett. 2001, 495, 191−196. (66) Meot-Ner, M.; Deakyne, C. A. Unconventional Ionic Hydrogen Bonds. 1. CHδ+···X Complexes of Quaternary Ions with n- and π-Donors. J. Am. Chem. Soc. 1985, 107, 469−474. (67) Desnoyers, J. E.; Pelletier, G. E.; Jolicoeur, C. Salting-In by Quaternary Ammonium Salts. Can. J. Chem. 1965, 43, 3232−3237.

2260


CationâÏ Interactions between Quaternary Ammonium Ions and

Recommend Documents

CationâÏ Interactions between Quaternary Ammonium Ions and