Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Phys Chem B. Author manuscript; available in PMC 2008 July 15.
Published in final edited form as:
PMCID: PMC2467438

Liquid-Structure Forces and Electrostatic Modulation of Biomolecular Interactions in Solution


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.

I. Introduction

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 molecules3 to the structure, dynamics, and thermodynamics of macromolecules.47 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.810 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 stress12,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.1618 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,2124 and their role in protein–ligand interactions and protein folding has been discussed.20,2527

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 media28,29 and theories of liquid structure.3033 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,3439 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 numerically40,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.

II. Electrostatics and Liquid Structure

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.4850 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 r can be written as28,51


where the number density ρ(r) represents the number of liquid molecules per unit volume at position r, and left angle bracketμright angle bracketis the statistical average of the permanent dipole moments of a liquid molecule at r. Ei(r) is the internal (or microscopic) field at the position of a liquid molecule. The orientations of the permanent dipoles are determined by the directing field52 Ed(r). Both internal and directing fields act on individual water molecules, but they are conceptually different from one another;52 expression for these fields will be given below. Equation 1 is the usual separation of the polarization into a component due to the induced dipoles, Pα (first term at the right-hand side), and a component due to the permanent dipoles, Pμ (second term). Thus, Pα is the only contribution to P in a nonpolar (μ = 0), polarizable (α ≠ 0) medium, while Pμ is the only contribution in a polar (μ ≠ 0), nonpolarizable (α = 0) medium. In both terms, ρ accounts for the number of liquid molecules locally contributing to PEquation 1 can be generalized to a multicomponent system.28

The polarization field is given, in general, by29 P(r) = ∫ χ(r,r′) · E(r′)dr′, 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,


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 nonlocality5357 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 r, and with a volume equal to that available to the molecule, is given by R(r) = f(r)μ*. In this expression, μ* is the total dipole moment of the molecule given by μ* = μ + αR(r), where the last term is the dipole induced by the reaction field itself. For a spherical cavity, f is a scalar quantity and, then, R(r) = μf(r)/[1 − αf(r)]. The scalar f can be calculated assuming that the cavity has a volume of 1/ρ(r) and is surrounded by an infinitely extended continuum with permittivity ε(r). This approximation has been discussed51,58 and can be considered a limit for the actual value of the local reaction field. In this case, f is given by28,52


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, Ei(r) – Ed(r) = left angle bracketR(r)right angle bracket. From this definition and the equation given above for R(r), the internal field can be written as


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


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


where E(r) [equivalent] |E(r)|, β [equivalent] 1/kBT (T is the absolute temperature; kB is the Boltzmann’s constant), L(x) = 1/tanh(x) − 1/x is the Langevin function which arises in the calculation28 of the canonical average left angle bracketμright angle bracket over the dipole orientations in eq 1.

In the Debye theory,59 Ei and Ed were not distinguished from one another. They were taken as equal to the Lorentz field,58 EL, which is defined as the electric field inside an ideal cavity surrounding the water molecule, and given by EL(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 corrections52 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 Ei and Ed, as in eqs 4 and 5. The LDS theory has been used in the past and forms the basis of the continuum electrostatics model44 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 atm60), an equation relating εand α is obtained,


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


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 structure30,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 G[ρ(r)] is defined and minimized with respect to ρ′(r). The function ρ′(r) that satisfies δG[ρ(r)]/δρ(r)=0 is the equilibrium, nonuniform spatial density distribution ρ(r) sought. This approach yields a formal expression for the liquid density in the form31


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′)dr′. 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


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


