Understanding Water Structure in an Ion-Pair Solvation Shell in the

Apr 22, 2019 - Krembil Research Institute, University Health Network , Toronto ... Department of Chemistry, University of Toronto, Toronto , Ontario M...
0 downloads 0 Views 2MB Size
Article Cite This: J. Phys. Chem. B XXXX, XXX, XXX−XXX

pubs.acs.org/JPCB

Understanding Water Structure in an Ion-Pair Solvation Shell in the Vicinity of a Water/Membrane Interface Myong In Oh,† Mayuri Gupta,† and Donald F. Weaver*,†,‡,§,∥ †

Krembil Research Institute, University Health Network, Toronto, Ontario M5T 0S8, Canada Department of Chemistry, University of Toronto, Toronto, Ontario M5S 3H6, Canada § Department of Medicine, University of Toronto, Toronto, Ontario M5G 2C4, Canada ∥ Department of Pharmaceutical Sciences, University of Toronto, Toronto, Ontario M5S 3M2, Canada ‡

J. Phys. Chem. B Downloaded from pubs.acs.org by RICE UNIV on 05/05/19. For personal use only.

S Supporting Information *

ABSTRACT: The anomalous properties of interfacial water at the surface of a lipid membrane and their implications on nearby chemical processes are well recognized. However, we have found that ion pairing thermodynamics may not be significantly affected by interfacial water in a classical, nonpolarizable force field. To trace the root cause of such a counterintuitive finding, we performed atomistic molecular dynamics simulations to explore the impact of polarizable interactions and characterize the hydration structure of a sodium chloride (NaCl) ion pair at the surface of a model lipid membrane and in a bulk phase. Our study reveals that the effect of the aqueous interface on the first solvation shell of the ion pair and thus on the ion pairing thermodynamics becomes more pronounced in the polarizable model, and that the free energy profile along the interionic distance cannot capture the difference in the degree of solvent participation in ion pairing at the water/membrane interface. This study also forms the basis for the future design of a reaction coordinate that takes the behavior of the interfacial water into account.



INTRODUCTION The lipid bilayer is a defining structural motif in living cells, found in both cellular compartments and as the delineating membrane. It is composed of amphiphilic phospholipids that self-assemble into two leaflets in the presence of water, with their hydrophilic heads exposed to the aqueous medium and fatty acid tails internalized to create a hydrophobic interior. Interactions between phospholipids and nearby water molecules create a unique biomolecular environment in the vicinity of the water/membrane interface. For example, the interfacial water is characterized by lower dielectric constants,1,2 lower translational and rotational entropies,3 higher lateral proton diffusivity,4,5 slower relaxation rates,3,6,7 and a mosaic of water orientation structures,8,9 as compared to bulk water. A comprehensive review on the recent experimental and theoretical advances in understanding the properties of various aqueous interfaces is well presented in ref 10. Recently, of particular importance is an understanding of membrane hydration in medicine. As membrane lipid composition differs greatly among cells of different tissue types and cell stages,11−14 diffusion magnetic resonance imaging (MRI) uses information on water diffusion near cellular membranes to map the circuitry of the human brain and to detect different stages of stroke, cancers, and metastasis.15−18 For example, malignant tissues are usually © XXXX American Chemical Society

marked by the reduction in water diffusion due to the underlying cell proliferation in the tumor.16 The anomalous properties of the interfacial water can, in turn, affect various chemical processes that take place at or near the membrane surface, such as membrane-assisted amyloid formation,19 selective ion binding by membrane proteins,20 ion-membrane interactions,21−25 and ion-pairing facilitated transmembrane transport of drugs and metals.26−30 Despite its practical significance, however, little or no theoretical work has been done to understand, at the atomistic level, the effect of highly organized water structure on solvation and on the complexation behavior of solutes at the membrane surface. In light of this requirement, we focus on simple ion pairing between sodium (Na+) and chloride (Cl−) ions, as ion pairing is the most fundamental and significant interaction in life. Although water-mediated ion−ion interactions in bulk water have been extensively explored,31−35 those that take place at the aqueous interface of a biological membrane still remain incompletely understood at the atomistic scale. We perform atomistic molecular dynamics (MD) simulations and compute Received: February 10, 2019 Revised: April 22, 2019 Published: April 22, 2019 A

DOI: 10.1021/acs.jpcb.9b01331 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B one-dimensional free energy profiles along the interionic distance for the simple ionic interaction and compare them when the ion pair is located near or far from the water/ membrane interface. We then characterize the change in hydration structure of the ion pair during ion association and dissociation that occur at different solvation environments, in order to elucidate how ion−ion interactions are affected by the presence of a water/membrane interface and suggest how reaction coordinates should be chosen for the interfacial ion pairing. We also address how the explicit inclusion of electronic polarization would influence ion pairing thermodynamics at the interface vs in bulk solution.



COMPUTATIONAL METHODS A 20 ns-long isobaric−isothermal (NPT) ensemble of the system was generated under 300 K and 1 atm by using the NAMD package version 2.13.36 We used a model membrane composed solely of 60 phosphatidylcholine (POPC, 1palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) molecules to reduce a computational burden37−39 and compare it with other theoretical works.6−8 The membrane was then solvated in an aqueous solution at infinite ion dilution that contains 7200 H2O molecules and a single pair of Na+ and Cl− ions. The CHARMM general force field (CGenFF)40−43 was employed to model the phospholipid molecules and the ions, and the solvent was explicitly represented using the modified TIP3P model.44 Every covalent bond involving a hydrogen atom was constrained and kept rigid using SETTLE45 for water and SHAKE46 for phospholipids. Long-range electrostatic interactions were treated using the particle mesh Ewald (PME) algorithm,47 while van der Waals interactions (modeled by the Lennard-Jones potential) were truncated at a cutoff distance of 12 Å. Langevin dynamics with a damping coefficient of 1 ps−1 was employed to thermostat the system. A Langevin piston Nosé−Hoover method48,49 was used to barostat the system, in which the pressure coupling is isotropic only in the x- and ydirections. All the systems were first equilibrated for 20 ns in total, and then the production runs of 20 ns were generated. The area per lipid was 67−70 Å2, and the membrane thickness was 54−59 Å, which were in good agreement with other works.8,50 MD trajectories were generated by using the velocity-Verlet algorithm51 with a time step of 2.0 fs. VMD version 1.9.352 was used for visualization. Scripts written inhouse were used to analyze our results. Utilizing Drude Prepper53 in CHARMM-GUI, we also prepared the same systems but with the Drude polarizable force field54 and the SWM4-NDP55 water model to assess the impact of polarizable interactions on ion pairing thermodynamics and solvation. In this study, 5 ns-long MD trajectories were generated with a time step of 1 fs after equilibration. The aqueous environment was divided into three distinct regions, based on the perpendicular distance (|z|) between the ion pair and the center of mass (COM) of the bilayer (see Figure S1). The ion pair was said to be near or at the water/ membrane interface when |z| ≤ 19 Å (green region in Figure 1a), in the intermediate region when 19 Å ≤ |z| ≤ 25 Å (red region), and in a bulk-like environment when |z| ≥ 25 Å (blue region). The stratification of the aqueous environment by distance was defined according to the electron density profile of water acquired from the Density Profile Tool.56 Umbrella sampling57 and the weighted histogram analysis method (WHAM)58 were performed to generate the potential of mean force (PMF) profiles for ion pairing thermodynamics

