Nonlinear response effects in continuum models of the hydration of ions

Hybrid Solvation Model with First Solvation Shell for Calculation of Solvation Free Energy

We present a hybrid solvation model with first solvation shell to calculate solvation free energies. This hybrid model combines the quantum mechanics and molecular mechanics methods with the analytical expression based on the Born solvation model to calculate solvation free energies. Based on calculated free energies of solvation and reaction profiles in gas phase, we set up a unified scheme to predict reaction profiles in solution. The predicted solvation free energies and reaction barriers are compared with experimental results for twenty bimolecular nucleophilic substitution reactions. These comparisons show that our hybrid solvation model can predict reliable solvation free energies and reaction barriers for chemical reactions of small molecules in aqueous solution.

Model for screened, charge-regulated electrostatics of an eye lens protein: Bovine gammaB-crystallin

We model screened, site-specific charge regulation of the eye lens protein bovine gammaB-crystallin (γB) and study the probability distributions of its proton occupancy patterns. Using a simplified dielectric model, we solve the linearized Poisson-Boltzmann equation to calculate a 54×54 work-of-charging matrix, each entry being the modeled voltage at a given titratable site, due to an elementary charge at another site. The matrix quantifies interactions within patches of sites, including γB charge pairs. We model intrinsic pK values that would occur hypothetically in the absence of other charges, with use of experimental data on the dependence of pK values on aqueous solution conditions, the dielectric model, and literature values. We use Monte Carlo simulations to calculate a model grand-canonical partition function that incorporates both the work-of-charging and the intrinsic pK values for isolated γB molecules and we calculate the probabilities of leading proton occupancy configurations, for 4<pH<8 and Debye screening lengths from 6 to 20 Å. We select the interior dielectric value to model γB titration data. At pH 7.1 and Debye length 6.0 Å, on a given γB molecule the predicted top occupancy pattern is present nearly 20% of the time, and 90% of the time one or another of the first 100 patterns will be present. Many of these occupancy patterns differ in net charge sign as well as in surface voltage profile. We illustrate how charge pattern probabilities deviate from the multinomial distribution that would result from use of effective pK values alone and estimate the extents to which γB charge pattern distributions broaden at lower pH and narrow as ionic strength is lowered. These results suggest that for accurate modeling of orientation-dependent γB-γB interactions, consideration of numerous pairs of proton occupancy patterns will be needed.

Origin of Asymmetric Solvation Effects for Ions in Water and Organic Solvents Investigated Using Molecular Dynamics Simulations: The Swain Acity-Basity Scale Revisited

The asymmetric solvation of ions can be defined as the tendency of a solvent to preferentially solvate anions over cations or cations over anions, at identical ionic charge magnitudes and effective sizes. Taking water as a reference, these effects are quantified experimentally for many solvents by the relative acity (A) and basity (B) parameters of the Swain scale. The goal of the present study is to investigate the asymmetric solvation of ions using molecular dynamics simulations, and to connect the results to this empirical scale. To this purpose, the charging free energies of alkali and halide ions, and of their hypothetical oppositely charged counterparts, are calculated in a variety of solvents. In a first set of calculations, artificial solvent models are considered that present either a charge or a shape asymmetry at the molecular level. The solvation asymmetry, probed by the difference in charging free energy between the two oppositely charged ions, is found to encompass a term quadratic in the ion charge, related to the different solvation structures around the anion and cation, and a term linear in the ion charge, related to the solvation structure around the uncharged ion-sized cavity. For these simple solvent models, the two terms are systematically counteracting each other, and it is argued that only the quadratic term should be retained when comparing the results of simulations involving physical solvents to experimental data. In a second set of calculations, 16 physical solvents are considered. The theoretical estimates for the acity A are found to correlate very well with the Swain parameters, whereas the correlation for B is very poor. Based on this observation, the Swain scale is reformulated into a new scale involving an asymmetry parameter Σ, positive for acitic solvents and negative for basitic ones, and a polarity parameter Π. This revised scale has the same predictive power as the original scale, but it characterizes asymmetry in an absolute sense, the atomistic simulations playing the role of an extra-thermodynamic assumption, and is optimally compatible with the simulation results. Considering the 55 solvents in the Swain set, it is observed that a moderate basity (Σ between -0.9 and -0.3, related to electronic polarization) represents the baseline for most solvents, while a highly variable acity (Σ between 0.0 and 3.0, related to hydrogen-bond donor capacity modulated by inductive effects) represents a landmark of protic solvents.