where ri = |rri|, and the sum runs over the N atoms of the solute. The distances σi mainly determine the positions of the density peaks around a molecule (N > 1) or the first peak of the radial distribution function, g(r) = ρ(r)/ρ0, of the liquid around a monovalent ion (N = 1). The energies ξi mainly determine the height of the density peaks. Comparison of eqs 9 and 10 indicates that Veff contains information on the exact solute–liquid potential energy, U, and on all of the intra-liquid interactions that can possibly be captured into the LJ functional form of eq 11 and its parameters. In particular, Veff contains the van der Waals solute–liquid interaction potential commonly used in classical mechanics force fields. It also contains information on the short-range electrostatics, particularly, that associated with the direct solute–water hydrogen-bond interactions in the hydration shells of polar and charged groups. In this paper, the choice of a LJ-type function as an effective solute–liquid interaction potential is based solely on practical considerations. LJ potentials have been used extensively in simulations of liquids,61 and their mathematical properties are well characterized.62 The only requirements for Veff are (i) to reproduce the short-range structure of the liquid surrounding a molecular solute and (ii) to provide for the short-range non-pairwise additivity of the resulting forces. These forces will operate in addition to the long-range electrostatic forces exerted by the essentially structureless bulk liquid. The adequacy of eqs 10 and 11 and the conditions to achieve these goals will be assessed and discussed in detail in the following sections.

III. Solvation of Alkali and Halide Ions

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.

III.1.Single Ions

For a monovalent ion, the parameter σ in eq 11 is separated into two parts, σ = RI + Rw, where RI is the radius of the ion in the liquid phase at infinite dilution, and Rw is an extension identified with the radius of a water molecule ([equivalent]1.38 Å). This definition is consistent with the way the ionic radii are calculated from experimental data.63,64 The radius RI 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 radii6567 given in Table 1. The reversible work of polarization of the dielectric medium is given by68

Experimental and Calculated Propertiesa of Alkali and Halide Ions in Waterb


where E0(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) = E0(r)/ε(r), with |E0(r)| = 1/r2 for the ion at r = 0, eq 12 yields


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 values60 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 ΔG0. The work of polarization ΔAe, calculated from eq 12, corresponds to a Helmholtz free energy of an isochoric thermodynamic process. The relation between both quantities is71,72


where ΔGnp is the nonpolar component of the free energy. This term can be decomposed as ΔGnp = ΔGcav + ΔGsr, where ΔGcav accounts for the reversible work of cavity formation, which can be described by scaled particle theories,69,70 and ΔGsr is a short-range solute–liquid interaction term. For charged species, ΔGnp makes a small contribution to ΔG0. Since the interest here is in the electrostatic component of ΔG0, the nonpolar term is approximated as ΔGnp ~ 1.325 kcal/mol for all of the ions.71,72 Partitioning the experimental ΔG0 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, ΔGH, which is used as a reference.72,73 Cluster–ion solvation data have recently been used74 to determine ΔGH; 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 term76 of ~1.9 kcal/mol has been added to the work of polarization ΔAe in eq 14 to relate it to the experimental quantity ΔGe0 [equivalent] ΔG0 − ΔGnp).

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 ∫r2g(r)dr, where the upper limit of integration is taken as Rm = RI + 2Rw, 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 Rm. 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 ΔAe as a function of the integration limit Ru in eq 13, showing that convergence is achieved only at Ru > ~30 Å. The contribution of the first hydration shell to ΔAe (denoted by ΔAe′ in Table 1) is similar in magnitude to the bulk contribution.

Figure 1
Permittivity profiles for a monovalent ion in water as a function of the distance to the center of the ion (eq 6): Na+ and Cl in a uniform (ρ(r) = ρ0) medium (dashed line); and Na+ (thin solid) and Cl (thick solid) in ...

III.2. Ion Pairs

III.2.1. Electrostatic Modulation

For an arbitrary static charge distribution, the relationship between the macroscopic and vacuum fields is given by68 E(r) = E0(r) −[nabla] P(r′)·[nabla]′|rr′|−1dr′, for |r| > 0. Solving this equation for E(r) (or the corresponding expression for time-dependent fields29) 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 E0(r) linearly related through E(r) = E0(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 earlier45,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, ΔAe(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,


The last two terms are self-energies; ΔAself,1(r) results from zeroing the charge on ion 2 (i.e., Cl in the Na+–Cl pair), while ΔAself,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, ΔAself 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 |ΔAself(r)| of an ion to decrease, as expected, because of the removal of polarizable medium near the charged particle. The term ΔAint(r) is calculated from eq 13 as the difference between the total ΔAe(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 ΔAe(r) with the distance (Figure 2A) is related to the interaction term ΔAint(r). The behavior of ΔAe(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 Fe,i = −dΔAe(r)/dri. The total electrostatic energy of the system is given by ET(r) = ΔAe(r) + ΔE0(r), where ΔE0(r) ) 1/r is the Coulomb energy in the vacuum.

Figure 2
(A) Electrostatic contribution of the free energy of hydration of ion pairs as a function of the interionic distance (eq 13): Na+–Cl(solid line); Na+–Na+ (thin dashed); and Cl–Cl (thick dashed). (B) Self-energy ...

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 RB is used; thus RB defines the lower limit of integration in eq 13. RB 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 ε(RB) for the five ions considered here. Following precedence, RB can be decomposed into two terms44,81,85,86 in the form RB = RI + δ, where δ is an extension of the ionic radius. Thus, RB 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 previously44,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, ΔAe ≈ (1/ε 0 − 1)/RB for a monovalent ion, with RB 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,8789 The model presented here provides a simple recipe to avoid such discontinuities in a physically meaningful way.

III.2.2. Liquid-Structure Forces

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 dF(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 by20,22,27


where U(r, r′) is the interaction potential energy between the atom and a liquid molecule, and dn(r′) = ρ(r′)dν′is the number of liquid molecules at r′. Here, the potential energy is given by U(r, r′) ≈ Veff(|rr′|) [cf. eq 11], and the density at r′is given by eq 10. The total liquid-structure force Fs,i(r) on atom i is then calculated by integrating dF(r, r′) over r′,


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, dF(r, r′) is calculated at each point r′ and then integrated in space to yield the total forces Fs,i(r) = Fs,i(r)zon ion i. From eq 17, the force Fs,2(r) on ion 2 is given by


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 Fs,1(r); Fs,2(r) = −Fs,1(r) in equilibrium). In eq 18, Fs,2 < 0 implies that the solvent induces an attractive force between the charges, while Fs,2 > 0 implies a repulsive force. The potential associated with Fs,2(r), that is, a potential of mean liquid-structure force type, is given by90


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 Veff 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 Fs,2(r) and [var phi]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 = rm > σ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 = rM < σ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 rm and the maximum at rM 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)9193 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].

Figure 3
Liquid-structure force (eq 18) (thick line; in kcal/mol/Å) and associated potential (eq 19) (thin; in kcal/mol) for the Na+–Clpair as a function of the interionic distance. The relative positions of the solvation spheres of the ...

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.

IV. Molecular Solutes

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, Fi, on an atom i of a molecule is divided into four terms


where Fc,i and FvdW,i are the total bare Coulomb and van der Waals forces on atom i due to all other atoms of the solute, respectively. Fe,i is the total electrostatic force exerted by the essentially structureless liquid of the bulk phase, and Fs,i is the total force exerted by the structured liquid surrounding the molecule. A simplified, pairwise model for Fe,i has been reported,94 and a summary is given below for completeness; the connection with the electrostatic model introduced in section II will be discussed. A simplified, non-pairwise model of Fs,i will be derived and discussed in detail in the following sections.

IV.1. Bulk Electrostatic Forces

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 (ΔAint) and the self-energy (ΔAself) 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.

IV.1.1. Single Charge Q

In this case, the electric potential in a polar/polarizable liquid is given by [var phi](r) = [var phi]0(r)/D(r) where [var phi] 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) = −[nabla][var phi](r). For a point charge, the relationship between both quantities is95 ε = D/[1 + (r/D)dD/dr]. 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,


where RB is the Born radius discussed in section III.2.1, which includes the saturation effect of the field. In this form, ΔAself can be evaluated with little computational cost provided that the functional form of D(r) is known and a prescription to determine RB is given. A sigmoidal function of D(r) that reproduces ε(r) in polar/polarizable liquids has been reported58,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 RB follows the same idea discussed in section II for ionic sources, that is, a charge radius RQ plus an extension δ (RB = RQ + δ).

IV.1.2. Two Interacting Charges Q1 and Q2

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



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 ΔAself,i with r is introduced through the effective Born radius RB,i(r) as a linear combination of two Born radii44,94



where R1B and R2B are the Born radii RB of eq 21 for each charge; ξ12 and ξ21 are the fractions of the solvent-accessible surface area of Q1 and Q2, respectively. The factors 1 − ξ12 and 1 − ξ21 are the fractions of the solvent-excluded surface area of each charge. The quantity R12 accounts for the effects on Q1 of the absence of polarizable medium in the region occupied by Q2; R21 accounts for the same effect on Q2 due to the presence of Q1. The linear combination of eq 23a,b is the simplest way to account for the dependence of ΔAself,i(r) versus r [cf. eq 22b] observed in the exact numerical solutions plotted in Figure 2B. The radii R12 and R21 control the hydration energy of the pair at short interparticle distances, and Figure 2B suggests that R12 > R1w and R21 > R2w. Therefore, in analogy with R1B and R2B, these Born radii are given by R12 = R1w + δ12 and R21 = R2w + δ21 (with δ > 0 in both cases44). This model reproduces the qualitative behavior of the self-energies [cf. Figure 2B] and the interaction term [cf. Figure 2A]. Quantitative agreement is attained with a suitable choice of parameters λ’s, δ12, and δ21; the fractions ξ12 and ξ21 can be calculated analytically.44,94

IV.1.3. N Interacting Charges {Qi}i = 1,N

For a molecular solute composed of N atoms, eq 22a,b can be generalized as44



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, RB,i is given by a generalization94 of eq 23, which accounts for the N − 1 atoms surrounding atom i; RB,i′is obtained as a particular case.44 For molecular solutes, the equivalent to RI or RQ is the covalent radius, RC, of the atom in the molecule.44 The total free energy of hydration of the molecule is given by eq 14. The nonpolar term scales roughly as the solvent-accessible surface area (SASA) of the molecule (i.e., ΔGnp = a + b SASA, where a and b are constants).

From eq 24, the total electrostatic energy ET of the hydrated molecule is given by44


As defined in eq 20, the electrostatic component of the solvent force on an atom i of the solute is then calculated as Fe,i = −[nabla]ETFc,i. A general expression has been derived in ref 94 and is given by


with the interaction and self-energy components given by



where rij = rirj and rij [equivalent] |rij|. Equations 2527 have been implemented into the CHARMM molecular dynamics program.99 In combination with an empirical correction for H-bond interactions in solution100 (see section IV.2.4), these equations define the SCP-based implicit solvent model (SCP–ISM). The model has been parametrized for peptides and proteins and used in previous studies, including structure calculations and dynamic simulations.83,94,101,102

IV.2. Liquid-Structure Forces

To derive a computationally efficient model for Fs,i, a strategy similar to the one summarized above for Fe,i is followed.

IV.2.1. Two Interacting Charges Q1 and Q2