Figure 1. (a) Division of the aqueous environment into three distinct regions based on the orthogonal distance from the COM of the lipid bilayer: green for water near the membrane surface, blue for bulk-like water, and red for water in between. The POPC molecules are colored in gray, and Na+ and Cl− ions are represented by blue and yellow spheres, respectively. (b) PMF profiles along the interionic distance (r Na+ ··· Cl−) for ion association in the three different regions of solvation, created with the nonpolarizable force field. Na+ and Cl− ions are represented by blue and green spheres, respectively. The oxygen and hydrogen atoms of a H2O molecule are colored in red and white, respectively. The first minimum of the PMF curve was set to zero.

and to characterize the structure of the first hydration shell of each ion as the interionic distance varies (hereafter, we write r = X Å to indicate that the interionic distance r Na+ ··· Cl− is constrained to X Å). A harmonic potential (with a force constant of 2.0 kcal/mol·Å2) was applied to constrain (1) the separation distance r between the cation and the anion for umbrella sampling and (2) the perpendicular distance |z| between the ion pair and the COM of the POPC bilayer. We generated nine different windows, spaced 0.5 Å apart, in each of which r is centered around a certain value ranging from 2 to 6 Å. The WHAM58 was then used to generate the onedimensional unbiased PMF curves for the ion−ion interaction. We also carried out additional MD simulations (with the nonpolarizable force field) of graphene solvated in water at infinite ion dilution for the purpose of comparison. The system contains a single layer of graphene (5 nm × 5 nm in size) fixed in location, 9316 H2O molecules, and a single pair of Na+ and Cl− ions. The same constraints were applied to manipulate (1) the interionic distance and (2) the perpendicular distance between the ion pair and the COM of the graphene layer. The ions were allowed to translate only on a plane parallel to the graphene (as in ref 56), and we chose the planes of |z| = 5 Å and |z| = 40 Å as near and bulk-like regions of solvation, respectively.



RESULTS AND DISCUSSION Ion Pairing Thermodynamics. Figure 1b shows three one-dimensional free energy profiles along the separation distance between the cation and the anion (r Na+ ··· Cl− or simply r), each of which differs only by the solvation zone of the ion pair. As typical of solvent-mediated interactions of ions of opposite charge,59 two energy minima appear at r = 2.7 Å and r = 5.1 Å, respectively, which are separated by an energy barrier of 2.7−2.8 kcal/mol (for the transition from CIP to SSIP) B

DOI: 10.1021/acs.jpcb.9b01331 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B located at r = 3.5 Å; the first minimum represents the state of a contact ion pair (CIP), where the ions are in direct contact with each other, and the second one corresponds to the state of a solvent-separated ion pair (SSIP), where one or more water molecules shared by both ions mediate the ion−ion interaction. Interestingly, the PMFs obtained from our simulations with the classical nonpolarizable force field show no or only a very small difference in ion pairing thermodynamics at the membrane surface vs in the bulk phase (i.e., the difference is less than 2.0% of the peak height). This is in contrast to a simple continuum dielectric description, which would suggest that the lower dielectric permittivity of water near the surface weakens screening of electrostatic interactions between oppositely charged ions, possibly leading to a stronger stability of the ion pair near the interface (an example of the dielectric effect on the stability of ion pairs can be found in the molecular dynamics studies on the PMF of a NaCl pair in supercritical water60). Since a lipid membrane consists mainly of neutral, zwitterionic, and anionic phospholipids, it carries an effective charge density varying from −0.002 to −0.3 C/m2,61 triggering charge-induced alignment of water molecules at the interface.62,63 Previous computational and experimental studies predict that, regardless of the degree of hydrophilicity of an aqueous interface, the dielectric constant of confined and interfacial water can be reduced by approximately an order of magnitude,1,2,64−66 as it forms layered structures and thus shows suppressed polarizability.2,10 We also demonstrated suppressed rotational motion of the water dipoles at the interface by quantifying their orientational relaxation times (Section S2 in the Supporting Information). Our finding is also in contrast to the findings that the watermediated ion−ion attractions at the water/vapor interface,67 at the water/1,2-dichloroethane interface,68 and at the ice/vapor interface69 are much stronger and longer ranged than that in bulk water. At the liquid/vapor interface, the enhanced ion− ion attractions are driven by the reduction of capillary deformations of the liquid surface in addition to the simple dielectric effect discussed above.67 The enhanced ion clustering at the ice/vapor interface originates from the minimization of liquidity of the ice surface created by ions.69 The discrepancy may stem from the ignorance of electronic polarization in our calculations. Recently, Dang et al.70 performed a detailed study of ion pairing thermodynamics and kinetics for Na+ and Cl− ions in the bulk and near the water liquid/vapor interface, using both classical (TIP4PEw71) and polarizable (Dang-Chang72) water models. It is appealing that they found that interfacial ion dissociation is significantly slower, as consistent with experiments, only when the effect of polarizability is explicitly included in their potential models.70 Motivated by the work of Dang et al.,70 we included the effect of polarizable interactions in our simulations by rebuilding the entire system with the Drude polarizable force field54 and the SWM4-NDP55 water model. The Drude oscillator model73 introduces a massless, charge-carrying virtual particle attached to each polarizable atom via a harmonic oscillator. Hence, electronic polarization is represented by the displacement of the auxiliary particle under the influence of an external electric field. As plainly shown in Figure 2, the CIP state of the interfacial ion pair is more stable than those of the noninterfacial ion pairs as manifested by the height of the energy barrier separating the CIP and SSIP states

Figure 2. PMF profiles along the interionic distance (r Na+ ··· Cl−) for ion association in the three different regions of solvation, created with the polarizable force field. The inset shows the stabilization of the CIP state due to the presence of phospholipids. The same color scheme was used as in Figure 1.

(i.e., ΔWCIP→SSIP = 1.75 kcal/mol for near vs 1.38 kcal/mol for intermediate and 1.28 kcal/mol for bulk). Although the use of the polarizable model significantly reduces the peak heights, it is consistent with the outcome from the work of Dang et al.70 that the significant degree of stabilization of the CIP state due to the interface (see the inset embedded in Figure 2) is only detected with the inclusion of the electronic polarizability in the model. The small difference in PMFs in Figure 1b may be ascribed as well to the incompleteness of the reaction coordinate to map the reaction pathway. Water influences the thermodynamics and kinetics of chemical processes by mediating interactions between solutes. However, a critical issue arises directly from the challenge in the selection of proper reaction coordinates that completely describe the mechanism of even a simple interaction involving, for instance, only two monovalent ions of opposite charge. In particular, the interionic distance alone, albeit the most common reaction coordinate, does not suffice to delineate the mechanism of ion pairing in aqueous solution, demanding a set of reaction coordinates that account for both solvent and solute degrees of freedom during the process.34,70,74−78 By way of illustration, Roy and his coworkers76−78 selected the ion−solvent coordination number (or the number of water molecules in the first hydration shell of an ion) as the reaction coordinate for the structural and dynamical characterization of ion solvation and, in combination with the ion pair distance, the development of the Marcus theory of ion pairing. The necessity of such a reaction coordinate becomes more pronounced when ions are located at aqueous interfaces where they show unusual solvation behavior.5,79−85 Therefore, it is essential to understand with full atomistic resolution how the hydration structure of an ion pair differs as it undergoes complexation near a membrane surface in comparison to ion pairing in the bulk aqueous phase. In the following section, we conduct a structural analysis of the first hydration shell of the ion pair with and without the inclusion of the polarization effect. Ion-Pair Hydration Structure without the Inclusion of Polarizable Interactions. Water Coordination Number and Angle. First, we determined the average water coordination number (⟨CN⟩) of each ion and how it is altered as the ion pair dissociates in different solvation regions by using the following expression C

DOI: 10.1021/acs.jpcb.9b01331 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B CN

=

∑ w ∈ FSS

the dipole that the two ions form (the ion pair can be treated as a dipole since the interionic distance is constrained to a defined length). Figure S3 shows the normalized distribution of the dipole orientation ϕ with respect to the unit normal vector of the membrane surface (i.e., z-axis). It is clearly shown that the Na+ ion dwells in the more inner region than the Cl− ion when the dipole interacts with POPC in the vicinity of the membrane surface. This is because the dipole aligns along the charged groups of the phospholipid molecule. On the contrary, the dipole orientation is equally distributed in every direction in bulk solution. Next, we examined the shape of the distribution of CN for each ion. We defined γ = AUCleft/AUCright (where AUC stands for the area under the curve) as the measure of distribution skewness, which is simply the ratio of the area under the curve from CN = 0 to the peak (AUCleft) to that from the peak to CN → ∞ (AUCright). If γ = 1, then the distribution is symmetric; otherwise, it shows a tendency to skew to the left (higher population of lower CN states) when γ > 1 or to the right (higher population of higher CN states) when γ < 1 (see Figure S4 for an example). The analysis of γ discloses two interesting trends, as shown in Figure 3b,d. First, when the ion pair moves toward the water/membrane interface, the distribution becomes less symmetric. Specifically, the Na+ ion is characterized by highly deviant γ values (higher than 1.5) when it settles at the interface (green). Second, regardless of the location of ion solvation, the distribution is much less symmetric for the Na+ ion as compared to its counterion, whose γ values are scattered around the horizontal purple dotted line (γ = 1). To examine specific orientational preferences of the water dipoles in the first solvation shell of the ion pair, we computed the normalized distribution of the angle (θ) between the dipole moment of an adjacent water molecule and the displacement vector that connects the central ion and the oxygen site of the water molecule (refer to the inset in Figure 4a). As an example, Figure 4a shows how water molecules are

Θ(g − rion − O) (1)

where Θ is the Heaviside step function, g is the location of the first local minimum in the radial distribution function of the water oxygen atoms around the ions, and rion−O is the distance between the central ion and the oxygen site. The summation is taken over all the water molecules (w) in the first solvation shell (FSS), and ⟨···⟩ indicates that the quantity is averaged over the trajectories. As shown in Figure 3a,c, ⟨CN⟩ increases abruptly at r = 3.5 Å, as consistent with the PMF, which corresponds to the

Figure 3. (a, c) Average water CN and (b, d) skewness parameter γ for the distribution of CN of Na+ and Cl− ions as their separation distance (r Na+ ··· Cl−) increases in the bulk-like (blue circles), intermediate (red squares), and interfacial (green stars) regions of solvation. The purple dotted lines indicate where γ = 1.

energy hill between CIP and SSIP. The change in ⟨CN⟩, Δ⟨CN⟩ = 1.2−1.8, occurs as one or two water molecules are engaged in the formation of water bridges between the cation and the anion (see also Figure 5 and the related discussion). Before and after this transitional slope, no significant disruption in the degree of hydration is observed, as indicated by the plateaus. A Cl− ion is associated with a higher ⟨CN⟩ because it is larger in size. While in the bulk-like environment, ⟨CN⟩ = 5.88 ± 0.03 and ⟨CN⟩ = 7.20 ± 0.04 for the Na+ and Cl− ions, respectively, and these values are in good agreement with diffraction data.86 We also note that the degree of ionic hydration decreases as the ion pair is located closer to the membrane surface due to phospholipid−ion interactions. In particular, the reduction in hydration is more pronounced for the Na+ ion (that is, CN Na+,bulk − CN Na+,near = 0.72 ± 0.23 vs CN Cl−,bulk − CN Cl−,near = 0.58 ± 0.10). This is attributed to the preferential interactions of the ions with specific segments of a POPC molecule. In other words, the hydrophilic headgroup of POPC consists of a positively charged choline group (located at the tip) and a negatively charged phosphate group (linking the choline group to the glycerol group and the hydrophobic tail), and hence, the Na+ ion is much less exposed to the solvent, since it is bound to the phosphate group. This becomes evident by analyzing the orientational preference of

Figure 4. (a) Normalized histograms of the angle (θ) between a water dipole and a Na+ ion (blue) or a Cl− ion (green) when the ion pair is separated by 6.0 Å in bulk solution. The inset shows how θ was defined. (b, c) Probability (P) of finding the most probable angle (θmax) for Cl− (left) and Na+ (right) as the interionic distance (r Na+ ··· Cl−) differs. The same color scheme was used as in Figure 1. D

DOI: 10.1021/acs.jpcb.9b01331 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B oriented around the Na+ ion (blue) and the Cl− ion (green) when the ion pair is in the dissociated state (i.e., r = 6.0 Å) and located in the bulk-like environment. This is in agreement with other computational studies.87,88 The maximum probability is found when the vectors are aligned in the same direction for the Na+ ion and when θ ≈ 130.7◦ for the Cl− ion. These two angles will be denoted as θmax hereafter. Though not shown, the same propensity is observed for the ion pair residing in the intermediate region and near the membrane surface; namely, the distribution shape and the value of θmax remain the same, irrespective of how far the ion pair is from the membrane. We also investigated how the frequency of the occurrence of θmax would differ as the ions are coupled in each solvation zone. To this end, we compared the values of P(cos θmax) for each ion as r Na+ ··· Cl− varies, while the ions reside in the same region of solvation, thereby creating Figure 4b,c for Cl− and Na+, respectively. The same curve pattern is found as in Figure 3a,c, including the existence of two plateaus separated by a sharp change of P at r = 3.5 Å. We found that although very small (the difference is less than 0.02), P(cos θmax) of Cl− becomes slightly higher when nearing the membrane surface (green curve in Figure 4b); whereas P(cos θmax) of Na+ is less sensitive to the solvation location (Figure 4c). Fractional Number of Water Bridges. We characterized the fractional number of water bridges (B), which are simultaneously bonding to both the cation and the anion, by measuring the same quantity as defined by Mullen et al.35 B=



min(f+ , f− )

Figure 5. (a) Typical snapshots showing the ion pair involving Na+ and Cl− ions and different numbers (B) of water bridges. (b−d) Probability (P) of finding the solvent state of B = 0, B = 1, and B = 2, respectively, as the interionic length (r Na+ ··· Cl−) changes. The same color scheme was used as in Figure 1.

the bulk-like water, red for the water in the intermediate region, and green for the interfacial water). First, the curves in each B-state are identical in shape regardless of the solvation location. Though relatively small (i.e., ΔP(B) ≤ 0.1), marginal enhancement in P(B) is detected when the ions have interactions with POPC molecules, implying some important role of phospholipids in the formation of interionic water bridges. Second, there exists a unique transition point r* (we omitted the subscript Na+···Cl− for convenience), after which (1) the order of the curves in terms of magnitude of P(B) is completely reversed, and (2) the difference in P(B) becomes less (for B = 0) or more noticeable (for B = 1 and B = 2); that is, r* ≈ 4.2 Å for B = 0, r* ≈ 5.1 Å for B = 1, r* ≈ 4.0 Å for B = 2, and r* ≈ 3.8 Å for B = 3 (not shown). The decrease in r* as the B-coordinate increases is attributable to the fact that the water bridges act as glue to “tighten” the ionic interaction. Third, irrespective of the solvation location, each B state is dominant in a specific range of r. This is more clearly seen in Figure 6, which shows the probability of finding a specific Bstate of the NaCl ion pair separated by r in the bulk region (the figures corresponding to the intermediate and near regions of solvation are not presented as they exhibit exactly the same shape). The ion pair with no water bridges (B = 0) is more probable than the other B-states when r ≤ 3.5 Å and r ≥ 5.5 Å; the ion pair with a single water bridge (B = 1) is the most prevalent when 4.5 Å ≤ r ≤ 5.5 Å; the state of the ion pair with a double water bridge (B = 2) is more populated when r ≈ 4.0 Å. This third trend is partly manifested in the PMF profiles (Figure 1b). For instance, B = 0 is the most favorable state when r ≤ 3.5 Å, which reflects the CIP state of the ion pair and when r ≥ 5.5 Å, which reflects its dissociated state. Likewise, B = 1 is the most probable state when 4.5 Å ≤ r ≤ 5.5 Å, as the system is in the SSIP state. The analysis of the B-coordinate raises a critical question regarding explicit inclusion of phospholipid headgroups in the set of reaction coordinates to describe the interfacial ion pairing process and its hydration structure. To gain an insight, we replaced the POPC membrane with graphene (as a simple representation of a hydrophobic analogue without any polar headgroups) and repeated the same procedure. Figure 7a

(2)

w ∈ FSS

where the summation is taken over all the water molecules (w) in the FSS, and f∈[0,1] is the continuous indicator function given by f± =

1 − tanh[a(r± − b±)] (3)

2 +