Charge-regulation phase transition on surface lattices of titratable sites adjacent to electrolyte solutions: An analog of the Ising antiferromagnet in a magnetic field

We report a charge-patterning phase transition on two-dimensional square lattices of titratable sites, here regarded as protonation sites, placed in a low-dielectric medium just below the planar interface between this medium and a salt solution. We calculate the work-of-charging matrix of the lattice with use of a linear Debye-Hückel model, as input to a grand-canonical partition function for the distribution of occupancy patterns. For a large range of parameter values, this model exhibits an approximate inverse cubic power-law decrease of the voltage produced by an individual charge, as a function of its in-lattice separation from neighboring titratable sites. Thus, the charge coupling voltage biases the local probabilities of proton binding as a function of the occupancy of sites for many neighbors beyond the nearest ones. We find that even in the presence of these longer-range interactions, the site couplings give rise to a phase transition in which the site occupancies exhibit an alternating, checkerboard pattern that is an analog of antiferromagnetic ordering. The overall strength W of this canonical charge coupling voltage, per unit charge, is a function of the Debye length, the charge depth, the Bjerrum length, and the dielectric coefficients of the medium and the solvent. The alternating occupancy transition occurs above a curve of thermodynamic critical points in the (pH-pK,W) plane, the curve representing a charge-regulation analog of variation of the Néel temperature of an Ising antiferromagnet as a function of an applied, uniform magnetic field. The analog of a uniform magnetic field in the antiferromagnet problem is a combination of pH-pK and W, and 1/W is the analog of the temperature in the antiferromagnet problem. We use Monte Carlo simulations to study the occupancy patterns of the titratable sites, including interactions out to the 37th nearest-neighbor category (a distance of √74 lattice constants), first validating simulations through comparison with exact and approximate results for the nearest-neighbor case. We then use the simulations to map the charge-patterning phase boundary in the (pH-pK,W) plane. The physical parameters that determine W provide a framework for identifying and designing real surfaces that could exhibit charge-patterning phase transitions.

Communication: modeling charge-sign asymmetric solvation free energies with nonlinear boundary conditions

We show that charge-sign-dependent asymmetric hydration can be modeled accurately using linear Poisson theory after replacing the standard electric-displacement boundary condition with a simple nonlinear boundary condition. Using a single multiplicative scaling factor to determine atomic radii from molecular dynamics Lennard-Jones parameters, the new model accurately reproduces MD free-energy calculations of hydration asymmetries for: (i) monatomic ions, (ii) titratable amino acids in both their protonated and unprotonated states, and (iii) the Mobley "bracelet" and "rod" test problems [D. L. Mobley, A. E. Barber II, C. J. Fennell, and K. A. Dill, "Charge asymmetries in hydration of polar solutes," J. Phys. Chem. B 112, 2405-2414 (2008)]. Remarkably, the model also justifies the use of linear response expressions for charging free energies. Our boundary-element method implementation demonstrates the ease with which other continuum-electrostatic solvers can be extended to include asymmetry.

Predicting hydration free energies with chemical accuracy: the SAMPL4 challenge

An implicit solvent model described by a non-simple dielectric medium is used for the prediction of hydration free energies on the dataset of 47 molecules in the SAMPL4 challenge. The solute is represented by a minimal parameter set model based on a new all atom force-field, named the liquid simulation force-field. The importance of a first solvation shell correction to the hydration free energy prediction is discussed and two different approaches are introduced to address it: either with an empirical correction to a few functional groups (alcohol, ether, ester, amines and aromatic nitrogen), or an ab initio correction based on the formation of a solute/explicit water complex. Both approaches give equally good predictions with an average unsigned error <1 kcal/mol. Chemical accuracy is obtained.

A continuum model of solvation energies including electrostatic, dispersion, and cavity contributions

Physically accurate continuum solvent models that can calculate solvation energies are crucial to explain and predict the behavior of solute particles in water. Here, we present such a model applied to small spherical ions and neutral atoms. It improves upon a basic Born electrostatic model by including a standard cavity energy and adding a dispersion component, consistent with the Born electrostatic energy and using the same cavity size parameter. We show that the well-known, puzzling differences between the solvation energies of ions of the same size is attributable to the neglected dispersion contribution. This depends on dynamic polarizability as well as size. Generally, a large cancellation exists between the cavity and dispersion contributions. This explains the surprising success of the Born model. The model accurately reproduces the solvation energies of the alkali halide ions, as well as the silver(I) and copper(I) ions with an error of 12 kJ mol(-1) (±3%). The solvation energy of the noble gases is also reproduced with an error of 2.6 kJ mol(-1) (±30%). No arbitrary fitting parameters are needed to achieve this. This model significantly improves our understanding of ionic solvation and forms a solid basis for the investigation of other ion-specific effects using a continuum solvent model.

Affine-response model of molecular solvation of ions: Accurate predictions of asymmetric charging free energies

Two mechanisms have been proposed to drive asymmetric solvent response to a solute charge: a static potential contribution similar to the liquid-vapor potential, and a steric contribution associated with a water molecule's structure and charge distribution. In this work, we use free-energy perturbation molecular-dynamics calculations in explicit water to show that these mechanisms act in complementary regimes; the large static potential (∼44 kJ/mol/e) dominates asymmetric response for deeply buried charges, and the steric contribution dominates for charges near the solute-solvent interface. Therefore, both mechanisms must be included in order to fully account for asymmetric solvation in general. Our calculations suggest that the steric contribution leads to a remarkable deviation from the popular "linear response" model in which the reaction potential changes linearly as a function of charge. In fact, the potential varies in a piecewise-linear fashion, i.e., with different proportionality constants depending on the sign of the charge. This discrepancy is significant even when the charge is completely buried, and holds for solutes larger than single atoms. Together, these mechanisms suggest that implicit-solvent models can be improved using a combination of affine response (an offset due to the static potential) and piecewise-linear response (due to the steric contribution).

Ion solvation in liquid mixtures: effects of solvent reorganization

Using field-theoretic techniques, we study the solvation of salt ions in liquid mixtures, accounting for the permanent and induced dipole moments, as well as the molecular volume of the species. With no adjustable parameters, we predict solvation energies in both single-component liquids and binary liquid mixtures that are in excellent agreement with experimental data. Our study shows that the solvation energy of an ion is largely determined by the local response of the permanent and induced dipoles, as well as the local solvent composition in the case of mixtures, and does not simply correlate with the bulk dielectric constant. In particular, we show that, in a binary mixture, it is possible for the component with the lower bulk dielectric constant but larger molecular polarizability to be enriched near the ion.

Computation of methodology-independent single-ion solvation properties from molecular simulations. III. Correction terms for the solvation free energies, enthalpies, entropies, heat capacities, volumes, compressibilities, and expansivities of solvated ions

The raw single-ion solvation free energies computed from atomistic (explicit-solvent) simulations are extremely sensitive to the boundary conditions (finite or periodic system, system or box size) and treatment of electrostatic interactions (Coulombic, lattice-sum, or cutoff-based) used during these simulations. However, as shown by Kastenholz and Hünenberger [J. Chem. Phys. 124, 224501 (2006)], correction terms can be derived for the effects of: (A) an incorrect solvent polarization around the ion and an incomplete or/and inexact interaction of the ion with the polarized solvent due to the use of an approximate (not strictly Coulombic) electrostatic scheme; (B) the finite-size or artificial periodicity of the simulated system; (C) an improper summation scheme to evaluate the potential at the ion site, and the possible presence of a polarized air-liquid interface or of a constraint of vanishing average electrostatic potential in the simulated system; and (D) an inaccurate dielectric permittivity of the employed solvent model. Comparison with standard experimental data also requires the inclusion of appropriate cavity-formation and standard-state correction terms. In the present study, this correction scheme is extended by: (i) providing simple approximate analytical expressions (empirically-fitted) for the correction terms that were evaluated numerically in the above scheme (continuum-electrostatics calculations); (ii) providing correction terms for derivative thermodynamic single-ion solvation properties (and corresponding partial molar variables in solution), namely, the enthalpy, entropy, isobaric heat capacity, volume, isothermal compressibility, and isobaric expansivity (including appropriate standard-state correction terms). The ability of the correction scheme to produce methodology-independent single-ion solvation free energies based on atomistic simulations is tested in the case of Na(+) hydration, and the nature and magnitude of the correction terms for derivative thermodynamic properties is assessed numerically.

