Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2467438

Formats

Article sections

- Abstract
- I. Introduction
- II. Electrostatics and Liquid Structure
- III. Solvation of Alkali and Halide Ions
- IV. Molecular Solutes
- V. Discussion and Conclusions
- References and Notes

Authors

Related links

J Phys Chem B. Author manuscript; available in PMC 2008 July 15.

Published in final edited form as:

PMCID: PMC2467438

NIHMSID: NIHMS55341

Sergio A. Hassan, Center for Molecular Modeling, DCB/CIT, National Institutes of Health, U.S. DHHS, Bethesda, Maryland 20892;

The publisher's final edited version of this article is available at J Phys Chem B

See other articles in PMC that cite the published article.

Molecular interactions in solution are controlled by the bulk medium and by the forces originating in the structured region of the solvent close to the solutes. In this paper, a model of electrostatic and liquid-structure forces for dynamics simulations of biomolecules is presented. The model introduces information on the microscopic nature of the liquid in the vicinity of polar and charged groups and the associated non-pairwise character of the forces, thus improving upon conventional continuum representations. The solvent is treated as a polar and polarizable medium, with dielectric properties described by an inhomogeneous version of the Onsager theory. This treatment leads to an effective position-dependent dielectric permittivity that incorporates saturation effects of the electric field and the spatial variation of the liquid density. The non-pairwise additivity of the liquid-structure forces is represented by centers of force located at specific points in the liquid phase. These out-of-the-solute centers are positioned at the peaks of liquid density and exert local, external forces on the atoms of the solute. The density is calculated from a barometric law, using a Lennard-Jones-type solute–liquid effective interaction potential. The conceptual aspects of the model and its exact numerical solutions are discussed for single alkali and halide ions and for ion-pair interactions. The practical aspects of the model and the simplifications introduced for efficient computation of forces in molecular solutes are discussed in the context of polar and charged amino acid dimers. The model reproduces the contact and solvent-separated minima and the desolvation barriers of intermolecular potentials of mean force of amino acid dimers, as observed in atomistic dynamics simulations. Possible refinements based on an improved treatment of molecular correlations are discussed.

Molecular interactions in solution are modulated by the aqueous environment, which may vary broadly in composition. The solvent controls most of the molecular properties in the cell,^{1}^{,}^{2} from chemical reactions of small molecules^{3} to the structure, dynamics, and thermodynamics of macromolecules.^{4}^{–}^{7} In particular, the solvent controls the dynamics of molecular encounters and association/dissociation mechanisms, thus partially dictating the kinetics of biochemical pathways. These processes usually involve the noncovalent binding of a small peptide or organic molecule, either a naturally occurring chemical species (e.g., hormones and neurotransmitters) or a synthetic product (e.g., drugs and molecular probes), to a protein (e.g., enzymes and transmembrane receptors). The statistical distribution of ligand binding modes and the tightness of a molecular complex depend on the forces exerted by the solvent, particularly at the protein–solvent interface. The accurate calculation and prediction of binding modes and free energies may have far-reaching implications in pharmacology and medicine.

Certain molecular processes in solution, such as protein unfolding and molecular dissociation, can be studied experimentally at the level of individual molecules.^{8}^{–}^{10} Atomic force microscopy and laser tweezers are commonly used in single-molecule pulling experiments. The free energy profiles extracted from these studies are shaped by the forces exerted by the aqueous environment. The interpretation of experimental observations and its relation to results from computer simulation are topics of current research.^{11} Direct forces between macromolecules can also be measured using osmotic stress^{12}^{,}^{13} and X-ray scattering experiments in molecular arrays. These techniques provide information on the forces at short intermolecular distances and on their modulation by the aqueous environment.^{14} Studies have suggested that the energy associated with the structuring of water at the molecular interfaces may dominate the interactions at close proximity.^{15}

Reliable calculations of molecular interactions in solution, interpretation of simulations data, and comparison with experimental data depend on the efficiency of the algorithms used for sampling the conformational space and the quality of the model describing the interactions. A substantial part of the work reported in the literature is concerned with the development of sampling techniques and other numerical methods of calculation.^{16}^{–}^{18} Comparatively, less effort has been devoted to improving models of molecular interactions, partially due to the conceptual difficulties in understanding the effects of the solvent at the atomic level.

Molecular interactions in pure water are modulated by bulk-phase electrostatics and by solvent-induced forces.^{19} Bulk electrostatic effects originate in the polarization and reorientation of water molecules in the bulk phase. Solvent-induced forces (SIF) result from the spatial rearrangement of the solvent molecules excluded from the region occupied by the solute.^{20}^{,}^{21} This rearrangement perturbs the hydrogen-bonding (HB) network of water in the solute’s hydration shells, thus generating net forces and torques that affect the solute’s equilibrium structure and dynamics. At the molecular level, hydrophobic and hydrophilic forces are special cases of SIF. The molecular origin of SIF has been studied using computer simulations,^{21}^{–}^{24} and their role in protein–ligand interactions and protein folding has been discussed.^{20}^{,}^{25}^{–}^{27}

A complete description of solvent effects requires both bulk electrostatic and solvent-induced forces to be accounted for. Theoretically robust approaches exist to describe such effects, including earlier and modern theories of dielectric media^{28}^{,}^{29} and theories of liquid structure.^{30}^{–}^{33} These approaches are mathematically and physically rigorous, thus making them ideal as a basis of simplified models for use in computer simulations of macromolecules. Important progress has been made in this direction in recent years,^{34}^{–}^{39} although the computational demands have limited the practical applications to relatively small solutes or model systems. Less rigorous but practical models to represent solvent effects in macromolecules are also available. These models have traditionally focused on solvent electrostatics but have largely neglected SIF, except for simple treatments of hydrophobic interactions. The performance of these empirical models depends mainly on the quality and the number of parameters defined. Most of these models have marginal connection with basic physics theory; so, they cannot be used as the basis for further theoretical developments. For the same reason, the interpretation of results obtained in computer simulations is usually limited, which is the most serious limitation, regardless of their predictive power. The merit of such models is mainly computational speed.

The above discussion highlights the main challenge in developing a model of solvent effects in biomolecules, specifically to find a reasonable balance between physical content and computational efficiency. Attempts at this have been made in the past to describe electrostatic effects of the solvent, for example, by treating the solute as a set of point charges and the solvent as a continuum medium and solving the Poisson–Boltzmann equation numerically^{40}^{,}^{41} or by treating the solute and the liquid as a lattice of point dipoles and solving the resulting equations self-consistently.^{42}^{,}^{43} These models have been used to provide insight into the electrostatic effects of the solvent on molecules of arbitrary shape. Because these approximations have their own practical and conceptual limitations, an alternative to the electrostatic problem has been sought in the dielectric theory of polar/polarizable liquids,^{44}^{,}^{45} which is further discussed in this paper. Physically meaningful, yet practical, solvent models have been used not only in the context of macromolecular electrostatics but also to represent electrostatic effects in quantum chemical calculations of small molecules.^{46}^{,}^{47}

In this paper, a model that combines both electrostatics and solvent-induced forces is presented and discussed in detail. The model is implemented and optimized for the calculation of forces in peptides and proteins. The paper is organized as follows. In section II, the basic theory is described, and the model is introduced. In section III, the model is solved numerically for single ions and ion pairs. These simple solutes are used to assess the suitability of the model for describing electrostatic and liquid-structure forces. In section IV, the application of the model is extended to molecular solutes. In this case, the exact numerical solution carried out in section III is abandoned for the sake of computational efficiency. Instead, a suitable approximation is derived based on the observations described in section III. This simplified model is optimized and applied to amino acid dimers for the calculation of intermolecular potentials of mean force. A summary and a discussion on possible improvements of the model are given in section V.

A solute immersed in a solvent perturbs the structure and dynamics of the liquid with respect to those in the non-perturbed state. These changes in the liquid generate a field that reacts back on the solute, thus modulating its gas-phase properties. A statistical mechanics treatment can formally describe the perturbations in the liquid, from which information on its dielectric response and structure can be obtained. The interest here is on the resulting effects of the perturbed liquid on the perturbing solute. In this paper, electrostatic and liquid-structure forces are treated separately but related to one another through the liquid density distribution.

The dielectric response of a medium can be studied statically or dynamically.^{29} The latter is needed to describe the temporal evolution of the polarization field **P**(**r**, t) upon changes in the solute (e.g., in proton and electron transfer, or chemical reactions in general). The dynamic response of dielectric media has been a topic of extensive experimental and theoretical investigations.^{48}^{–}^{50} Only static fields will be considered here. In a polar and polarizable liquid composed of molecules with isotropic polarizability, α, and permanent dipole moments of magnitude *μ* = |** μ**|, the static polarization field at a position

(1)

where the number density *ρ*(**r**) represents the number of liquid molecules per unit volume at position **r**, and ** μ**is the statistical average of the permanent dipole moments of a liquid molecule at

The polarization field is given, in general, by^{29} **P**(**r**) = *∫ χ*(**r**,**r**′) · **E**(**r**′)d**r**′, where *χ*(**r**, **r**′) is the susceptibility tensor, and **E**(**r**′) is the macroscopic electric field at position **r**′ (i.e., the electric field described by Maxwell’s equations inside continuous media). A phenomenological, yet general molecular theory of dielectric response based on this nonlocal relationship has been derived for spatially and temporally varying electric fields.^{29} If the field varies smoothly with the distance (i.e., if the wavelength is large compared to the characteristic molecular length scale), a local form of **P**(**r**) can be used,