To gain insight into the dependence of Fs,2(r) on r shown in Figure 3, the integrand of eq 17 is analyzed. A force density function (η,z;r) for a fixed r can be defined from eq 17 as Fs,2(r)=(η,z;r)dηdz. Figure 4 shows a gray-scale plot of (η,z;r) 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 (Fs,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 (Fs,2 = 0). This follows from the functional form of Veff 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 (Fs,2 = 0), but with the density distribution controlled by both ions.

Figure 4
Gray-scale representation of the liquid-structure force density function (η,z;r) (integrand of eq 18) for Na+–Cl, as a function of the interionic distance. Black: regions of the liquid inducing a repulsive force between ...
Figure 5
(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, (η ...

The main characteristics of Fs,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, |Fs,2(r)| becomes smaller for all r, as expected, but both the minimum at rm and the maximum at rM remain well defined. This suggests that the modulation of Fs,2(r) versus r originates mainly in regions of high liquid density. Because of the functional form of Veff [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, (η,z;r), discussed above suggest a simple model for liquid-structure forces in a system composed of two charges. In this model, the total force Fs,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 Fs,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, (η,z), on an isolated ion shown in Figure 5B and by the behavior of (η,z;r) discussed in section III.2.4 and illustrated in Figure 4, on ion pairs. An external shell of thickness de is formed by two hemispherical shells; Ve+ exerts a repulsive force, and Ve exerts an attractive force. Similarly, the internal shell of thickness di is formed by the hemispherical shells Vi+ and Vi (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 Ve+, Ve, Vi+, and Vi 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 Ve and Vi+ change according to the extent of overlap of the solvation spheres (of radii σ1 and σ2). These overlapping regions are labeled Ve and Vi. If λe > 0 and λi > 0 are parameters representing the force per unit volume exerted by each shell, then


where the explicit dependence on r was indicated in Ve(r) and Vi+(r). Since the volumes of the hemispherical shells are Ve+ = Ve + Ve and Vi = Vi + Vi+, the force is given by


where the volumes Ve(r) and Vi(r) can be calculated as the difference of two volumes, Ve(r) = V1e(r) − V12(r) and Vi(r) = V12(r) − V1i(r), where V1e(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 + de). V12(r) is the overlapping volume of the solvation spheres of ion 1 and ion 2, and V1i(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 σ2di). These volumes can then be written as


where the index l is e, 2, or i, V1=4πσ12/3, and θ(x) is the Heaviside function defining the values of r for which overlap between the spheres exist. The functions x1l and y1l are given by x1l = (r1 + r2) − r and y1l = r − (r2r1), with r1 = σ1 and r2 = σ2 + de (for l = e), r2 = σ2 (l = 2), and r2 = σ2di (l = i). The volumes V1l of the overlapping spheres are continuous functions of r given by 3V1l(r) = πa2(3r2a) + πb2(3r1b) [i.e., the overlapping volume of two spheres of radii r1 and r2 with center-to-center distance r (see inset of Figure 5)]; the distances a and b are given by a=r2(r2+r22r12)/2r and b=r1(r2+r12r22)/2r, with r1 and r2 defined as x1l and y1l. Introducing these expressions into eq 29 yields Fs,2(r) as a continuous function of r that can be evaluated rapidly. The precise form of Fs,2(r) depends on the values of the force density parameters λe and λi, and the thicknesses, de and di, of the outer and inner shells, respectively. Figure 6 shows Fs,2(r) for two arbitrary sets of parameters (with λe < λi and de > di), which can be compared to the plot in Figure 3. The minimum at rm is due to the forces exerted by the peaks at r > σ1 + σ2, which are not discussed in this section (see section IV.2.4). The overall shape of the curve is reproduced but with 4 orders of magnitude faster computation than the numerical integration of eq 17. These practical considerations justify introducing the set of four parameters (λe, λi, de, and di) in the model. These parameters control the location (rR), the height (ΔFR), and the width (ΔrR) of the repulsive force (the precise form of Fs,2 at a short interionic distance is of no practical interest because of the interionic core repulsion). The basis of this model is the competition of forces exerted by the four hemispherical shells around the ions. The strengths of the forces exerted by each shell and the regions where each shell operates are controlled by λ and d, respectively.

Figure 6
Liquid-structure force between two ions at r < σ1 + σ2, calculated with the double-shell model (eq 29) as a function of the interionic distance (force in kcal/mol/Å). The parameters of the model control the height (Δ ...

IV.2.2. N Interacting Charges {Qi}i=1,N

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.

IV.2.2.1. Single Amino Acids

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/cm3). 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)dt, with N(r, t)(r = Σ δ−ri(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 {rf}f=1,n, consistent with direct H-bond interactions with polar protons or acceptor atoms of the solute’s functional group.19 The set of values of g(r) at the positions of the peaks are denoted by {gf}f=1,n, with gf [equivalent] g(rf).

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 < gf < 3.4). In Lys+ (B), the peaks are in a tetrahedral arrangement with respect to the central nitrogen (2.8 < gf < 3.3). In His+ (C), the peaks are in the plane of the ring (gf 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 < gf < 3.4). For the polar amino acids, the values of gf 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 < gf < 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 His0 (G) ring (1.8 < gf < 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 < gf < 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.

Figure 7
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 Veff(r) of eq 11. In principle, each atom i of the solute can be assigned a parameter, ξi and σi, to be optimized. Here, this general optimization is simplified by defining a subset of atom types: polar hydrogen (ξH and σH) and acceptor (ξA and σA) atoms in neutral functional groups, polar hydrogen (ξH′ and σH′) and acceptor (ξA′ and σA′) atoms in charged functional groups, and all other atoms (ξ0 and σ0). With this definition, there are four or six parameters per molecule, depending on the amino acid type. The optimization then consists of finding the values of the ξ’s and σ’s that reproduced the positions {rf} and densities {gf} of the peaks obtained from the MD simulation. The solute–liquid interaction potential is optimized for each molecule separately. The optimization is carried out using simulated annealing Monte Carlo in the space of the parameters.103 A function S is defined as


where rp and gp are the position and density of a peak of ρ(r) calculated from eq 13, using the potential Veff(r) of eq 16, which contains the set {ξ, σ} to be optimized. The peaks are formally defined as the set of m points {rp}p=1,m, where ρ(r) has local maxima. The basic assumption is that the static number density ρ(r) defined in eq 11 can be equated to the number density δN(r)/δV obtained as a time average from the liquid dynamics. The function S is minimized using a Boltzmann-like factor, exp(−S/T), with the dimensionless variable T decreasing in a logarithmic schedule. Trial moves are accepted or rejected according to a Metropolis criterion.104 For arbitrary values of {ξ, σ}, the number, m, of peaks of ρ(r) is, in general, different than n, so the following criterion is used. If m > n in a trial move, the subset of n peaks used in eq 31 is that which minimizes the distance defined by the square root of the first sum in eq 31; if m < n, then the trial move is rejected. To identify a peak of ρ(r), a two-step focusing protocol is used for computational efficiency. Thus, spherical grids of radii rg and cell size (Δr, Δθ, Δ[var phi]) (in spherical coordinates) centered at each of the polar hydrogen and acceptor atoms of the solute are defined. Local maxima of g(r) are found within each grid by directly comparing the value of g(rijk) at a point (i, j, k)[equivalent] (ri, θj, [var phi]k) of the grid with the values of g(ri′,j′,k) at each of the 26 neighboring points (i′, j′, k′) = (i ± a, j ± b, k ± c), with a, b, c equal to −1, 0, +1 (not all 0). A local maximum at (i, j, k) is defined by the condition g(rijk) > g(ri′,j′,k) for all (ijk′). Once a maximum of g(r) is found at rijk in this grid, a new fine-grained homogeneous grid of cells size (Δr′, Δθ′, Δ[var phi]′) is built around rijk, with ri−1 < r < ri+1, θj−1 < θ < θj+1, and [var phi]k−1 < [var phi] < [var phi]k+1. A search of a local maximum at each point (i, j, k) in this new grid is repeated as before. Once a new maximum is found at a point rijk in the finer grid, the peak is defined by its coordinates rp [equivalent] rijk and its density gp [equivalent] g(rijk); these are the values used in eq 31. The uncertainty in the locations of the peaks calculated from the MD simulations is 0.65 Å (see ref 19 for details), while the partitions of the spherical grids are chosen such that the uncertainty in the location of a peak is 0.05 Å or less. If two peaks at positions rp′and rp″are separated by a distance |rp′rp″ | < 1.38 Å, then they are replaced by a single peak located at rp = (gp′rp′+ gp″ rp″)/(gp′+ gp″), with density gp = (gp′+ gp″)/2.

Figure 7 shows the location of the peaks after optimization of the potential Veff (white dots). In all cases, the distributions of {rp} closely follow those of {rf}, although some deviations are observed in Asp (in particular, the distances of {rp} to the oxygen atoms in the –COO group tend to be shorter than those for {rf}). 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 {rp} with respect to {rf}, averaged over the eight molecules, is 0.8 Å. The values of {gp} are also well reproduced, yielding an average rmsd with respect to {gf} of 0.4. These results show that a barometric law with an optimized Lennard-Jones-type potential, Veff(r), yields satisfactory results despite the relatively small set of parameters used.

IV.2.2.2. Amino Acid Pairs

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), [var phi](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.105107 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, rcm, the desolvation barrier or transition state, rts, and the solvent-separated minimum, rss, along with the corresponding values of the potentials, that is, [var phi]cm [equivalent] [var phi](rcm), [var phi]ts [equivalent] [var phi](rts), and [var phi]ss [equivalent] [var phi](rss). 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, empirically100 accounts for the values of [var phi]cm obtained in atomistic dynamic simulations of the liquid.90 This empiricism is needed because, if the correct value of [var phi]cm is sought and only an electrostatic model is used, the correction for the missing SIF effects on [var phi]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 [var phi]cm at rcm. 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, Fs,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, Fs,i can be written as


where r is the position of atom i in the solute, n is the total number of peaks in the liquid, and N(rj) is the number of water molecules at the position rj of the j-th peak. The number of water molecules at rj is N(rj)=V(rj)ρ(rj), where ρ(rj) is the density at the peak [cf. eq 10], and V(rj) is a local volume which accounts for the spatial extent of the peak in the three-dimensional space. Accounting for V(rj) is conceptually necessary because the value of the density at a peak carries no information on the density distribution around the peak; two points of equal density might be surrounded by a different amount of liquid, owing to different packing of the solute atoms around each point. This information is automatically captured when the integral in eq 17 is carried out over the space but is removed when the integral is replaced by the sum over peaks in eq 32. To account for molecular packing in a computationally efficient way, a contact model, akin to that used in ref 94, may prove useful. For simplicity, the position dependence of the volume will be ignored, and V(rj)=V for all peaks. A correction, δFs,i, will be added later [cf. eq 34] to partially compensate for this approximation. Regardless of the functional form of V(rj), the total force Fs = ΣFs,i and torque τ = Σri [equivalent] Fs,i exerted by each peak is zero in equilibrium because the gradient of the potential at the position of the peaks is zero by definition.

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


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 Veff. For nonpolar solutes, the force, Fs,i, calculated from eq 32 is strictly zero, but it is different than zero if eq 17 is integrated numerically. Therefore, a nonpolar force, Fnp,i = −[nabla]ΔGnp, is added to Fs,i of eq 32 as a correction term [cf. section IV.1] to represent the forces in the limit of nonpolar solutes. For nonpolar solutes, Fnp,i is the only contribution to the total SIF (the simplest approximation to hydrophobic forces commonly used in molecular simulations). For polar and charged solutes, both Fnp,i and Fs,i contribute to the total SIF. The force given by eq 32 can be viewed as originating in the structure of the liquid, induced in the presence of polar and charged groups. If a more sophisticated treatment of c(2) in eq 9 was available, eq 32 could be extended to include all of the peaks surrounding the solute. Then, the force Fnp,i would be included in the total SIF, and there would be no need to resort to the empirical term ΔGnp to represent such forces.

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, V=1A˚3, 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 described90 (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 Fs(r) [equivalent] Fs(rR (thin line) and the associated potential [var phi]s(r) (thick line) for the dimer. R 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: Fs(r) = [Fs,A(r) − Fs,D(r)]/2, where Fs,A and Fs,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 Fs,i of eq 33 over all atoms i of the monomer. Fs(r) < 0 implies attraction between the monomers. The peaks associated with the long-range attractive force are marked with arrows in Figure 8D. These peaks merge at r = rm, yielding only one peak located between the proton and the acceptor, which is marked with an arrow in Figure 8C. The location of this peak is consistent with water molecules bridging the proton and the acceptor atom through H-bonds, as discussed,19 and is responsible for the minimum of Fs(r) at r = rm. The effect of these peaks is identical to the effects of those discussed qualitatively in section III for the case of ion pairs at r > σ1 + σ2 [cf. Figure 3]. In that case, these peaks produce the long-range attraction between the ions and are associated with the minimum of Fs(r) at rm in Figure 3. Dynamics simulations show that new peaks appear between the monomers as they dissociate beyond the solvent-separated minimum distance, which further modulates the intermolecular potential.19 The model developed here does not account for liquid–liquid interactions, and thus, no new maxima and minima exist for r > rm associated with the interactions between peaks.

Figure 8
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 9
Forces and potentials for the Asp–Arg+ dimer. (A) Total, non-pairwise intermolecular force Fs(r) (thick line) exerted by the peaks of liquid density (eq 33) and associated potential [var phi]s(r) (thin, as a function of the proton-acceptor ...

Figure 9A also shows a maximum of Fs(r) at r = rM. This maximum is similar to that in Figure 3 for the case of ion pairs. As discussed above, Fs(r) at rM 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 V(rj). This problem is simplified here, and a correction, δFs,i, to Fs,i is introduced, which enhances the repulsive forces around r = rM in a computationally efficient way. This correction is based on eq 29 and applies only to acceptor atoms entering the solvation sphere of a polar proton. Therefore, the liquid-structure force on atom i is composed of non-pairwise contributions, Fs,i, generated by the density peaks and given by eq 33 and of proton-acceptor pairwise contributions, δFs,i, given by eq 29. If i is an acceptor and j is a polar proton, eq 29 yields


where rij = rirj and rij [equivalent]| rij|; λe, λi, de, and di are the parameters for i; NH is the total number of polar protons in the solute; and the force on a proton j is given by eq 34 after changing rij to −rij and NH to NA (total number of acceptor atoms). Figure 9B shows the correction δFs(r) [equivalent] δFs(rR (thin line) and the associated potential δεs(r) (thick line) of the dimer. As in Figure 9A, δFs(r) = [δFs,A(r) − δFs,D(r)]/2, where δFs,A and δFs,D are the total force, δFs, on the acceptor and donor monomer, respectively, evaluated at their centers-of-mass. The most important feature displayed in Figure 9B is the added repulsive force around r = rM (parameters not optimized). This force complements the repulsive force of the peaks in Figure 4A; the form of δFs(r) at small values of r is of no practical interest.

The total intermolecular force and associated potential can now be calculated. The resulting force, Fi, on atom i is given by eq 20. Fc,i and FvdW,i are calculated with the all-atom representation109 of the CHARMM force field.99 Fe,i is given by eq 26. Fs,i is given by eq 33; a nonpolar force correction, Fnp,i, is added. δFs,i is given by eq 34. Figure 9C shows the intermolecular potential calculated from the total force F(r) = Fc(r) + FvdW(r) + Fe(r) + Fs(r) (thick line) and the PMF obtained from the MD simulations reported in ref 90 (thin). The main features of the PMF can be reproduced within the statistical errors, although rss is slightly shorter (~4.1 vs ~4.6 Å). This is due to the small number of parameters used to optimize Veff(r) for Asp-, which results in small deviations of the positions of the peaks with respect to the results from the simulation [cf.Figure 7].

V. Discussion and Conclusions

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 Veff. 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 approximation29,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.112114 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, pK 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 pK shifts have been used44,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.

References and Notes

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]