Free energy simulations of a GTPase: GTP and GDP binding to archaeal initiation factor 2

Archaeal initiation factor 2 (aIF2) is a protein involved in the initiation of protein biosynthesis. In its GTP-bound, “ON” conformation, aIF2 binds an initiator tRNA and carries it to the ribosome. In its GDP-bound, “OFF” conformation, it dissociates from tRNA. To understand the specific binding of GTP and GDP and its dependence on the ON or OFF conformational state of aIF2, molecular dynamics free energy simulations (MDFE) are a tool of choice. However, the validity of the computed free energies depends on the simulation model, including the force field and the boundary conditions, and on the extent of conformational sampling in the simulations. aIF2 and other GTPases present specific difficulties; in particular, the nucleotide ligand coordinates a divalent Mg²⁺ ion, which can polarize the electronic distribution of its environment. Thus, a force field with an explicit treatment of electronic polarizability could be necessary, rather than a simpler, fixed charge force field. Here, we begin by comparing a fixed charge force field to quantum chemical calculations and experiment for Mg²⁺:phosphate binding in solution, with the force field giving large errors. Next, we consider GTP and GDP bound to aIF2 and we compare two fixed charge force fields to the recent, polarizable, AMOEBA force field, extended here in a simple, approximate manner to include GTP. We focus on a quantity that approximates the free energy to change GTP into GDP. Despite the errors seen for Mg²⁺:phosphate binding in solution, we observe a substantial cancellation of errors when we compare the free energy change in the protein to that in solution, or when we compare the protein ON and OFF states. Finally, we have used the fixed charge force field to perform MDFE simulations and alchemically transform GTP into GDP in the protein and in solution. With a total of about 200 ns of molecular dynamics, we obtain good convergence and a reasonable statistical uncertainty, comparable to the force field uncertainty, and somewhat lower than the predicted GTP/GDP binding free energy differences. The sign and magnitudes of the differences can thus be interpreted at a semiquantitative level, and are found to be consistent with the experimental binding preferences of ON- and OFF-aIF2.

Electrostatic solvation energy for two oppositely charged ions in a solvated protein system: salt bridges can stabilize proteins

Born-type electrostatic continuum methods have been an indispensable ingredient in a variety of implicit-solvent methods that reduce computational effort by orders of magnitude compared to explicit-solvent MD simulations and thus enable treatment using larger systems and/or longer times. An analysis of the limitations and failures of the Born approaches serves as a guide for fundamental improvements without diminishing the importance of prior works. One of the major limitations of the Born theory is the lack of a liquidlike description of the response of solvent dipoles to the electrostatic field of the solute and the changes therein, a feature contained in the continuum Langevin-Debye (LD) model applied here to investigate how Coulombic interactions depend on the location of charges relative to the protein/water boundary. This physically more realistic LD model is applied to study the stability of salt bridges. When compared head to head using the same (independently measurable) physical parameters (radii, dielectric constants, etc.), the LD model is in good agreement with observations, whereas the Born model is grossly in error. Our calculations also suggest that a salt bridge on the protein's surface can be stabilizing when the charge separation is < or =4 A.

Fluctuation in electrolyte solutions: the self energy

We address the issue of the self energy of the mobile ions in electrolyte solutions within a general Gaussian renormalized fluctuation theory using a field-theoretic approach. We introduce the Born radii of the ions in the form of a charge distribution allowing for different Born radii between the cations and anions. The model thus automatically yields a theory free of divergences and accounts for the solvation of the ions at the level of continuous dielectric media. In an inhomogeneous dielectric medium, the self energy is in general position dependent and differences in the self energy between cations and anions can give rise to local charge separation in a macroscopically neutral system. Treating the Born radius a as a smallness parameter, we show that the self energy can be split into an O(a(-1)) nonuniversal contribution and an O(a0) universal contribution that depends only on the ion concentration, valency, and the spatially varying dielectric constant. For a weakly inhomogeneous dielectric medium, the nonuniversal part of the self energy is shown to have the form of the Born energy with the local dielectric constant. This self energy is incorporated into the Poisson-Boltzmann equation as a simple means of including this local fluctuation effect in a mean-field theory. We illustrate the phenomenon of charge separation by considering cations and anions of difference sizes and valencies in a periodic dielectric medium.

Influence of nonlinear electrostatics on transfer energies between liquid phases: charge burial is far less expensive than Born model

The widely used Born model describes the electrostatic response of continuous media using static dielectric constants. However, when applied to a liquid environment, a comparison of Born model predictions with experimental values (e.g., transfer free energies and pK(a) shifts) found that agreement is only achieved by using physically unrealistic dielectric constants for proteins, lipids, etc., and/or equally unrealistic atomic radii. This leads to questions concerning the physical origins for this failure of the Born model. We partially resolve this question by applying the Langevin-Debye (LD) model of a continuous distribution of point, polarizable dipoles, a model that contains an added dependence of the electrostatic response on the solvent's optical dielectric constant and both gas- and liquid-phase dipole moments, features absent in the Born model to which the LD model reduces for weak fields. The LD model is applied to simple representations of three biologically relevant systems: (i) globular proteins, (ii) lipid bilayers, and (iii) membrane proteins. The linear Born treatment greatly overestimates both the self-energy and the transfer free energy from water to hydrophobic environments (e.g., a protein interior). By using the experimental dielectric constant, the energy cost of charge burial in either globular or membrane proteins of the Born model is reduced by almost 50% with the nonlinear theory as is the pK(a) shift, and the shifts agree well with experimental trends.

Solvation effect on conformations of 1,2:dimethoxyethane: charge-dependent nonlinear response in implicit solvent models

The physical content of and, in particular, the nonlinear contributions from the Langevin-Debye model are illustrated using two applications. First, we provide an improvement in the Langevin-Debye model currently used in some implicit solvent models for computer simulations of solvation free energies of small organic molecules, as well as of biomolecular folding and binding. The analysis is based on the implementation of a charge-dependent Langevin-Debye (qLD) model that is modified by subsequent corrections due to Onsager and Kirkwood. Second, the physical content of the model is elucidated by discussing the general treatment within the LD model of the self-energy of a charge submerged in a dielectric medium for three different limiting conditions and by considering the nonlinear response of the medium. The modified qLD model is used to refine an implicit solvent model (previously applied to protein dynamics). The predictions of the modified implicit solvent model are compared with those from explicit solvent molecular dynamics simulations for the equilibrium conformational populations of 1,2-dimethoxyethane (DME), which is the shortest ether molecule to reproduce the local conformational properties of polyethylene oxide, a polymer with tremendous technological importance and a wide variety of applications. Because the conformational population preferences of DME change dramatically upon solvation, DME is a good test case to validate our modified qLD model. The present analysis of the modified qLD model provides the motivation and tools for studying a wide variety of other interesting systems with heterogeneous dielectric properties and spatial anisotropy.

Coupling Quantum Monte Carlo to a Nonlinear Polarizable Continuum Model for Spherical Solutes

Starting from the nonlinear dielectric response model of Sandberg and Edholm, we derive an analytical expression of the polarization contribution to the solvation free energy in terms of the electronic density of the solute and the dielectric properties of the solvent. The solvent inhomogeneity is taken into account with the use of a smooth switching function whose spacial variation is established on the basis of how the solvent is arranged around the solute. An explicit form of a local potential representing the solvent effect on the solute is thus obtained by functional analysis. This effective potential can be combined with density functional or quantum chemical methods for the quantum mechanical treatment of the solute. Here, we use quantum Monte Carlo techniques for the solute and apply the method to the hydration of atomic ions finding very good agreement with experimental data.