where r± is the distance between Na (or Cl ) and the oxygen atom (or the closest hydrogen atom), and b± is the location of the first local minimum in the radial distribution function of Na+ (b+ = 3.23 Å) or Cl− (b− = 2.88 Å). The parameter a was set to 3, as in the study of Mullen et al.35 Figure 5a contains typical snapshots of the ion pair involving different numbers of water bridges mediating the ion−ion interaction. We studied B since it was found to be an important component of solvent coordinates for ion pairing.34,35 Yonetani34 computed the solvent-coordinate free-energy landscapes and the dissociation rate constants for four different ion pairs (M+···Cl−, where M+ = Li+, Na+, K+, and Cs+) using MD simulations and the adaptive biasing force (ABF) algorithm. He found that the formation of water bridges is a key step in ion-pair dissociation, as it lowers the dissociation energy barrier.34 We examined three different solvent-bridge states, namely, B = 0, 1, and 2, because the states of B = 3 and higher are of great rarity (for example, P(B = 3) is less than or around 0.01). For the analysis, we followed the same procedure aforementioned; that is, computing first the normalized distribution of B given r Na+ ··· Cl−, followed by extracting P(B = 0), P(B = 1), and P(B = 2) for each range of r Na+ ··· Cl− and each solvation region to plot Figure 5b−d. Interestingly, we captured three distinct trends that are consistently displayed by each set of the three curves (blue for E

DOI: 10.1021/acs.jpcb.9b01331 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

as it nears the interface is slightly larger (i.e., ΔP(B) ≤ 0.15). Third, in terms of magnitude of P(B), the curves are in reverse order to those in Figure 5. In other words, P(B) increases before the transition point for B = 0 and after the transition point for B = 1 and B = 2, as the ion pair draws near the membrane surface. On the other hand, P(B) rises after the transition point for B = 0, before the transition point for B = 1, and over the entire range of r∈[2 Å,6 Å] for B = 2, as the ion pair comes close to the graphene surface. Therefore, this finding supports our argument about the pivotal role of phospholipids in the formation of interionic water bridges, and as a consequence, the design of the reaction coordinates for ion pairing at the lipid/water interface should take their headgroups into consideration. Effect of the Inclusion of Polarizable Interactions on Ion-Pair Hydration Structure. For the change in ⟨CN⟩ as r increases, we found the same persistent trend, even in the polarizable model; that is, the curve shape is characterized by the abrupt transition at r = 3.5 Å and the two plateaus before and after the transition point. On average, however, one fewer water molecule is engaged with the first hydration shell of each ion in the polarizable model. This finding is in conformity with the theoretical study of Dang and his colleagues.70 The difference in ⟨CN⟩ between the ion-paired and dissociated states is slightly less; that is, Δ⟨CN⟩ = 1.0−1.5 in the polarizable model as compared to Δ⟨CN⟩ = 1.2−1.8 in the nonpolarizable model. As in the nonpolarizable model, the cation undergoes a larger reduction in the degree of hydration than the anion, as the ion pair comes closer to the membrane surface. However, the cation loses two water molecules ( CN⟩Na+,bulk − CN⟩Na+,near = 1.88 ± 0.21) in the polarizable model, whereas it discards a single or no water molecule ( CN⟩Na+,bulk − CN⟩Na+,near = 0.72 ± 0.23) in the nonpolarizable model. For the anion, the quantity does not alter, to a great extent, with the inclusion of the polarization effect. Therefore, we attribute the higher stability of the CIP state in the polarizable model to the stronger desolvation of the Na+ ion when the ion pair moves to the water/membrane interface. The values of θmax for both ions do not change in the polarizable model; even so, the inclusion of polarizable interactions resulted in two substantial changes, particularly for the Na+ ion, as follows. First, the probability distribution of cos θ becomes bimodal (see the inset in Figure 8a), and we observed strong bimodality arises when the interionic separation is small and the ion pair is close to the interface. Second, the increase in P(cos θmax), as the ion pair moves from the bulk phase to the interface, is approximately 10−20% (see Figure 8a), which is far larger than that found in the nonpolarizable model. These findings are clear indications that the interactions with phospholipids greatly affect the hydration structure of the ion pair, which would in turn exert an impact on ion pairing thermodynamics. These findings are less obvious in the nonpolarizable model. As in the nonpolarizable model, the curve shapes in each Bstate are alike regardless of the solvation location, and ΔP(B) ≤ 0.1 for all the B-states. No difference in the overall shape of the curve is detected between the polarizable and nonpolarizable models, but P(B = 0) becomes much lower and P(B = 2) becomes much higher when r ≤ 3.5 Å. In the polarizable model, slight enhancement in P(B) of the interfacial ion pair is not evidently seen for all the B-states

Figure 6. Two-dimensional probability distribution as a function of the interionic distance (r Na+ ··· Cl−) and the fractional number of water bridges (B) for NaCl ion pairing in the bulk region of solvation.

Figure 7. (a) PMF profiles along the interionic distance (r Na+ ··· Cl−) for ion pairing in the two different regions of solvation (near vs bulk). (b−d) Probability (P) of finding the solvent states of B = 0, B = 1, and B = 2, respectively, as the interionic distance changes. The same color scheme was used as in Figure 1.

shows two PMF curves for a NaCl ion pair at the aqueous interface of graphene (green) and in the bulk environment (blue). Here, we spotted a minute difference in magnitude of the energy barriers (≈0.2 kcal/mol corresponding to ≈7.5% of the peak height). Also, the second minimum occurs at a slightly shorter interionic distance, that is, r = 4.5 Å at the interface (as compared to r = 5.2 Å in the bulk region). Figure 7b−d shows the normalized probability distributions of the ion pair having each B-state when the ions are separated by r in different regions of solvation. Here, we detected the same curve shapes, as in the case of the POPC membrane, but there are several notable differences. First, the location of the transition point is different in B = 1 (occurring at a much shorter interionic distance r = 3.3 Å) and absent in B = 2. Second, the change in P(B) that the ion pair would experience F

DOI: 10.1021/acs.jpcb.9b01331 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

solvent participation coupled to phospholipid headgroups in ion pairing. There still remain many open questions to be answered, including “How are the solvent molecules beyond the FSS involved in ion pairing at an aqueous interface?”; “What is the effect of lipid composition on the thermodynamics and kinetics of an ionic interaction?”; and “Do different types of an ion pair display the same solvation and complexation behavior as Na+Cl− at the interface?”. The generalization of our study to a series of other ion pairs requires the use of different (polarizable) force fields and water models and exploration of multiple reaction coordinates. Since our results are based on one-dimensional PMFs, the generation of multidimensional PMFs (e.g., two-dimensional PMF as a function of ion-pair separation distance and the ion−water CN as in the work of Roy et al.76) will improve our understanding of interfacial ion pairing, and thus, we are currently working toward this direction.

Figure 8. (a) Probability (P) of finding the most probable angle (θmax) for the Na+ ion as the interionic distance (r Na+ ··· Cl−) differs, calculated with the polarizable force field. The inset shows the (normalized) probability distribution of the angle (θ) between a water dipole and the Na+ ion when the interfacial ion pair is separated by 3.0 Å. (b−d) Probability (P) of finding the solvent state of B = 0, B = 1, and B = 2, respectively, as the interionic length (r Na+ ··· Cl−) changes, calculated with the polarizable force field. The same color scheme was used as in Figures 1 and 4.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.9b01331. Figures illustrating the electron density profile of water, rotational relaxation time of water, normalized distribution of ion-pair dipole orientation, and normalized distribution of water coordination numbers (PDF)

(Figure 8b−d), and therefore, the unique transition point r* is absent.





CONCLUSIONS We carried out a detailed analysis of the thermodynamics and solvation of NaCl ion pairing at the surface of a lipid membrane and in a bulk phase to derive an insight into the influence of the presence of the aqueous interface on the simple ion−ion interaction. In the classical, nonpolarizable force field, we observed relatively little consequence of water molecule ordering at the surface of a POPC membrane on the thermodynamics of NaCl ion pairing. To track the origin of such a counterintuitive finding, we examined the effect of electronic polarizability on the interfacial ion pairing and characterized the hydration structure of the ion pair as the interionic separation length and the nearness to the membrane surface vary. A significant degree of stabilization of the CIP state at the interface is found when the polarization effect is included in our simulations. The quantitative analysis of the three coordinates of solvent (namely, the water coordination number of Na+ and Cl− ions, the orientation of water dipoles around the ions, and the number of water bridges between the ions) reveals that, in the nonpolarizable model, there is a minute change in hydration of the ion pair as it approaches to the membrane surface. In the polarizable model, on the other hand, we observed a larger decrease in the degree of water coordination and a larger change in the distribution of water orientation in the first solvation shell of the cation, which contribute appreciably to the stabilization of the paired state. Furthermore, we remark that all the critical changes in solvent structure are not mirrored in the PMF curves, indicating that the critical role of solvent is not well captured by using the conventional reaction coordinate (i.e., interionic distance), as consistent with other theoretical works.74,75 We anticipate that the extension of this study will shed more light on the design of the reaction coordinate. Our study suggests the design of a set of reaction coordinates that account for

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Myong In Oh: 0000-0003-1145-5142 Donald F. Weaver: 0000-0003-2042-6220 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors would like to thank the Krembil Foundation for financial support. D.F.W. acknowledges salary support from a Tier 1 Canada Research Chair. The authors acknowledge SciNet and Compute Canada for providing computing facilities.



ABBREVIATIONS AUC, area under the curve; CGenFF, CHARMM general force field; CIP, contact ion pair; CN, coordination number; COM, center of mass; FSS, first solvation shell; MD, molecular dynamics; MRI, magnetic resonance imaging; PME, particle mesh Ewald; PMF, potential of mean force; POPC, phosphatidylcholine; SSIP, solvent-separated ion pair; WHAM, weighted histogram analysis method



REFERENCES

(1) Cherepanov, D. A.; Feniouk, B. A.; Junge, W.; Mulkidjanian, A. Y. Low Dielectric Permittivity of Water at the Membrane Interface: Effect on the Energy Coupling Mechanism in Biological Membranes. Biophys. J. 2003, 85, 1307−1316. (2) Fumagalli, L.; Esfandiar, A.; Fabregas, R.; Hu, S.; Ares, P.; Janardanan, A.; Yang, Q.; Radha, B.; Taniguchi, T.; Watanabe, K.; et al. Anomalously Low Dielectric Constant of Confined Water. Science 2018, 360, 1339−1342.

G

DOI: 10.1021/acs.jpcb.9b01331 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (3) Debnath, A.; Mukherjee, B.; Ayappa, K. G.; Maiti, P. K.; Lin, S.T. Entropy and Dynamics of Water in Hydration Layers of a Bilayer. J. Chem. Phys. 2010, 133, 174704. (4) Teissie, J.; Prats, M.; LeMassu, A.; Stewart, L. C.; Kates, M. Lateral Proton Conduction in Monolayers of Phospholipids from Extreme Halophiles. Biochemistry 1990, 29, 59−65. (5) Weichselbaum, E.; Ö sterbauer, M.; Knyazev, D. G.; Batishchev, O. V.; Akimov, S. A.; Hai Nguyen, T.; Zhang, C.; Knör, G.; Agmon, N.; Carloni, P.; et al. Origin of Proton Affinity to Membrane/Water Interfaces. Sci. Rep. 2017, 7, 4553. (6) Róg, T.; Murzyn, K.; Pasenkiewicz-Gierula, M. The Dynamics of Water at the Phospholipid Bilayer Surface: A Molecular Dynamics Simulation Study. Chem. Phys. Lett. 2002, 352, 323−327. (7) Murzyn, K.; Zhao, W.; Karttunen, M.; Kurdziel, M.; Róg, T. Dynamics of Water at Membrane Surfaces: Effect of Headgroup Structure. Biointerphases 2006, 1, 98−105. (8) Re, S.; Nishima, W.; Tahara, T.; Sugita, Y. Mosaic of Water Orientation Structures at a Neutral Zwitterionic Lipid/Water Interface Revealed by Molecular Dynamics Simulations. J. Phys. Chem. Lett. 2014, 5, 4343−4348. (9) Chen, X.; Hua, W.; Huang, Z.; Allen, H. C. Interfacial Water Structure Associated with Phospholipid Membranes Studied by Phase-Sensitive Vibrational Sum Frequency Generation Spectroscopy. J. Am. Chem. Soc. 2010, 132, 11336−11342. (10) Björneholm, O.; Hansen, M. H.; Hodgson, A.; Liu, L.-M.; Limmer, D. T.; Michaelides, A.; Pedevilla, P.; Rossmeisl, J.; Shen, H.; Tocci, G.; et al. Water at Interfaces. Chem. Rev. 2016, 116, 7698− 7726. (11) Ingólfsson, H. I.; Carpenter, T. S.; Bhatia, H.; Bremer, P.-T.; Marrink, S. J.; Lightstone, F. C. Computational Lipidomics of the Neuronal Plasma Membrane. Biophys. J. 2017, 113, 2271−2280. (12) van Meer, G. Cellular Lipidomics. EMBO J. 2005, 24, 3159− 3165. (13) van Meer, G.; Voelker, D. R.; Feigenson, G. W. Membrane Lipids: Where They Are and How They Behave. Nat. Rev. Mol. Cell Biol. 2008, 9, 112−124. (14) Holthuis, J. C. M.; Menon, A. K. Lipid Landscapes and Pipelines in Membrane Homeostasis. Nature 2014, 510, 48−57. (15) Le Bihan, D. Diffusion MRI: What Water Tells Us About the Brain. EMBO Mol. Med. 2014, 6, 569−573. (16) Le Bihan, D.; Iima, M. Diffusion Magnetic Resonance Imaging: What Water Tells Us About Biological Tissues. PLoS Biol. 2015, 13, No. e1002203. (17) Bai, R.; Stewart, C. V.; Plenz, D.; Basser, P. J. Assessing the Sensitivity of Diffusion MRI to Detect Neuronal Activity Directly. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, E1728−E1737. (18) Tsurugizawa, T.; Ciobanu, L.; Le Bihan, D. Water Diffusion in Brain Cortex Closely Tracks Underlying Neuronal Activity. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 11636−11641. (19) Qiang, W.; Yau, W.-M.; Schulte, J. Fibrillation of Β Amyloid Peptides in the Presence of Phospholipid Bilayers and the Consequent Membrane Disruption. Biochim. Biophys. Acta, Biomembr. 2015, 1848, 266−276. (20) Zhekova, H. R.; Ngo, V.; da Silva, M. C.; Salahub, D.; Noskov, S. Selective Ion Binding and Transport by Membrane Proteins − a Computational Perspective. Coord. Chem. Rev. 2017, 345, 108−136. (21) Yang, J.; Bonomi, M.; Calero, C.; Martí, J. Free Energy Landscapes of Sodium Ions Bound to DMPC−Cholesterol Membrane Surfaces at Infinite Dilution. Phys. Chem. Chem. Phys. 2016, 18, 9036−9041. (22) Berkowitz, M. L.; Vácha, R. Aqueous Solutions at the Interface with Phospholipid Bilayers. Acc. Chem. Res. 2012, 45, 74−82. (23) Melcrová, A.; Pokorna, S.; Pullanchery, S.; Kohagen, M.; Jurkiewicz, P.; Hof, M.; Jungwirth, P.; Cremer, P. S.; Cwiklik, L. The Complex Nature of Calcium Cation Interactions with Phospholipid Bilayers. Sci. Rep. 2016, 6, 38035. (24) Knecht, V.; Klasczyk, B. Specific Binding of Chloride Ions to Lipid Vesicles and Implications at Molecular Scale. Biophys. J. 2013, 104, 818−824.

(25) Reif, M. M.; Kallies, C.; Knecht, V. Effect of Sodium and Chloride Binding on a Lecithin Bilayer. A Molecular Dynamics Study. Membranes 2017, 7, 5. (26) Neubert, R. Ion Pair Transport across Membranes. Pharm. Res. 1989, 06, 743−747. (27) Miller, J. M.; Dahan, A.; Gupta, D.; Varghese, S.; Amidon, G. L. Enabling the Intestinal Absorption of Highly Polar Antiviral Agents: Ion-Pair Facilitated Membrane Permeation of Zanamivir Heptyl Ester and Guanidino Oseltamivir. Mol. Pharmaceutics 2010, 7, 1223−1234. (28) Karakus, M.; Alpoguz, H. K.; Kaya, A.; Acar, N.; Görgülü, A. O.; Arslan, M. A Kinetic Study of Mercury(II) Transport through a Membrane Assisted by New Transport Reagent. Chem. Cent. J. 2011, 5, 43. (29) Song, I. S.; Choi, M. K.; Shim, W. S.; Shim, C. K. Transport of Organic Cationic Drugs: Effect of Ion-Pair Formation with Bile Salts on the Biliary Excretion and Pharmacokinetics. Pharmacol. Ther. 2013, 138, 142−154. (30) Zhao, H.; Liu, C.; Quan, P.; Wan, X.; Shen, M.; Fang, L. Mechanism Study on Ion-Pair Complexes Controlling Skin Permeability: Effect of Ion-Pair Dissociation in the Viable Epidermis on Transdermal Permeation of Bisoprolol. Int. J. Pharm. 2017, 532, 29−36. (31) Fennell, C. J.; Bizjak, A.; Vlachy, V.; Dill, K. A. Ion Pairing in Molecular Simulations of Aqueous Alkali Halide Solutions. J. Phys. Chem. B 2009, 113, 6782−6791. (32) Zangi, R. Attraction between Like-Charged Monovalent Ions. J. Chem. Phys. 2012, 136, 184501. (33) Pratt, L. R.; Hummer, G.; Garciá, A. E. Ion Pair Potentials-ofMean-Force in Water. Biophys. Chem. 1994, 51, 147−165. (34) Yonetani, Y. Distinct Dissociation Kinetics between Ion Pairs: Solvent-Coordinate Free-Energy Landscape Analysis. J. Chem. Phys. 2015, 143, 044506. (35) Mullen, R. G.; Shea, J.-E.; Peters, B. Transmission Coefficients, Committors, and Solvent Coordinates in Ion-Pair Dissociation. J. Chem. Theory Comput. 2014, 10, 659−667. (36) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (37) Carpenter, T. S.; Kirshner, D. A.; Lau, E. Y.; Wong, S. E.; Nilmeier, J. P.; Lightstone, F. C. A Method to Predict Blood-Brain Barrier Permeability of Drug-Like Compounds Using Molecular Dynamics Simulations. Biophys. J. 2014, 107, 630−641. (38) Simons, K.; Vaz, W. L. C. Model Systems, Lipid Rafts, and Cell Membranes. Annu. Rev. Biophys. Biomol. Struct. 2004, 33 (1), 269− 295. (39) Martinez-Seara, H.; Róg, T. In Biomolecular Simulations: Methods and Protocols; Monticelli, L., Salonen, E., Eds.; Humana Press: Totowa, NJ, 2013; pp 407−429. (40) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2009, 31, 671−690. (41) Vanommeslaeghe, K.; MacKerell, A. D. Automation of the CHARMM General Force Field (CGenFF) I: Bond Perception and Atom Typing. J. Chem. Inf. Model. 2012, 52, 3144−3154. (42) Vanommeslaeghe, K.; Raman, E. P.; MacKerell, A. D. Automation of the CHARMM General Force Field (CGenFF) II: Assignment of Bonded Parameters and Partial Atomic Charges. J. Chem. Inf. Model. 2012, 52, 3155−3168. (43) Yu, W.; He, X.; Vanommeslaeghe, K.; MacKerell, A. D., Jr. Extension of the CHARMM General Force Field to SulfonylContaining Compounds and Its Utility in Biomolecular Simulations. J. Comput. Chem. 2012, 33, 2451−2468. (44) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. H

DOI: 10.1021/acs.jpcb.9b01331 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (45) Miyamoto, S.; Kollman, P. A. Settle: An Analytical Version of the Shake and Rattle Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952−962. (46) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (47) Ewald, P. P. Die Berechnung Optischer Und Elektrostatischer Gitterpotentiale. Ann. Phys. (Berlin, Ger.) 1921, 369, 253−287. (48) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. Constant Pressure Molecular Dynamics Simulation: The Langevin Piston Method. J. Chem. Phys. 1995, 103, 4613−4621. (49) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (50) Kučerka, N.; Nieh, M.-P.; Katsaras, J. Fluid Phase Lipid Areas and Bilayer Thicknesses of Commonly Used Phosphatidylcholines as a Function of Temperature. Biochim. Biophys. Acta, Biomembr. 2011, 1808, 2761−2771. (51) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: New York, U.S.A., 1989. (52) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (53) Jo, S.; Kim, T.; Iyer, V. G.; Im, W. CHARMM-GUI: A WebBased Graphical User Interface for CHARMM. J. Comput. Chem. 2008, 29, 1859−1865. (54) Jiang, W.; Hardy, D. J.; Phillips, J. C.; MacKerell, A. D.; Schulten, K.; Roux, B. High-Performance Scalable Molecular Dynamics Simulations of a Polarizable Force Field Based on Classical Drude Oscillators in NAMD. J. Phys. Chem. Lett. 2011, 2, 87−92. (55) Lamoureux, G.; Harder, E.; Vorobyov, I. V.; Roux, B.; MacKerell, A. D. A Polarizable Model of Water for Molecular Dynamics Simulations of Biomolecules. Chem. Phys. Lett. 2006, 418, 245−249. (56) Giorgino, T. Computing 1-D Atomic Densities in Macromolecular Simulations: The Density Profile Tool for VMD. Comput. Phys. Commun. 2014, 185, 317−322. (57) Torrie, G. M.; Valleau, J. P. Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling. J. Comput. Phys. 1977, 23, 187−199. (58) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011−1021. (59) Timko, J.; Bucher, D.; Kuyucak, S. Dissociation of NaCl in Water from ab Initio Molecular Dynamics Simulations. J. Chem. Phys. 2010, 132, 114510. (60) Gao, J. Simulation of the Na+Cl- Ion Pair in Supercritical Water. J. Phys. Chem. 1994, 98, 6049−6053. (61) Pekker, M.; Shneider, M. N. Interaction between Electrolyte Ions and the Surface of a Cell Lipid Membrane. J. Phys. Chem. Biophys. 2015, 5, 1−7. (62) Cheng, J.-X.; Pautot, S.; Weitz, D. A.; Xie, X. S. Ordering of Water Molecules between Phospholipid Bilayers Visualized by Coherent Anti-Stokes Raman Scattering Microscopy. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 9826−9830. (63) Dreier, L. B.; Nagata, Y.; Lutz, H.; Gonella, G.; Hunger, J.; Backus, E. H. G.; Bonn, M. Saturation of Charge-Induced Water Alignment at Model Membrane Surfaces. Sci. Adv. 2018, 4, No. eaap7415. (64) Teschke, O.; Ceotto, G.; de Souza, E. F. Interfacial Water Dielectric-Permittivity-Profile Measurements Using Atomic Force Microscopy. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2001, 64, 011605. (65) Senapati, S.; Chandra, A. Dielectric Constant of Water Confined in a Nanocavity. J. Phys. Chem. B 2001, 105, 5106−5109. (66) Meneses-Juárez, E.; Rivas-Silva, J. F.; González-Melchor, M. Static Dielectric Constant of Water within a Bilayer Using Recent Water Models: A Molecular Dynamics Study. J. Phys.: Condens. Matter 2018, 30, 195001.

(67) Venkateshwaran, V.; Vembanur, S.; Garde, S. Water-Mediated Ion−Ion Interactions Are Enhanced at the Water Vapor−Liquid Interface. Proc. Natl. Acad. Sci. U. S. A. 2014, 111, 8729−8734. (68) Schweighofer, K.; Benjamin, I. Ion Pairing and Dissociation at Liquid/Liquid Interfaces: Molecular Dynamics and Continuum Models. J. Chem. Phys. 2000, 112, 1474−1482. (69) Hudait, A.; Allen, M. T.; Molinero, V. Sink or Swim: Ions and Organics at the Ice−Air Interface. J. Am. Chem. Soc. 2017, 139, 10095−10103. (70) Dang, L. X.; Schenter, G. K.; Wick, C. D. Rate Theory of Ion Pairing at the Water Liquid−Vapor Interface. J. Phys. Chem. C 2017, 121, 10018−10026. (71) Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. Development of an Improved Four-Site Water Model for Biomolecular Simulations: TIP4P-Ew. J. Chem. Phys. 2004, 120, 9665−9678. (72) Dang, L. X.; Chang, T.-M. Many-Body Interactions in Liquid Methanol and Its Liquid/Vapor Interface: A Molecular Dynamics Study. J. Chem. Phys. 2003, 119, 9851−9857. (73) Lamoureux, G.; Roux, B. T. Modeling Induced Polarization with Classical Drude Oscillators: Theory and Molecular Dynamics Simulation Algorithm. J. Chem. Phys. 2003, 119, 3025−3039. (74) Geissler, P. L.; Dellago, C.; Chandler, D. Kinetic Pathways of Ion Pair Dissociation in Water. J. Phys. Chem. B 1999, 103, 3706− 3710. (75) Ballard, A. J.; Dellago, C. Toward the Mechanism of Ionic Dissociation in Water. J. Phys. Chem. B 2012, 116, 13490−13497. (76) Roy, S.; Baer, M. D.; Mundy, C. J.; Schenter, G. K. Marcus Theory of Ion-Pairing. J. Chem. Theory Comput. 2017, 13, 3470− 3477. (77) Roy, S.; Baer, M. D.; Mundy, C. J.; Schenter, G. K. Reaction Rate Theory in Coordination Number Space: An Application to Ion Solvation. J. Phys. Chem. C 2016, 120, 7597−7605. (78) Roy, S.; Bryantsev, V. S. Finding Order in the Disordered Hydration Shell of Rapidly Exchanging Water Molecules around the Heaviest Alkali Cs+ and Fr+. J. Phys. Chem. B 2018, 122, 12067− 12076. (79) McCaffrey, D. L.; Nguyen, S. C.; Cox, S. J.; Weller, H.; Alivisatos, A. P.; Geissler, P. L.; Saykally, R. J. Mechanism of Ion Adsorption to Aqueous Interfaces: Graphene/Water vs. Air/Water. Proc. Natl. Acad. Sci. U. S. A. 2017, 114, 13369−13373. (80) Vazdar, M.; Pluhařová, E.; Mason, P. E.; Vácha, R.; Jungwirth, P. Ions at Hydrophobic Aqueous Interfaces: Molecular Dynamics with Effective Polarization. J. Phys. Chem. Lett. 2012, 3, 2087−2091. (81) Jungwirth, P.; Winter, B. Ions at Aqueous Interfaces: From Water Surface to Hydrated Proteins. Annu. Rev. Phys. Chem. 2008, 59, 343−366. (82) Bauer, B. A.; Ou, S.; Patel, S. Solvation Structure and Energetics of Single Ions at the Aqueous Liquid−Vapor Interface. Chem. Phys. Lett. 2012, 527, 22−26. (83) Chang, T.-M.; Dang, L. X. Recent Advances in Molecular Simulations of Ion Solvation at Liquid Interfaces. Chem. Rev. 2006, 106, 1305−1322. (84) Ou, S.; Hu, Y.; Patel, S.; Wan, H. Spherical Monovalent Ions at Aqueous Liquid−Vapor Interfaces: Interfacial Stability and Induced Interface Fluctuations. J. Phys. Chem. B 2013, 117, 11732−11742. (85) Caleman, C.; Hub, J. S.; van Maaren, P. J.; van der Spoel, D. Atomistic Simulation of Ion Solvation in Water Explains Surface Preference of Halides. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 6838− 6842. (86) Mancinelli, R.; Botti, A.; Bruni, F.; Ricci, M. A.; Soper, A. K. Hydration of Sodium, Potassium, and Chloride Ions in Solution and the Concept of Structure Maker/Breaker. J. Phys. Chem. B 2007, 111, 13570−13577. (87) Belch, A. C.; Berkowitz, M.; McCammon, J. A. Solvation Structure of a Sodium Chloride Ion Pair in Water. J. Am. Chem. Soc. 1986, 108, 1755−1761. I

DOI: 10.1021/acs.jpcb.9b01331 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (88) Ding, Y.; Hassanali, A. A.; Parrinello, M. Anomalous Water Diffusion in Salt Solutions. Proc. Natl. Acad. Sci. U. S. A. 2014, 111, 3310−3315.

J

DOI: 10.1021/acs.jpcb.9b01331 J. Phys. Chem. B XXXX, XXX, XXX−XXX