(2)

which defines *ε*(**r**) as the local, static dielectric permittivity of the medium, assumed here to be a scalar quantity. The approximation of eq 2 worsens as the distance to a solute decreases because of the short-wavelength components of the electric field. However, corrections for nonlocality^{53}^{–}^{57} are beyond the scope of this paper, and eq 2 will be used regardless of the proximity to the solute.

Expressions of the internal and directing field as a function of the macroscopic field are derived here in the context of the Onsager theory.^{28}^{,}^{52} The instantaneous reaction field, **R**(**r**), acting on a single water molecule with a fixed orientation of ** μ**, located at the center of a spherical cavity positioned at

(3)

By definition, the difference between the internal and directing fields acting on a single water molecule is the average reaction field calculated over all of the orientations of the dipole, that is, **E**_{i}(**r**) – **E**_{d}(**r**) = **R**(**r**). From this definition and the equation given above for **R**(**r**), the internal field can be written as

(4)

Following refs ^{28} and ^{52}, the directing field satisfies the equation **E**_{d}(**r**) = **E**_{c}(**r**) + *f*(**r**)α**E**_{d}(**r**), where **E**_{c} is the cavity field and *f*α**E**_{d} is the reaction field at the position of the molecule, originating in the dipole α**E**_{d} induced by the directing field itself. **E**_{c} is defined as the total field at the center of a physical cavity surrounding the molecule and is given^{28} by **E**_{c}(**r**) = 3*ε*(**r**)**E**(**r**)/[2*ε*(**r**) + 1]. From these two equations, an expression for the directing field is obtained, which is given by

(5)

Introducing eqs 4 and 5 into eq 1 and using eq 2 yields a closed equation for *ε*(**r**) in the form

(6)