Interlayer expansion and mechanisms of anion sorption of Na-montmorillonite modified by cetylpyridinium chloride: a Monte Carlo study

To study the change of interlayer structure of a Wyoming-type Na-montmorillonite as a result of the replacement of interlayer Na+ ions by cetylpyridinium (CP+) ions, a series of NPT Monte Carlo simulations of the clay mineral with different contents of CP+, Na+, Cl- ions and water in its interlayer space is carried out. In agreement with conclusions from experimental studies, the simulations show that the CP+ ions form monomolecular, bimolecular, and pseudotrimolecular layers with increasing interlayer contents. Calculated potential energies reveal that clay-organic interactions are stronger than organic-organic interactions in CP+-modified montmorillonite, which is in conformity with observations of earlier thermogravimetric experiments. The simulation results indicate that the pseudotrimolecular arrangement of CP+ ions is a prerequisite for the experimentally observed interlayer sorption of inorganic anions. Furthermore, in the interlayer space with a pseudotrimolecular layer, chloride ions favor the formation of pairs with inorganic rather than organic cations. On the basis of these findings and available experimental data, we propose that the interlayer sorption of inorganic anions from the pore space of an organically modified montmorillonite may occur not only in pairs with organic cations, as suggested earlier, but also in pairs with inorganic cations, which represents a so-far unconsidered and maybe more important mechanism of anion sorption on clay minerals.

Ginzburg-Landau theory of solvation in polar fluids: Ion distribution around an interface

We present a Ginzburg-Landau theory of solvation of ions in polar binary mixtures. The solvation free energy arising from the ion-dipole interaction can strongly depend on the composition and the ion species. Most crucial in phase separation is then the difference in the solvation free energy between the two phases, which is the origin of the Galvani potential difference known in electrochemistry. We also take into account an image potential acting on each ion, which arises from inhomogeneity in the dielectric constant and is important close to an interface at very small ion densities. Including these solvation and image interactions, we calculate the ion distributions and the electric potential around an interface with finite thickness. In particular, on approaching the critical point, the ion density difference between the two phases becomes milder. The critical temperature itself is much shifted even by a small amount of ions. We examine the surface tension in the presence of ions in various cases.

Ion-induced nucleation in polar one-component fluids

We present a Ginzburg-Landau theory of ion-induced nucleation in a gas phase of polar one-component fluids, where a liquid droplet grows with an ion at its center. By calculating the density profile around an ion, we show that the solvation free energy is larger in gas than in liquid at the same temperature on the coexistence curve. This difference much reduces the nucleation barrier in a metastable gas.

Swelling of Na∕Mg-montmorillonites and hydration of interlayer cations: A Monte Carlo study

While the swelling behavior of laboratory-prepared homoionic montmorillonites has been studied extensively in numerous experimental and simulation works, far less attention has been given to much more abundant natural montmorillonites, containing a mix of monovalent and/or bivalent cations in interlayer spaces. We carried out a series of Monte Carlo simulations in order to investigate the reasons for the remarkable difference of experimental swelling patterns of a natural Na-rich/Mg-poor montmorillonite and a homoionic Na-montmorillonite. The simulations reproduced the swelling pattern of a natural montmorillonite, suggesting a mechanism of its hydration different from that of the homoionic montmorillonite. We also found that the differences in size and hydration energy of Mg2+ and Na+ ions have strong implications for the structure and the internal energy of interlayer water. This leads to a difference in the layer spacings of the simulated Mg- and Na-montmorillonites as large as approximately 2.1 A at lower water contents.

Solvation effects in near-critical binary mixtures

A Ginzburg-Landau theory is presented to investigate solvation effects in near-critical polar fluid binary mixtures. Concentration dependence of the dielectric constant gives rise to a shell region around a charged particle within which solvation occurs preferentially. As the critical point is approached, the concentration has a long-range Ornstein-Zernike tail representing strong critical electrostriction. If salt is added, strong coupling arises among the critical fluctuations and the ions. The structure factors of the critical fluctuations and the charge density are calculated and the phase transition behavior is discussed.