where *E*(**r**) |**E**(**r**)|, *β* 1/*k*_{B}*T* (*T* is the absolute temperature; *k*_{B} is the Boltzmann’s constant), *L*(*x*) = 1/tanh(*x*) − 1/*x* is the Langevin function which arises in the calculation^{28} of the canonical average ** μ** over the dipole orientations in eq 1.

In the Debye theory,^{59} **E**_{i} and **E**_{d} were not distinguished from one another. They were taken as equal to the Lorentz field,^{58} **E**_{L}, which is defined as the electric field inside an ideal cavity surrounding the water molecule, and given by **E**_{L}(**r**) = [*ε*(**r**) + 2]**E**(**r**)/3. The Debye theory was used earlier to derive an expression similar to eq 6 for uniform liquids based on the Lorentz–Debye–Sack (LDS) approximation.^{44}^{,}^{45}^{,}^{51}^{,}^{58} Reaction field corrections^{52} were introduced at a later stage in the derivation. To deal with nonuniform liquids, however, it is convenient to introduce the reaction field and its dependence on *ρ* formally into **E**_{i} and **E**_{d}, as in eqs 4 and 5. The LDS theory has been used in the past and forms the basis of the continuum electrostatics model^{44} reviewed in section IV.1.

Setting *μ* = 0 in eq 6 leads to the dielectric permittivity of a nonpolar, polarizable liquid or to the optical (high-frequency) permittivity *ε*_{∞}of a polar and polarizable liquid. For a uniform liquid, *ρ*(**r**) = *ρ*_{0} = 1/*ν*, where *ν* is the molecular volume (*ν* ~ 29.8 Å^{3} for water, estimated at 25 °C and 1 atm^{60}), an equation relating *ε*_{∞}and α is obtained,

(7)

In the limit *E* →0 (non-perturbed liquid, or *r* →∞), eqs 6 and 7 lead to the Onsager equation,^{51}^{,}^{52}^{,}^{58}

(8)

where *ε*_{0} is the static dielectric permittivity of the liquid in the bulk phase.

Equation 6 can be used to calculate *ε*(**r**) once **E**(**r**) and *ρ*(**r**) are known. These quantities (**E** and *ρ*) are determined by the geometry and the charge distribution of the solute. Theories of liquid structure^{30}^{,}^{31} provide the proper framework to obtain an expression for *ρ*(**r**) in the case of a general solute–liquid interaction potential energy *U*(**r**). A particular approach to this problem is based on a variational principle, where a free energy functional
is defined and minimized with respect to *ρ*′(**r**). The function *ρ*′(**r**) that satisfies
is the equilibrium, nonuniform spatial density distribution *ρ*(**r**) sought. This approach yields a formal expression for the liquid density in the form^{31}

(9)

where *η* depends on the temperature and the mass of the liquid molecules. The functional *c*[*ρ*(**r**);**r**] in eq 9 is the single-particle direct correlation function, which is related to the Orstein–Zernike two-particle direct correlation function *c*^{(2)} through *c*[*ρ*(**r**);**r**] = *∫ρ*(**r**′)*c*^{(2)}(**r**,**r**′)d**r**′. In eq 9, the term -*c*[*ρ*(**r**);**r**]/*β* can be viewed as an effective potential which is determined by the interactions between the liquid molecules. Much of the theoretical work on the theory of liquid structure has been aimed at finding approximate solutions of eq 9 or simplified expressions for *c*^{(2)}(**r**,**r**′) for practical calculations. The simplest approximation is to neglect pair correlations altogether; therefore, the local spatial density is given by a barometric law,^{30}^{,}^{31}

(10)

This functional form will be used here. With this approximation, *V*_{eff}(**r**) is an effective solute–liquid interaction potential energy. A Lennard-Jones (LJ) type function will be used here,

(11)

where *r _{i}* = |

The model introduced in section II is now applied to single ions and ion pairs. For these simple solutes, the free energies and forces can be calculated exactly by numerical integration.

For a monovalent ion, the parameter *σ* in eq 11 is separated into two parts, *σ* = *R*_{I} + *R*_{w}, where *R*_{I} is the radius of the ion in the liquid phase at infinite dilution, and *R*_{w} is an extension identified with the radius of a water molecule (1.38 Å). This definition is consistent with the way the ionic radii are calculated from experimental data.^{63}^{,}^{64} The radius *R*_{I} can be determined experimentally using X-ray or neutron scattering, or through computer simulations. Simulations permit a straightforward representation of the ideal infinite dilution, but their reliability is limited by the quality of the force field. On the other hand, experiments may produce better estimates of ionic radii, but the salt concentrations used may corrupt the interpretation of the results due to the proximity of the ions in the solution. As a compromise, the ionic radii used here were determined by neutron scattering experiments at relatively low salt concentrations. To illustrate, three alkali metal and two halide ions are considered: Li^{+}, Na^{+}, K^{+}, Cl^{−}, and F^{−}, with ionic radii^{65}^{–}^{67} given in Table 1. The reversible work of polarization of the dielectric medium is given by^{68}

(12)

where **E**_{0}(**r**) is the electric field of the solute at position **r** in the absence of the dielectric medium (vacuum field), and the integration is carried out over the entire space; **D**(**r**) = *ε*(**r**)**E**(**r**) is the displacement field. Using eq 2 and **E**(r) = **E**_{0}(*r*)/*ε*(*r*), with |**E**_{0}(*r*)| = 1/*r*^{2} for the ion at *r* = 0, eq 12 yields

(13)

where *ε*(*r*) is obtained from eq 6. The density *ρ*(**r**) in eq 6 is given by eq 10 with the potential of eq 11; the value of *μ* is given by eq 8 with *ε*_{0} = 78.39 and *ε*_{∞}= 1.47 (the values^{60} corresponding to pure water at 1 atm and 25 °C). Here, the polarizability, α, is given by eq 7, where *ρ*(*r*) ≤ *ρ*_{0}, and by the corresponding density-dependent expression (i.e., with *ρ*(*r*) replacing 1/*ν*) where *ρ*(*r*) > *ρ*_{0}. This assumption is physically reasonable since α reaches its finite, gas-phase value as the density decreases, leading to *ε*_{∞}→1 and *ε* →1 as *ρ* →0. On the other hand, as *ρ*(*r*) increases, α decreases in a similar proportion, keeping *ε*_{∞} approximately constant.^{28}

The value of *ξ* in eq 11 for each ion is determined from the standard Gibbs free energy of hydration Δ*G*^{0}. The work of polarization Δ*A*_{e}, calculated from eq 12, corresponds to a Helmholtz free energy of an isochoric thermodynamic process. The relation between both quantities is^{71}^{,}^{72}

(14)

where Δ*G*_{np} is the nonpolar component of the free energy. This term can be decomposed as Δ*G*_{np} = Δ*G*_{cav} + Δ*G*_{sr}, where Δ*G*_{cav} accounts for the reversible work of cavity formation, which can be described by scaled particle theories,^{69}^{,}^{70} and Δ*G*_{sr} is a short-range solute–liquid interaction term. For charged species, Δ*G*_{np} makes a small contribution to Δ*G*^{0}. Since the interest here is in the electrostatic component of Δ*G*^{0}, the nonpolar term is approximated as Δ*G*_{np} ~ 1.325 kcal/mol for all of the ions.^{71}^{,}^{72} Partitioning the experimental Δ*G*^{0} of a neutral salt into components assigned to the individual ions requires knowing the absolute value of the Gibbs free energy of hydration of a proton, Δ*G*_{H}, which is used as a reference.^{72}^{,}^{73} Cluster–ion solvation data have recently been used^{74} to determine Δ*G*_{H}; these provide a more reliable estimate than traditional methods. This value has been used to calculate the hydration energies of alkali and halide ions,^{74}^{,}^{75} and their values are given in Table 1 (a standard conversion term^{76} of ~1.9 kcal/mol has been added to the work of polarization Δ*A*_{e} in eq 14 to relate it to the experimental quantity Δ*G*_{e}0 Δ*G*^{0} − Δ*G*_{np}).

With the above specifications, the integral of eq 13 can be carried out numerically, and *ξ* is chosen for each ion as to satisfy eq 14. The results are presented in Table 1. This procedure is similar to that used in the context of classical force fields to parametrize ion–water van der Waals interactions.^{77} Figure 1 shows the permittivity profiles for Na^{+} and Cl^{−} as a function of the distance to the center of the ion. The figure illustrates the saturation effect of the electric field in the spatially uniform liquid [*ρ*(**r**) = *ρ*_{0}] (dotted line) and the additional modulation of *ε*(*r*) due to a nonuniform liquid density (solid lines). A discussion on the treatment of the ion/liquid interface and the sensitivity of the results to the form of *ε*(*r*) is given in section V. If correlations were introduced, the effective *ε*(*r*) would display an oscillatory behavior with the distance,^{57}^{,}^{78} instead of the monotonic increase displayed in Figure 1. In the absence of correlation, the first peak of the radial distribution function *g*(*r*) is shallow (not shown), and neither a second peak nor a first minimum is observed. Therefore, the coordination number *κ* is calculated as *κ* = *ρ*_{0} *∫r*^{2}*g*(*r*)d*r*, where the upper limit of integration is taken as *R*_{m} = *R*_{I} + 2*R*_{w}, which defines the first hydration shell around the ion. This upper limit would correspond to the position of the first minimum of *g*(*r*), which is usually smaller than *R*_{m}. The coordination numbers obtained with the values of *ξ* calculated above are presented in Table 1, showing reasonable agreement with experimental and theoretical estimates.^{63} The inset of Figure 1 shows Δ*A*_{e} as a function of the integration limit *R*_{u} in eq 13, showing that convergence is achieved only at *R*_{u} > ~30 Å. The contribution of the first hydration shell to Δ*A*_{e} (denoted by Δ*A*_{e}′ in Table 1) is similar in magnitude to the bulk contribution.

For an arbitrary static charge distribution, the relationship between the macroscopic and vacuum fields is given by^{68} **E**(**r**) = **E**_{0}(**r**) − *∫***P**(**r**′)·′|**r** − **r**′|^{−1}d**r**′, for |**r**| > 0. Solving this equation for **E**(**r**) (or the corresponding expression for time-dependent fields^{29}) is the basic goal of a dielectric theory, but it is a difficult task, in general.^{29}^{,}^{79}^{,}^{80} Only for a point charge are **E**(*r*) and **E**_{0}(*r*) linearly related through **E**(**r**) = **E**_{0}(**r**)/*ε*(**r**). In the general case, however, an effective local dielectric *ε*(**r**) can be defined such that the fields satisfy the same linear relation. The approximation introduced here and used earlier^{45}^{,}^{51} is that the effective local permittivity also satisfies eq 6. In this case, eq 6 can be solved for two interacting ions at infinite dilution, and the electrostatic component of the hydration energy, Δ*A*_{e}(*r*), of the pair can be obtained by numerical integration of eq 12, where now *r* is the distance separating the ions. For the purpose of discussion, the pairs Na^{+}–Cl^{−}, Na^{+}–Na^{+}, and Cl^{−}–Cl^{−} will be considered. The values of *ξ* and *σ* found above for single ions are transferred to the corresponding ions in the pair, regardless of the interionic distance, which is a simplification. The hydration energy of each pair can be divided into three components,

(15)

The last two terms are self-energies; Δ*A*_{self,1}(*r*) results from zeroing the charge on ion 2 (i.e., Cl^{−} in the Na^{+}–Cl^{−} pair), while Δ*A*_{self,2}(*r*) results from zeroing the charge on ion 1. The values of *ξ* and *σ* of an ion are retained when its charge is set to zero, which is an assumption. Thus, Δ*A*_{self} measures the change in the hydration energy of an ion a distance *r* from a solvent-excluding cavity. The numerical solution of eq 12, for each value of *r*, is plotted in Figure 2A and 2B. The proximity of a solvent-excluding cavity causes |Δ*A*_{self}(*r*)| of an ion to decrease, as expected, because of the removal of polarizable medium near the charged particle. The term Δ*A*_{int}(*r*) is calculated from eq 13 as the difference between the total Δ*A*_{e}(*r*) and the sum of the self-energies. Figure 2B shows that the self-energy terms decay sharply with the interionic distance and are practically constant for *r* > *σ*_{1} + *σ*_{2}. In contrast, the slow decay of Δ*A*_{e}(*r*) with the distance (Figure 2A) is related to the interaction term Δ*A*_{int}(*r*). The behavior of Δ*A*_{e}(*r*) at small interionic distances, or its limit as *r* →0, is of no physical or practical interest, unless the size of the particle resulting from the limiting process is also specified. This is so because the charge alone does not determine the particle’s hydration energy; information on its size is also needed, which is determined by the electron structure.^{81} The electrostatic component of the solvent force on ion *i* is then *F*_{e,}* _{i}* = −dΔ

The inhomogeneity of the system discussed above arises from the saturation effects of the field and the nonuniform density distribution. Saturation effects are manifested in smaller values of the local dielectric permittivity *ε*(**r**) in the vicinity of charges. The density distribution *ρ*(**r**) of the liquid further modulates *ε*(**r**). If the liquid is assumed to be spatially uniform but saturation effects are still retained, the density modulation of the polarization field is removed, and the integral of eq 12 diverges for solutes containing point charges. To avoid this divergence, a region containing the charges can be excluded from the integration domain. For a single ion, a spherical region of radius *R*_{B} is used; thus *R*_{B} defines the lower limit of integration in eq 13. *R*_{B} is closely related to the Born radius82 commonly used in macroscopic electrostatics but with saturation effects explicitly included.^{44}^{,}^{51}^{,}^{81}^{,}^{83}^{,}^{84} Table 1 shows the Born radii and the values of the permittivity *ε*(*R*_{B}) for the five ions considered here. Following precedence, *R*_{B} can be decomposed into two terms^{44}^{,}^{81}^{,}^{85}^{,}^{86} in the form *R*_{B} = *R*_{I} + *δ*, where *δ* is an extension of the ionic radius. Thus, *R*_{B} is viewed as the size of a cavity formed by the ion in the dielectric medium. Table 1 shows that *δ* is positive and larger for cations than for anions; its value depends on the ion type and not only on its charge. This result has been discussed previously^{44}^{,}^{81}^{,}^{85}^{,}^{86} and interpreted as a manifestation of the different orientations of water molecules in the hydration shells of ions.

An additional approximation to the uniform liquid is the neglect of saturation effects. This leads to the usual Born approximation in macroscopic electrostatics, that is, Δ*A*_{e} ≈ (1/*ε* _{0} − 1)/*R*_{B} for a monovalent ion, with *R*_{B} now being the standard Born radius.^{82} The practical and conceptual problems of introducing a solute/liquid dielectric discontinuity when dealing with microscopic solutes have been recognized and discussed.^{83}^{,}^{84}^{,}^{87}^{–}^{89} The model presented here provides a simple recipe to avoid such discontinuities in a physically meaningful way.

The structure of the liquid is represented by the spatial variations of the liquid density [cf. eq 9]. In addition to bulk electrostatics, the liquid structure generates forces that further modulate interactions. The differential force d**F**(**r**, **r**′) on atom i of the solute at position **r** exerted by the liquid molecules within a volume d*ν*′ at position **r**′ is given by^{20}^{,}^{22}^{,}^{27}

(16)

where *U*(**r**, **r**′) is the interaction potential energy between the atom and a liquid molecule, and d*n*(**r**′) = *ρ*(**r**′)d*ν*′is the number of liquid molecules at **r**′. Here, the potential energy is given by *U*(**r**, **r**′) ≈ *V*_{eff}(|**r** − **r**′|) [cf. eq 11], and the density at **r**′is given by eq 10. The total liquid-structure force **F**_{s,i}(**r**) on atom i is then calculated by integrating d**F**(**r**, **r**′) over **r**′,

(17)

The integration is carried out over the space occupied by the liquid, or it can be extended to the entire space provided that *ρ*(**r**) formally accounts for the region inaccessible to the liquid, for example, another molecule, a lipid bilayer, or an idealized confining boundary, such as a container or solid surface.

The liquid-structure forces on a pair of interacting ions will be analyzed first. In the discussion below, ion 1 is located at the origin of a cylindrical coordinate system (*η*, *θ*, *z*), and ion 2 is at a distance *r* from ion 1 in the +*z* direction. For a given value of *r*, d**F**(*r*, **r**′) is calculated at each point **r**′ and then integrated in space to yield the total forces **F**_{s,}* _{i}*(

(18)

after integrating over the azimuthal angle *θ*. The function *ρ*(*η*, *z; r*) denotes the (*η*, *z*)-dependent density for an interionic separation *r* (a similar equation can be obtained for *F*_{s,1}(*r*); *F*_{s,2}(*r*) = −*F*_{s,1}(*r*) in equilibrium). In eq 18, *F*_{s,2} < 0 implies that the solvent induces an attractive force between the charges, while *F*_{s,2} > 0 implies a repulsive force. The potential associated with *F*_{s,2}(*r*), that is, a potential of mean liquid-structure force type, is given by^{90}

(19)

To solve eqs 18 and 19, the values of *σ*_{1}, *σ*_{2}, *ξ*_{1}, and *ξ*_{2} must be specified. The values reported in Table 1 were obtained from the requirement that the dipole density of the liquid around single ions produce the correct polarization field, hence, the correct values of the hydration energies. Because the same functional form of *V*_{eff} is also used to represent the liquid-structure forces on interacting ions, the values in Table 1 are revisited. To this end, eq 18 was solved numerically for generic values of *σ*’s and *ξ*’s, as to cover all possible cases of relative ion sizes and ion–liquid interaction strengths. Typical results for *F*_{s,2}(*r*) and _{s,2}(*r*) are shown in Figure 3 for Na^{+}–Cl^{−}; the values *ξ*_{1}, *ξ*_{2}, *σ*_{1}, and *σ*_{2} are taken directly from Table 1, but the energies are rescaled as to preserve the relative liquid–ion interaction strength, *ξ*_{1}/*ξ*_{2}. The values used in Figure 3 are *ξ*_{1} = 0.3 and *ξ*_{2} = 1.5 kcal/mol (ion 2 is Cl^{−}). At a large interionic separation, the force is attractive. As the distance decreases, the force reaches a minimum at *r* = *r*_{m} > *σ*_{1} + *σ*_{2}. It vanishes at *r*≈ *σ*_{1} + *σ*_{2} (point of contact of the solvation spheres) and becomes repulsive. It then reaches a maximum at *r* = *r*_{M} < *σ*_{1} + *σ*_{2} and starts decreasing again as the ions move closer to each other. In general, the force displays a second (positive or negative) minimum at shorter distances. Depending on the parameters, the force might show a second maximum as the distance decreases even further. The force vanishes at *r* = 0, as expected. For small values of *ξ*_{1} and *ξ*_{2}, the liquid becomes more uniform, although the minimum at *r*_{m} and the maximum at *r*_{M} can still be identified. These observations suggest that the forces associated with the desolvation barriers and the solvent-separated minima of interionic potentials of mean force (PMF)^{91}^{–}^{93} can be reproduced with a suitable choice of *σ*’s and *ξ*’s. However, a complete representation of the effects of the liquid structure on the PMF would require one to account explicitly for the interactions between liquid molecules. The approach developed here captures the effects that can be embedded in an effective liquid–solute potential [cf. eq 11].

In section IV.2, an analysis of the liquid density as a function of the interionic distance and the induced effects leading to the modulation of the liquid-structure force shown in Figure 3 will be presented.

The model presented in section II and solved numerically in section III for the specific case of ions and ion pairs can be applied to molecular solutes as well. Computationally efficient grid-based methods can be used to integrate eqs 6, 12, and 17 in arbitrary geometries. However, simplifications are still needed because the computational task may be too demanding in simulations of macromolecules or other biological applications. Such simplifications are presented in this section.

The total noncovalent force, **F**_{i}, on an atom i of a molecule is divided into four terms

(20)

where **F**_{c,}* _{i}* and

A continuum model of solvent electrostatic effects based on a uniform version of eq 6 has been reported.^{44}^{,}^{45} The model is based on a superposition of screened Coulomb potentials (SCP) and describes the interaction (Δ*A*_{int}) and the self-energy (Δ*A*_{self}) terms of the hydration energy of a molecule. To obtain expressions for these terms in molecular solutes, the hydration energies of one charge and of two interacting charges were first analyzed. The insight from these studies provided the basis to derive an expression for the hydration energy of *N* interacting charges. These steps have been discussed in previous publications,^{44}^{,}^{45}^{,}^{83} and a summary is given below.

In this case, the electric potential in a polar/polarizable liquid is given by (*r*) = _{0}(*r*)/*D*(*r*) where _{0}(*r*) is the potential in the vacuum at a distance *r* from the charge. The scalar function *D*(*r*) is the screening function, which is related to *ε*(*r*) through the definition of electric potential, **E**(**r**) = −(**r**). For a point charge, the relationship between both quantities is^{95} *ε* = *D*/[1 + (*r*/*D*)d*D*/d*r*]. This equation can be solved for *D*(*r*) once the dielectric permittivity *ε*(*r*) is known from experimental data or from theoretical considerations.^{44}^{,}^{45} In the existing SCP model, *ε*(*r*) is given by eq 6 in ref ^{44}, which is based on the LDS theory. If the liquid density is not uniform, eq 6 should be used instead (this case will not be addressed here). The hydration energy of the charge is given by eq 13,

(21)

where *R*_{B} is the Born radius discussed in section III.2.1, which includes the saturation effect of the field. In this form, Δ*A*_{self} can be evaluated with little computational cost provided that the functional form of *D*(*r*) is known and a prescription to determine *R*_{B} is given. A sigmoidal function of *D*(*r*) that reproduces *ε*(*r*) in polar/polarizable liquids has been reported^{58}^{,}^{96} and discussed.^{44}^{,}^{45} This function depends on a single parameter *λ* that controls the rate of increase of *D*(*r*) with the distance to the charge.^{45} In the SCP model, the single parameter *λ* embeds the physical quantities of the system contained in the uniform version of eq 6: *ε*_{0}, *ε*_{∞}, *T*, α, *μ*, *ν*, and *Q*. The assignment of *R*_{B} follows the same idea discussed in section II for ionic sources, that is, a charge radius *R*_{Q} plus an extension *δ* (*R*_{B} = *R*_{Q} + *δ*).

For two charges in a polar/polarizable liquid, the expressions for the interaction energy and the self-energy terms of eq 15 are

(22a)

(22b)

where *i* = 1 or 2, and *D*(*r*) in eq 22a differs quantitatively from those in eq 22b but has the same sigmoidal behavior discussed above for a single charge.^{45} The screening functions in the self-energies are the same as those in eq 21, so the hydration energy of each charge is recovered as *r* →∞. The dependence of Δ*A*_{self,}* _{i}* with

(23a)

(23b)

where *R*_{1B} and *R*_{2B} are the Born radii *R*_{B} of eq 21 for each charge; *ξ*_{12} and *ξ*_{21} are the fractions of the solvent-accessible surface area of Q_{1} and Q_{2}, respectively. The factors 1 − *ξ*_{12} and 1 − *ξ*_{21} are the fractions of the solvent-excluded surface area of each charge. The quantity *R*_{12} accounts for the effects on Q_{1} of the absence of polarizable medium in the region occupied by Q_{2}; *R*_{21} accounts for the same effect on Q_{2} due to the presence of Q_{1}. The linear combination of eq 23a,b is the simplest way to account for the dependence of Δ*A*_{self,}* _{i}*(

For a molecular solute composed of *N* atoms, eq 22a,b can be generalized as^{44}

(24a)

(24b)

where *D* is the screening functions of the fully hydrated solute [the same as in eq 22], and *D*′ is the screening functions of the gas-phase solute. In general, *D*′> 1 since a molecule in vacuum is itself a dielectric medium (if *D*′ = 1, eq 24 is a trivial generalization of eq 22). Models have been proposed that represent a molecule in vacuum as a dielectric composed of charges and polar/polarizable dipoles.^{97}^{,}^{98} These models are based on the theory of dielectric media discussed in section II and, thus, have a direct relation with the model introduced in this paper. A connection between the electron structure of a molecule and the screening functions in solution has also been discussed.^{45} In the SCP continuum model, *R*_{B,}* _{i}* is given by a generalization

From eq 24, the total electrostatic energy *E*_{T} of the hydrated molecule is given by^{44}

(25)

As defined in eq 20, the electrostatic component of the solvent force on an atom *i* of the solute is then calculated as **F**_{e,}* _{i}* = −

(26)

with the interaction and self-energy components given by

(27a)

(27b)

where **r*** _{ij}* =

To derive a computationally efficient model for **F**_{s,}* _{i}*, a strategy similar to the one summarized above for

To gain insight into the dependence of *F*_{s,2}(*r*) on *r* shown in Figure 3, the integrand of eq 17 is analyzed. A force density function
for a fixed *r* can be defined from eq 17 as
. Figure 4 shows a gray-scale plot of
on a plane containing the pair of ions Na^{+}–Cl^{−} at three representative interionic distances, 2, 5, and 7 Å; the black represents regions of the liquid eliciting repulsive interionic forces (*F*_{s,2} > 0), while the white shows the regions of attractive forces. An isolated ion (*r* →∞; left panel in Figure 5) is surrounded by four distinct regions, generating forces that exactly cancel each other in equilibrium (*F*_{s,2} = 0). This follows from the functional form of *V*_{eff} after recognizing that a liquid molecule located at a distance *r*′< *σ* from the center of the ion repels the ion, while a molecule at *r*′> *σ* attracts the ion. As the distance between the two ions decreases, the balance of forces on the ions breaks down due to the spatial distortion of *ρ*(**r**) (along the *z*-axis): At *r* = 7 Å (Figure 4A), the distortion of the density in the attractive and repulsive regions is such that overall attraction is favored (a in Figure 3). This net attractive force originates mainly in a higher liquid density developed in the attractive region between the ions, relative to that in the repulsive region. It can be shown that for *r* > *σ*_{1} + *σ*_{2}, there are two peaks of *ρ*(**r**) on the *z*-axis and in the attractive region between the ions. The positions of the peaks and their densities depend on the value of the *σ*’s and *ξ*’s assigned to both ions. As the ions approach each other so do the density peaks, and they finally coalesce when the solvation spheres start overlapping at *r* ~ *σ*_{1} + *σ*_{2} (5.6 Å in this case; b in Figure 3). These density peaks will play a central role in the model of liquid-structure forces in molecular solutes developed in section IV.2.2.2. For *r* < *σ*_{1} + *σ*_{2}, the local maxima of the density between the ions disappear, becoming a continuum of saddle points corresponding to the intersection line of the solvation spheres. At *r* = 5 Å (Figure 4B), the density in the attractive region between the ions has been removed away from the *z*-axis, leaving an unbalanced repulsive force (c in Figure 3). At *r* = 2 Å (Figure 4C), high liquid density begins to develop in the attractive region at the left of ion 1, while the density in the repulsive region between the ions has been completely removed; thus, an attractive net force results (d in Figure 3). Spherical symmetry is regained at *r* = 0 (*F*_{s,2} = 0), but with the density distribution controlled by both ions.

(A) Schematic representation of the regions defining the double-shell model of liquid-structure forces for ion pairs at *r* < *σ*_{1} + *σ*_{2}. (B) Gray-scale representation of the liquid-structure force density function,
, for an isolated **...**

The main characteristics of *F*_{s,2}(*r*) versus *r*, plotted in Figure 3, are preserved if the domain of the integral in eq 17 is restricted to regions where *ρ* > *ρ*_{th}, with *ρ*_{th} > 1 being an arbitrary threshold. As *ρ*_{th} increases, |*F*_{s,2}(*r*)| becomes smaller for all *r*, as expected, but both the minimum at *r*_{m} and the maximum at *r*_{M} remain well defined. This suggests that the modulation of *F*_{s,2}(*r*) versus *r* originates mainly in regions of high liquid density. Because of the functional form of *V*_{eff} [cf. eq 11], these regions are located near the solute and provide the short-range character of the liquid-structure forces. If the regions of high density are identified, the integration domain in eq 17 shrinks, and the calculation of the force becomes computationally tractable, even in complex solutes. These observations and the analysis of the force density,
, discussed above suggest a simple model for liquid-structure forces in a system composed of two charges. In this model, the total force *F*_{s,i} on an ion will be associated only with the peaks (for *r* > *σ*_{1} + *σ*_{2}) and to the continuous region of stationary points around the ions (for *r* <*σ*_{1} + *σ*_{2}). The forces exerted by the peaks will be discussed in the more general case of molecular solutes (section IV.2.2 below); so, the focus here is only on the forces exerted by the high-density regions at *r* < *σ*_{1} + *σ*_{2}. These regions are unique to the cylindrical symmetry of an ion pair and are generally not present in molecules. For completeness, however, a model is derived here to account for such effects. This model will play a central role when explicit ions are included as part of the molecular solute (e.g., as counterions) in dynamic simulations using a continuum solvent model (to be reported).

To derive a model of *F*_{s,i} for *r* < *σ*_{1} + *σ*_{2}, the solvent surrounding ion 2 is divided into the regions shown in Figure 5A. This partition is suggested by the force density distribution,
, on an isolated ion shown in Figure 5B and by the behavior of
discussed in section III.2.4 and illustrated in Figure 4, on ion pairs. An external shell of thickness *d*_{e} is formed by two hemispherical shells; *V*_{e}^{+} exerts a repulsive force, and *V*_{e}^{−} exerts an attractive force. Similarly, the internal shell of thickness *d*_{i} is formed by the hemispherical shells *V*_{i}^{+} and *V*_{i}^{−} (the index e/i indicates external/internal shell, while the signs +/− indicate the direction of the induced force on the *z*-axis, that is, repulsive/attractive). In the discussion below *V*_{e}^{+}, *V*_{e}^{−}, *V*_{i}^{+}, and *V*_{i}^{−} will denote not only the corresponding regions but also their volumes. If ion 1 is centered at *r* = 0 and the center of ion 2 is moved from *r* = 0 in the +*z* direction, the volumes *V*_{e}^{−} and *V*_{i}^{+} change according to the extent of overlap of the solvation spheres (of radii *σ*_{1} and *σ*_{2}). These overlapping regions are labeled V_{e} and V_{i}. If *λ*_{e} > 0 and *λ*_{i} > 0 are parameters representing the force per unit volume exerted by each shell, then

(28)

where the explicit dependence on *r* was indicated in *V*_{e}^{−}(*r*) and *V*_{i}+(*r*). Since the volumes of the hemispherical shells are *V*_{e}^{+} = *V*_{e} + *V*_{e}^{−} and *V*_{i}^{−} = *V*_{i} + *V*_{i}^{+}, the force is given by

(29)

where the volumes *V*_{e}(*r*) and *V*_{i}(*r*) can be calculated as the difference of two volumes, *V*_{e}(*r*) = *V*_{1e}(*r*) − *V*_{12}(*r*) and *V*_{i}(*r*) = *V*_{12}(*r*) − *V*_{1i}(*r*), where *V*_{1e}(*r*) is the volume of the region defined by the intersection of the solvation sphere of ion 1 and the sphere defined by the outermost surface of the external shell around ion 2 (of radius *σ*_{2} + *d*_{e}). *V*_{12}(*r*) is the overlapping volume of the solvation spheres of ion 1 and ion 2, and *V*_{1i}(*r*) is the overlapping volume of the solvation sphere of ion 1 and the sphere defined by the innermost surface of the internal shell around ion 2 (of radius *σ*_{2} − *d*_{i}). These volumes can then be written as

(30)

where the index *l* is e, 2, or i,
, and *θ*(*x*) is the Heaviside function defining the values of *r* for which overlap between the spheres exist. The functions *x*_{1}* _{l}* and

In analogy with the case of ions, two kinds of systems will be considered, (1) single amino acids and (2) amino acid dimers. Single amino acids are used to gain insight into the structure of the surrounding liquid; amino acid dimers are used to investigate the restructuring of the liquid upon dimer dissociation and its effects on the intermolecular forces.

The focus of this section is on the forces induced by the structure of the liquid surrounding polar and charged groups. For these groups, the barometric law of eq 10 is expected to be a better approximation than for groups interacting more weakly with the solvent. In the case of nonpolar solutes, a more explicit treatment of *c*^{(2)} in eq 9 is needed because the effects related to liquid–liquid interactions may not be properly captured in an effective liquid–solute potential. Formally, however, the treatment of liquid-structure forces described in this paper is general and can be applied to polar and nonpolar solutes, provided that a suitable density function *ρ*(**r**) is available.

Experimental information on the structure of liquids can be obtained from radial distribution functions, for which X-ray or neutron scattering methods are commonly used. Because of the practical difficulties, however, experimental studies have been limited mainly to pure liquids or relatively small solutes (e.g., ions). Computer simulations can overcome these difficulties and provide a great deal of information on the physics of hydration. Despite their obvious limitations, computational methods have a unique advantage over experimental techniques in that the microscopic structure of the liquid around the solutes, its spatial restructuring upon changes of a solute conformation, and the forces elicited by the liquid can be studied in atomistic detail. Molecular dynamics (MD) simulations are used here to elucidate the structure of water surrounding single polar/charged amino acid side chains. The computational setup of the simulations is the same as described in ref ^{90}. A classical, nonpolarizable, rigid three-point water model has been used to represent the liquid. The MD simulations were carried out at a temperature of 25 °C; the volume and number of water molecules in the simulation box were kept constant, with a density *ρ*_{0} ≈ 0.03325 Å^{−3} (corresponding to ~0.993 g/cm^{3}). Four charged (Arg^{+}, Lys^{+}, His^{+}, and Asp^{−}) and four polar, net-neutral (His, Asn, Tyr, Ser) amino acids were modeled as described.^{90} Spatial distribution functions *g*(**r**) = *ρ*^{−1}(**r**) = *ρ*_{0} ^{−1} *δN*(**r**)*δV*^{−1} were calculated for each solute, where *δN*(**r**) is *ρ*_{0} the average number of water molecules within an element of volume *δV* located at position **r**, given by *δN*(**r**) = *τ*^{−1}*∫N*(**r**, *t*)d*t*, with N(**r**, *t*)(**r** = Σ δ−**r**_{i}(*t*)), and the sum runs over all of the water molecules; *δ*(*x*) is 1 for *x* = 0 (i.e., if the oxygen atom of water i is within *δV* at time *t*) and 0 otherwise; *τ* is the total simulation time, which is long enough to ensure convergence of the distribution functions. The positions of the peaks of *ρ*(**r**) around each molecule were calculated as described in ref ^{19}. The focus here is only on the first peaks, defined as the set of points {**r**_{f}}_{f=1,}* _{n}*, consistent with direct H-bond interactions with polar protons or acceptor atoms of the solute’s functional group.

The positions of the peaks around each molecule are shown in Figure 7 (red dots). The peaks distributions can be understood in terms of the spatial symmetry of the functional groups. For Arg^{+} (A), they appear in the plane of the ring, near the hydrogen atoms, in the direction of the N–H bonds (2.8 < *g*_{f} < 3.4). In Lys^{+} (B), the peaks are in a tetrahedral arrangement with respect to the central nitrogen (2.8 < *g*_{f} < 3.3). In His^{+} (C), the peaks are in the plane of the ring (*g*_{f} of 2.8 and 3.7); in both cases, the peaks appear in the direction of the N–H bonds. In Asp^{−} (D), three peaks are located in a tetrahedral arrangement around each of the oxygen atom in the –COO^{−} group; an additional peak is shared by both oxygen atoms (1.9 < *g*_{f} < 3.4). For the polar amino acids, the values of *g*_{f} are, in general, smaller than those for charged molecules. In Asn (H) three peaks are located around the oxygen atom, with a tetrahedral arrangement similar to those around the oxygen atoms in Asp^{−}. Two other peaks are observed near the protons in the direction of each N–H bonds (1.7 < *g*_{f} < 2.6). The removal of a proton in His^{+} causes one of the peaks to split into two of lower densities, one at each side of the His^{0} (G) ring (1.8 < *g*_{f} < 2.5). The positions of the peaks in Tyr (F) and Ser (E) are similar, as expected from the geometry of their functional groups (the double bond character of the –CO– bond in Tyr has no direct effect on the results of the classical simulations). In both cases, two peaks are observed close to the oxygen atom, and a third peak is bonded to the polar proton in the direction of the O–H bond. For both molecules, 1.8 < *g*_{f} < 3.6, where the larger value corresponds to the peak close to the proton. The peak distributions are elicited by the geometry of the side chains only and contain no information on the quantum mechanical nature of the solute or water molecules.

Peaks of liquid density surrounding charged and net-neutral polar amino acids. Red dots: first peaks calculated from a time average of water dynamics in an atomistic simulation of the liquid. White dots: peaks calculated numerically from the liquid density **...**

Parametrization. The densities and positions of the peaks identified above are used to optimize the solute–liquid potential *V*_{eff}(**r**) of eq 11. In principle, each atom *i* of the solute can be assigned a parameter, *ξ _{i}* and

(31)

where **r**_{p} and *g*_{p} are the position and density of a peak of *ρ*(**r**) calculated from eq 13, using the potential *V*_{eff}(**r**) of eq 16, which contains the set {*ξ*, *σ*} to be optimized. The peaks are formally defined as the set of *m* points {**r**_{p}}_{p=1,}* _{m}*, where

Figure 7 shows the location of the peaks after optimization of the potential *V*_{eff} (white dots). In all cases, the distributions of {**r**_{p}} closely follow those of {**r**_{f}}, although some deviations are observed in Asp^{−} (in particular, the distances of {**r**_{p}} to the oxygen atoms in the –COO^{−} group tend to be shorter than those for {**r**_{f}}). These deviations are not surprising since, for Asp^{−}, 28 scalar quantities are optimized with only 4 parameters. Additional peaks may also appear in some cases. In His^{+}, two more peaks are observed at each side of the ring’s plane on a line passing through the ring’s geometric center; a similar peak is observed at one side of the ring in Arg^{+}. The root-mean-square deviation (rmsd) of {**r**_{p}} with respect to {**r**_{f}}, averaged over the eight molecules, is 0.8 Å. The values of {*g*_{p}} are also well reproduced, yielding an average rmsd with respect to {*g*_{f}} of 0.4. These results show that a barometric law with an optimized Lennard-Jones-type potential, *V*_{eff}(**r**), yields satisfactory results despite the relatively small set of parameters used.

A systematic study of the intermolecular H-bond interactions between polar/charged amino acid side chain dimers in pure water has been reported.^{90} Those simulations were aimed at (i) quantifying the strength of the interactions in solution and obtaining information on the minima and barriers of the intermolecular potentials of mean force (PMF), (*r*), and (ii) understanding the structural and dynamical properties of the liquid surrounding the solutes and the molecular origin of the forces exerted by the liquid. Point (i) was reported in ref ^{90} where detailed information on the PMF, a quantity which is not accessible experimentally, was obtained. Techniques are currently being developed to experimentally explore the energy landscapes of molecules in solution and provide more direct information on the PMF.^{105}^{–}^{107} Point (ii) was partially addressed in ref ^{19}, a study that provided insight into the liquid structure and its effects on the intermolecular interactions. The theoretical model presented in this paper is based on the observations reported in refs ^{19} and ^{90}, which can be viewed as the experimental component of this work.

The PMF of a dimer is characterized by the positions of the contact minimum, *r*_{cm}, the desolvation barrier or transition state, *r*_{ts}, and the solvent-separated minimum, *r*_{ss}, along with the corresponding values of the potentials, that is, _{cm} (*r*_{cm}), _{ts} (*r*_{ts}), and _{ss} (*r*_{ss}). A simplified yet complete solvent model requires these minima and barriers to be accounted for. The minima and the desolvation barriers are determined not only by the electrostatic forces exerted by the bulk liquid but also by the SIF originating in the microscopic nature of the solvent. The SCP–ISM, which is based on the SCP continuum electrostatics model reviewed in section IV.1, empirically^{100} accounts for the values of _{cm} obtained in atomistic dynamic simulations of the liquid.^{90} This empiricism is needed because, if the correct value of _{cm} is sought and only an electrostatic model is used, the correction for the missing SIF effects on _{cm} must be incorporated in some way. Some of the problems observed in macromolecular simulations seem to stem from a failure to appreciate this. Barriers and minima of intermolecular potentials have been observed with a macroscopic electrostatic model and a suitable choice of the solute/solvent dielectric interface.^{108} In the SCP–ISM, the correction is made by adjusting the Born radii [cf. eq 23] of protons shared by the donor and the acceptor atoms in hydrogen bonds. This empirical correction also reproduces barriers and additional minima that are qualitatively similar to the desolvation barriers and solvent-separated minima of the PMF (not shown). However, such additional modulation is a byproduct of the approach employed to control the values of _{cm} at *r*_{cm}. A physically complete solvent model should account explicitly for the SIF responsible for the form of the PMF; in such a case, the empiricism discussed above can be removed.

In analogy with the case of ion pairs, the liquid-structure forces are here assumed to originate mainly in regions of high liquid density. In molecular solutes of arbitrary shape, such regions correspond to peaks located in the vicinity of the molecule. In this case, the liquid-structure force, **F**_{s,}* _{i}*, associated with the peaks is responsible for the non-pairwise additivity of the liquid forces and is discussed in detail in this section. From eq 17,

(32)

where **r** is the position of atom i in the solute, *n* is the total number of peaks in the liquid, and *N*(**r*** _{j}*) is the number of water molecules at the position

Introducing eq 11 into eq 32, the total force exerted by the peaks on atom i of the solute is given by

(33)

There is no conceptual need to restrict the summation in eq 33 to peaks located around polar and charged groups because eq 17 is general. However, as discussed in section IV.2.3, possible limitations arise from the barometric law and the effective potential *V*_{eff}. For nonpolar solutes, the force, **F**_{s,}* _{i}*, calculated from eq 32 is strictly zero, but it is different than zero if eq 17 is integrated numerically. Therefore, a nonpolar force,

The intermolecular force and associated potential can now be calculated for the amino acid dimers studied in ref ^{90}. The values of the *ξ*’s and *σ*’s obtained for single amino acids are transferred to the case of the dimers. For illustration purpose,
, but, in general, it should be considered a parameter of the model, as discussed below. For the purpose of discussion, only one representative dimer (Asp^{−}–Arg^{+}) is considered here; a detailed analysis of all of the dimers studied in ref ^{90} will be reported elsewhere. The initial setup of the dimer conformation and the protocol for the calculation of forces and potentials are as described^{90} (except for the distance step, which is set here to Δ*r* = 0.02 Å).

The restructuring of the liquid as the dimer dissociates from close contact is illustrated in Figure 8. The peaks of *ρ*(**r**) are shown for representative intermolecular distances. The redistribution of peaks closely follows the changes observed in the dynamics simulations.^{19} Figure 9A shows the component *F*_{s}(*r*) **F**_{s}(*r*)· (thin line) and the associated potential _{s}(*r*) (thick line) for the dimer. is a unit vector in the direction of the line connecting the centers-of-mass of the two monomers. The value *r* is the proton-acceptor distance, and the intermolecular force is calculated as in ref ^{90}: **F**_{s}(*r*) = [**F**_{s,A}(*r*) − **F**_{s,D}(*r*)]/2, where **F**_{s,A} and **F**_{s,D} are the total force on the acceptor and donor monomer, respectively, evaluated at their centers-of-mass. The total force on a monomer is obtained from the summation of **F**_{s,}* _{i}* of eq 33 over all atoms i of the monomer.

Restructuring of the liquid as the Asp^{−}–Arg^{+} dimer dissociates from close contact. The density peaks’ redistribution is shown for representative proton-acceptor distances (in Å): (A) *r* = 1.8; (B) *r* = 3.2; (C) *r* = 4.6; the **...**

Figure 9A also shows a maximum of *F*_{s}(*r*) at *r* = *r*_{M}. This maximum is similar to that in Figure 3 for the case of ion pairs. As discussed above, *F*_{s}(*r*) at *r*_{M} contains only the information embedded in the density at the positions of the peaks, not in the amount of liquid surrounding the peaks. When the polar proton and the acceptor atom of the dimer separate from each other, it creates a situation qualitatively similar to that in Figure 4 for two ions. In an ion pair, the peaks that exist at *r* > *σ*_{1} + *σ*_{2} are replaced by a continuum of saddle points at *r* < *σ*_{1} + *σ*_{2}. A molecular solute lacks the symmetry of the ion pair; therefore, the peaks that exist at *r* > *σ*_{H′}+ *σ*_{A′}in the dimer do not completely disappear at *r* < *σ*_{H′}+ *σ*_{A′}. However, the spatial extent of the peaks at *r* < *σ*_{H′}+ *σ*_{A′} still generates a repulsive force that complements the repulsive force calculated from the peaks densities. A formal way of taking these effects into account is through a suitable treatment of the local volume
. This problem is simplified here, and a correction, *δ***F**_{s,}* _{i}*, to

(34)

where **r**_{i}* _{j}* =

The total intermolecular force and associated potential can now be calculated. The resulting force, **F*** _{i}*, on atom

Molecular interactions in solution are controlled by bulk electrostatic forces and by solvent-induced forces. A complete description of solvent effects requires both to be accounted for. Efforts in developing models for use in computer simulations of macromolecules have focused almost exclusively on electrostatic forces. In contrast, the forces exerted by the microscopic nature of the structured solvent at the solute/liquid interface have been traditionally neglected. In this paper, a model has been developed to account for both effects, thus improving upon continuum representations.

The model of electrostatics developed here was based on the dielectric theory of polar and polarizable fluids. Both the polarizability and the permanent dipole moment of the water molecules were accounted for. The reaction field was explicitly introduced into the internal and directing electric fields. The model captures the saturation effects of the field and the spatial dependence of the liquid density. The density was represented by a barometric law with an effective Lennard-Jones-type solute–liquid interaction potential *V*_{eff}. A closed equation for the effective, position-dependent static dielectric permittivity, *ε*, was derived for arbitrary solutes and solved numerically for ions and ion pairs. Experimental free energies of hydration were used to calibrate the ion–liquid interaction potential. The potential thus derived was then used to calculate the free energy of hydration of ion pairs as a function of the interionic distance. The forces exerted by the structure of the liquid on ion pairs were then calculated. The results showed that the forces associated with the desolvation barriers and the solvent-separated minima observed in interionic potentials of mean force can be reproduced with a suitable choice of the effective potential parameters.

The results summarized above were obtained by numerical integration of eqs 6, 12, and 17. However, the numerical solution of the forces may be too demanding in biological applications of macromolecular simulations, and thus a simplified model for bulk electrostatic and liquid-structure forces was proposed. The electrostatic model has been developed previously, and its main features were reviewed. The model for solvent-induced forces was introduced here and analyzed in detail. The most important characteristic of the model is its non-pairwise nature, typical of the short-range liquid structure forces. The non-pairwise additivity of the liquid forces arises from pairwise forces between the solute and the peaks. The out-of-the-solute centers then exert local forces and torques on the solute. The solute–liquid interaction potential was optimized to provide a good representation of the liquid structure surrounding single charged and net-neutral polar amino acids. The model reproduces the spatial distribution of peaks obtained from atomistic dynamics simulations. Only peaks consistent with water molecules in direct H-bond interactions with acceptor and donor groups were considered. The fields created by these groups are strong enough to control the structure and dynamics of the surrounding water. Therefore, the effects of water–water correlations could be neglected, to a reasonable extent, and their effect on the local density could be captured in an effective water–solute interaction potential. The model is not expected to provide a good representation of liquid structure around nonpolar solutes, unless a more explicit treatment of water–water correlation is introduced. To partially compensate for this approximation in the limit of nonpolar molecules, a force obtained from a nonpolar solvation term was used.

Improvements of the model may follow two directions: to account for dipole–dipole correlations in the calculation of the electric polarization, hence, in the effective dielectric permittivity of eq 6; and to account for water–water correlation in an expression of the direct-pair correlation function, hence, in the calculation of the liquid density of eq 9. A Kirkwood–Fröhlich approximation^{29}^{,}^{110}^{,}^{111} may be the first correction to the Onsager model used here. The effects of dipole–dipole short-range interactions in the saturation effects of the field have been discussed.^{112}^{–}^{114} These corrections show that the onset of saturation starts at smaller values of the electric fields.^{112} This implies that *ε* reaches its bulk static value farther away from the source when compared to results from the Onsager approximation. The structured liquid in the vicinity of a solute may reinforce such effects. This illustrates a more general problem of molecular electrostatics, namely, the sensitivity of the results to the treatment of the solute/liquid interface. Reaction field effects and electrostriction tend to increase the values of the dielectric permittivity close to the solute,^{44}^{,}^{45}^{,}^{51}^{,}^{115} but bulk dipole–dipole correlation tends to decrease it,^{112} while correlations closer to the source may decrease it further. The effect of the liquid density is yet another factor to be considered, as discussed here and elsewhere.^{57} The coexistence of these competing effects makes it difficult to predict the correct behavior of *ε* near a solute, and the distance at which bulk liquid dielectric properties are recovered. This topic can be studied with atomistic dynamics simulations of the liquid in the proximity of a solute.^{116} To compare the results with the prediction of the model developed here, a polarizable force field would be required.^{117}^{,}^{118} On the other hand, experimental data provide some insight into the local electric field in proteins and protein–liquid interfaces.^{119}^{,}^{120} In particular, p*K* shifts of ionizable groups are sensitive to local protein electrostatics and can be used as benchmarks to validate and optimize solvent models.^{121} Data on p*K* shifts have been used^{44}^{,}^{58} to guide the optimization of the screening functions of the SCP–ISM reviewed in section IV.1.

The author thanks Peter Steinbach for reading the manuscript and for valuable discussions. This study was supported by the Intramural Research Program of the NIH, U.S. Department of Health and Human Services, Bethesda, Maryland.

1. Timasheff SM. Annu Rev Biophys Biomol Struct. 1993;22:67. [PubMed]

2. Yancey PH, Clark ME, Hand SC. Science. 1982;217:1214. [PubMed]

3. Reichardt C. Solvents and Solvent Effects in Organic Chemistry. 3. Wiley-VCH; Weinheim, Germany: 2002.

4. Karplus K. Acc Chem Res. 2002;35:321. [PubMed]

5. Karplus M, McCammon JA. Nat Struct Biol. 2002;9:646. [PubMed]

6. Makhatadze GI, Privalov PL. J Mol Biol. 1993;232:639. [PubMed]

7. Privalov PL, Makhatadze GI. J Mol Biol. 1993;232:660. [PubMed]

8. Carrion-Vazquez M, Oberhauser AF, Fowler SB, Marszalek PE, Broedel SE, Clarke J, Fernandez JM. Proc Natl Acad Sci USA. 1999;96:3694. [PubMed]

9. Cui Y, Bustamante C. Proc Natl Acad Sci USA. 2000;97:127. [PubMed]

10. Florin EL, Moy VT, Gaub HE. Science. 1994;264:415. [PubMed]

11. Dudko OK, Hummer G, Szabo A. Phys Rev Lett. 2006;96:108101. [PubMed]

12. Parsegian VA, Rand RP, Fuller NL, Rau DC. Methods Enzymol. 1986;127:400. [PubMed]

13. Scatchard G. J Am Chem Soc. 1946;68:2315.

14. Rau DC, Parsegian VA. Science. 1990;249:1278. [PubMed]

15. Stanley C, Rau DC. Biophys J. 2006;91:912. [PubMed]

16. Bennett CH. J Comput Phys. 1976;22:245.

17. Jarzynski C. Phys Rev Lett. 1997;78:2690.

18. Bolhuis PG, Dellago C, Chandler D. Faraday Discuss. 1998;110:421.

19. Hassan SA. J Phys Chem B. 2005;109:21989. [PMC free article] [PubMed]

20. Ben-Naim A. J Phys Chem. 1990;94:6893.

21. Bruge F, Fornilli SL, Malenkov GG, Palma-Vittorelli MB, Palma MU. Chem Phys Lett. 1996;254:283.

22. Durell SR, Brooks BR, Ben-Naim A. J Phys Chem. 1994;98:2198.

23. Bruge F, Cottone G, Fornili SL. Chem Phys Lett. 1995;235:402.

24. Martorana V, Corongiu G, Palma MU. Chem Phys Lett. 1996;254:292.

25. Ben-Naim A, Ting KL, Jernigan RL. Biopolymers. 1990;29:901. [PubMed]

26. Ben-Naim A. Biopolymers. 1990;29:567. [PubMed]

27. Ben-Naim A. J Chem Phys. 1990;93:8196.

28. Bottcher CJF. Theory of Electric Polarisation. Elsevier; Amsterdam, The Netherlands: 1952.

29. Madden P, Kivelson D. Adv Chem Phys. 1984;56:467.

30. Hansen J-P, McDonald IR. Theory of Simple Liquids. 2. Elsevier; London: 1986.

31. Evans R. Adv Phys. 1979;28:143.

32. Lebowitz JL, Percus JK. J Math Phys. 1963;4:116.

33. Chandler D, McCoy JD, Singer SJ. J Chem Phys. 1986;85:5971.

34. Ramirez R, Gebauer R, Mareschal M, Borgis D. Phys Rev E. 2002;66:31206.

35. Imai T, Kovalenko A, Hirata F. Chem Phys Lett. 2004;395:1.

36. Lum K, Chandler D, Weeks JD. J Phys Chem B. 1999;103:4570.

37. TenWolde PR, Chandler D. Proc Natl Acad Sci USA. 2002;99:6539. [PubMed]

38. Beglov D, Roux B. J Chem Phys. 1995;103:360.

39. Beglov D, Roux B. J Phys Chem B. 1997;101:7821.

40. Gilson MK, Davis ME, Luty BA, McCammon JA. J Phys Chem. 1993;97:3591.

41. Honig B, Nicholls A. Science. 1995;268:1144. [PubMed]

42. Warshel A, Russell ST. Q Rev Biophys. 1984;17:283. [PubMed]

43. Warshel A, Aqvist J. Annu Rev Biophys Biophys Chem. 1991;20:267. [PubMed]

44. Hassan SA, Guarnieri F, Mehler EL. J Phys Chem B. 2000;104:6478.

45. Hassan SA, Mehler EL. Int J Quantum Chem. 2005;102:986.

46. Tomasi J, Persico M. Chem Rev. 1994;94:2027.

47. Tomasi J, Mennucci B, Cammi R. Chem Rev. 2005;105:2999. [PubMed]

48. Maroncelli M. J Mol Liq. 1993;57:1.

49. Powles JG, Williams ML, Evans WAB. Mol Phys. 1989;66:1107.

50. Bertolini D, Tani A. Mol Phys. 1992;75:1065.

51. Ehrenson S. J Comput Chem. 1989;10:77.

52. Onsager L. J Am Chem Soc. 1936;58:1486.

53. Dogonadze RR, Kornyshev AA. J Chem Soc, Faraday Trans. 1974;70:1121.

54. Bopp PA, Kornyshev AA, Sutmann G. Phys Rev Lett. 1996;76:1280. [PubMed]

55. Kornyshev AA. Electrochim Acta. 1981;26:1.

56. Nienhuis G, Deutch JM. J Chem Phys. 1971;55:4213.

57. Chandra A, Bagchi B. J Phys Chem. 1989;93:6996.

58. Mehler EL. The Lorentz-Debye-Sack Theory and Dielectric Screening of Electrostatic Effects in Proteins and Nucleic Acids. In: Murray JS, Sen K, editors. Molecular Electrostatic Potential: Concepts and Applications. Elsevier Science; Amsterdam, The Netherlands: 1996. pp. 3–371.

59. Debye P. Polar Molecules. Dover; New York: 1929.

60. Weast RC, editor. CRC Handbook of Chemistry and Physics. Vol. 59. CRC Press; Boca Raton, FL: 1978.

61. Allen MP, Tildesley DJ. Computer Simulation of Liquids. Clarendon Press; Oxford, UK: 1987.

62. Broderix K, Bhattacharya KK, Cavagna A, Zippelius A, Giardina I. Phys Rev Lett. 2000;85:5360. [PubMed]

63. Marcus Y. Chem Rev. 1988;88:1475.

64. Conway BE. J Solution Chem. 1999;28:163.

65. Narten AH. J Phys Chem. 1970;74:765.

66. Ohtomo N, Arakawa K. Bull Chem Soc Jpn. 1979;52:2755.

67. Ohtomo N, Arakawa K. Bull Chem Soc Jpn. 1980;53:1789.

68. Jackson JD. Classical Electrodynamics. Wiley; New York: 1975.

69. Pierotti RA. Chem Rev. 1976;76:717.

70. Ashbaugh HS, Pratt LR. Rev Mod Phys. 2006;78:159.

71. Noyes RM. J Chem Phys. 1962;84:513.

72. Noyes RM. J Am Chem Soc. 1964;86:971.

73. Conway BE. J Solution Chem. 1978;7:721.

74. Tissandier MD, Cowen KA, Feng WY, Gundlach E, Cohen MH, Earhart AD, Coe JV. J Phys Chem A. 1998;102:7787.

75. Pliego JR, Riveros JM. Chem Phys Lett. 2000;332:597.

76. Marcus Y. Pure Appl Chem. 1987;59:1093.

77. Aqvist J. J Phys Chem. 1990;94:8021.

78. Chan DYC, Mirtchell DJ, Ninham BW. J Chem Phys. 1979;70:2946.

79. Pollock EL, Alder BJ, Pratt LR. Proc Natl Acad Sci USA. 1980;77:49. [PubMed]

80. Pollock EL, Alder BJ, Pratt LR. Chem Phys Lett. 1981;78:201.

81. Bucher M, Porter TL. J Phys Chem. 1986;90:3406.

82. Born M. Z Phys. 1920;1:45.

83. Hassan SA, Mehler EL. Proteins: Struct, Funct, Genet. 2002;47:45. [PubMed]

84. Ehrenson S. J Phys Chem. 1987;91:1868.

85. Latimer WM, Rodebush WH. J Am Chem Soc. 1920;42:1419.

86. Rashin AA, Honig B. J Phys Chem. 1985;89:5588.

87. Block H, Walker SM. Chem Phys Lett. 1973;19:363.

88. Davis ME, McCammon JA. J Comput Chem. 1991;12:909.

89. Abraham MH, Liszi J, Meszaros L. J Chem Phys. 1979;70:2491.

90. Hassan SA. J Phys Chem B. 2004;108:19501.

91. Pettitt BM, Rossky PJ. J Chem Phys. 1986;84:5836.

92. Guardia E, Rey R, Padro JA. Chem Phys. 1991;155:187.

93. Dang LX, Pettitt BM, Rossky PJ. J Chem Phys. 1992;96:4046.

94. Hassan SA, Mehler EL, Zhang D, Weinstein H. Proteins: Struct, Funct, Genet. 2003;51:109. [PubMed]

95. Mehler EL. J Phys Chem. 1996;100:16006.

96. Mehler EL, Eichele E. Biochemistry. 1984;23:3887.

97. Voges D, Karshikoff A. J Chem Phys. 1998;108:2219.

98. Song X. J Chem Phys. 2002;116:9359.

99. Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M. J Comput Chem. 1983;4:187.

100. Hassan SA, Guarnieri F, Mehler EL. J Phys Chem B. 2000;104:6490.

101. Hassan SA, Mehler EL. Int J Quantum Chem. 2001;83:193.

102. Li X, Hassan SA, Mehler EL. Proteins: Struct, Funct, Genet. 2005;60:464. [PMC free article] [PubMed]

103. Kirkpatrick S, Gelatt CD, Jr, Vecchi MP. Science. 1983;220:671. [PubMed]

104. Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. J Chem Phys. 1953;21:1087.

105. Finney JL, Bowron DT. Philos Trans R Soc London, Ser A. 2005;363:469.

106. Soper AK. Chem Phys. 1996;202:295.

107. Tavernier HL, Fayer MD. J Chem Phys. 2001;114:4552.

108. Swanson JMJ, Mongan J, McCammon JA. J Phys Chem B. 2005;109:14769. [PubMed]

109. MacKerell AD, Jr, Bashford D, Bellott M, Dunbrack RL, Jr, Evanseck JD, Field MJ, Fischer S, Gao J, Guo H, Ha S, Joseph-McCarthy D, Kuchnir L, Kuczera K, Lau FTK, Mattos C, Michnick S, Ngo T, Nguyen DT, Prodhom B, Reiher IWE, Roux B, Schlenkrich M, Smith JC, Stote R, Straub J, Watanabe M, Wiorkiewicz-Kuczera J, Yin D, Karplus M. J Phys Chem B. 1998;102:3586.

110. Kirkwood JG. J Chem Phys. 1939;7:911.

111. Frohlich H. Theory of Dielectrics. Clarendon; Oxford, UK: 1958.

112. Booth F. J Chem Phys. 1951;19:391. errata, ibidem 19.

113. Grahame DC. J Chem Phys. 1950;18:903.

114. Laidler KJ, Pegis C. Proc R Soc London. 1957;241:80.

115. Carnie SL, Stell G. J Chem Phys. 1982;77:1017.

116. Lamm G, Pack GR. J Phys Chem B. 1997;101:959.

117. Lamoureux G, Harder E, Vorobyov IV, Roux B, MacKerell AD. Chem Phys Lett. 2005;418:241.

118. Rick SW, Stuart SJ, Berne BJ. J Chem Phys. 1994;101:6141.

119. Cohen BE, McAnaney TB, Park ES, Jan YN, Boxer SG, Jan LY. Science. 2002;296:1700. [PubMed]

120. Mehler EL, Fuxreiter M, Simon I, Garcia-Moreno EB. Proteins: Struct, Funct, Genet. 2002;48:283. [PubMed]

121. Schutz CN, Warshel A. Proteins: Struct, Funct, Genet. 2001;44:400. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's Canada Institute for Scientific and Technical Information in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |