138

Roux

that ∂Uuv/∂λ in Eq. (10) plays the role of a generalized thermodynamic force similar to that of ∂U/∂xi in Eq. (8).

III. DECOMPOSITION OF THE FREE ENERGY

Intermolecular forces are dominated by short-range harsh repulsive interactions arising from Pauli’s exclusion principle, van der Waals attractive forces arising from quantum dispersion, and long-range electrostatic interactions arising from the nonuniform charge distribution. It is convenient to express the potential energy Uuv(X, Y) as a sum of electrostatic contributions and the remaining nonpolar (nonelectrostatic) contributions,

Uuv(X, Y) U (np)uv (X, Y) U (elec)uv (X, Y)

(11)

Although such a representation of the microscopic non-bonded interactions does not follow directly from a quantum mechanical description of the Born–Oppenheimer energy surface, it is commonly used in most force fields for computer simulations of biomolecules (e.g., AMBER [5], CHARMM [6], OPLS [7]). The separation of the solute–solvent interactions in Eq. (11) is useful for decomposing the reversible work that defines the function W(X). The total free energy of a solute in a fixed configuration X may be expressed rigorously as the reversible thermodynamic work needed to construct the system in a step- by-step process. In a first step, the nonpolar solute–solvent interactions are switched ‘‘on’’ in the absence of any solute–solvent electrostatic interactions; in a second step, the solute– solvent electrostatic interactions are switched ‘‘on’’ in the presence of the solute–solvent nonpolar interactions. The solute is kept in the fixed configuration X throughout the whole process, and the intramolecular potential energy does not vary during this process. By construction, the total PMF is

W(X) Uu(X) ∆W (np)(X) ∆W (elec)(X)	(12)
where the nonpolar solvation contribution is
exp{ ∆W (np)(X)/k BT} ∫ dY exp{ [Uvv(Y) U (np)uv (X, Y)]/k BT}	(13)
∫ dY exp{ Uvv(Y)/k BT}
and the electrostatic solvation contribution is

∆ ∫ dY exp{ [U (Y) U (np)(X, Y) U (elec)(X, Y)]/k T}

exp{ W (elec)(X)/k BT} vv uv uv B

∫ dY exp{ [Uvv(Y) U (np)uv (X, Y)]/k BT}

(14)

Combining Eqs. (12)–(14) yields Eq. (7) directly. Although such a free energy decomposition is path-dependent [8], it provides a useful and rigorous framework for understanding the nature of solvation and for constructing suitable approximations to the nonpolar and electrostatic free energy contributions.

In the following sections, we describe an implicit solvent model based on this free energy decomposition that is widely used in biophysics. It consists in representing the nonpolar free energy contributions on the basis of the solvent-accessible surface area

Implicit Solvent Models

139

(SASA), a concept introduced by Lee and Richards [9], and the electrostatic free energy contribution on the basis of the Poisson–Boltzmann (PB) equation of macroscopic electrostatics, an idea that goes back to Born [10], Debye and Hu¨ckel [11], Kirkwood [12], and Onsager [13]. The combination of these two approximations forms the SASA/PB implicit solvent model. In the next section we analyze the microscopic significance of the nonpolar and electrostatic free energy contributions and describe the SASA/PB implicit solvent model.

A. Nonpolar Free Energy Contribution

To clarify the significance of ∆W (np), let us first consider the special case of a nonpolar molecule solvated in liquid water. We assume that the electrostatic free energy contribution is negligible. Typically, the solute–solvent van der Waals dispersion interactions are relatively weak and the nonpolar free energy contribution is dominated by the reversible work needed to displace the solvent molecules to accommodate the short-range harsh repulsive solute–solvent interaction. For this reason, ∆W (np) is often refered to as the ‘‘free energy of cavity formation.’’ The reversible thermodynamic work corresponding to this process is positive and unfavorable. It gives rise to two aspects of the hydrophobic effect: hydrophobic solvation and hydrophobic interaction [14]. The former phenomenon is responsible for the poor solubility of nonpolar molecules in water; the latter accounts for the propensity of nonpolar molecules to cluster and form aggregates in water.

Modern understanding of the hydrophobic effect attributes it primarily to a decrease in the number of hydrogen bonds that can be achieved by the water molecules when they are near a nonpolar surface. This view is confirmed by computer simulations of nonpolar solutes in water [15]. To a first approximation, the magnitude of the free energy associated with the nonpolar contribution can thus be considered to be proportional to the number of solvent molecules in the first solvation shell. This idea leads to a convenient and attractive approximation that is used extensively in biophysical applications [9,16–18]. It consists in assuming that the nonpolar free energy contribution is directly related to the SASA [9],

∆W (np)(X) γv tot(X)

(15)

where γv has the dimension of a surface tension and tot(X) is the configuration-dependent SASA (note that both polar and nonpolar chemical groups must be included in the SASA for a correct estimate of ∆W (np)). As pointed out by Tanford [19], there should be a close relationship between the macroscopic oil–water surface tension, interfacial free energies, and the magnitude of the hydrophobic effect. However, in practical applications, the surface tension γv is usually adjusted empirically to reproduce the solvation free energy of

˚ 2

alkane molecules in water [18]. Its value is typically around 20–30 cal/(mol A ), whereas

˚ 2

the macroscopic oil–water surface tension is around 70 cal/(mol A ) [19]. The difference between the optimal parameter γv for alkanes and the true macroscopic surface tension for oil/water interfaces reflects the influence of the microscopic length scale and the crudeness of the SASA model. A simple statistical mechanical approach describing the free energy of inserting hard spheres in water, called scaled particle theory (SPT), provides an important conceptual basis for understanding some of the limitations of SASA models [20–22]. It is clear that the SASA does not provide an ultimate representation of the nonpolar contribution to the solvation free energy. Other theories based on cavity distributions in liquid water [23,24] and long-range perturbation of water structure near large

140

Roux

obstacles [25] are currently being explored. A quantitative description of the hydrophobic effect remains a central problem in theoretical chemical physics and biophysics.

B. Electrostatic Free Energy Contribution

The electrostatic free energy contribution in Eq. (14) may be expressed as a thermodynamic integration corresponding to a reversible process between two states of the system: no solute–solvent electrostatic interactions (λ 0) and full electrostatic solute–solvent interactions (λ 1). The electrostatic free energy has a particularly simple form if the thermodynamic parameter λ corresponds to a scaling of the solute charges, i.e., U (elec)uv (X, Y; λ) λU (elec)uv (X, Y), and the coupling is linear,

∆W elec(X) 1 dλ U (elec)uv (λ) (16) 0

For this reason, the quantity ∆W (elec)(X) is often called the ‘‘charging free energy.’’ If one assumes that the solvent responds linearly to the charge of the solute, then U (elec)uv (λ) is proportional to λ and the charging free energy can be written as

∆W elec(X) 0	1	dλ qiΦrf(xi; λ)	1	qiΦrf(xi; λ 1)	(17)
			2
		i		i

where Φrf(xi; λ 1) is the solvent field acting on the ith solute atomic charge located at position xi in reaction to the presence of all the solute charges (in the following, the coupling parameter λ will be omitted for the sake of simplicity). The ‘‘reaction field’’ is thus the electrostatic potential exerted on the solute by the solvent that it has polarized. The assumption of linear response implies that ∆W elec (1/2) U (elec)uv , a relationship, that is often observed in calculations based on simulations with explicit solvent [26,27]. The factor 1/2 is a characteristic signature of linear solvent response.

The dominant effects giving rise to the charging free energy are often modeled on the basis of classical continuum electrostatics. This approximation, in which the polar solvent is represented as a structureless continuum dielectric medium, was originally pioneered by Born in 1920 to calculate the hydration free energy of spherical ions [10]. It was later extended by Kirkwood [12] and Onsager [13] for the treatment of arbitrary charge distributions inside a spherical cavity. Nowadays, the treatment of solutes of arbitrary shape is possible with the use of powerful computers and numerical methods. In many cases, this is an excellent approximation. The classical electrostatics approach is remarkably successful in reproducing the electrostatic contribution to the solvation free energy of small solutes [26,28] or amino acids [27], as shown by comparisons to free energy simulations with explicit solvent. Applications to biophysical systems are reviewed in Refs. 29 and 30. Because of its importance the next section is devoted completely to this topic.

IV. CLASSICAL CONTINUUM ELECTROSTATICS

A. The Poisson Equation for Macroscopic Media

The continuum electrostatic approximation is based on the assumption that the solvent polarization density of the solvent at a position r in space is linearly related to the total local electric field at that position. The Poisson equation for macroscopic continuum media

Implicit Solvent Models

141

follows from those assumptions about the local and linear electrostatic response of the solvent [31]:

[ε(r) φ(r)] 4πρu(r)

(18)

where φ(r), ρu(r), and ε(r) are the electrostatic potential, the charge density of the solute, and the position-dependent dielectric constant at the point r, respectively. The Poisson equation (18) can be solved numerically by mapping the system onto a discrete grid and using a finite-difference relaxation algorithm [32,33]. Several programs are available for computing the electrostatic potential using this approach, e.g., DelPhi [33,34], UHBD [35], and the PBEQ module [27,36] incorporated in the simulation program CHARMM [37]. Alternatively, one can use an approach based on finite elements distributed at the dielectric boundary (the boundary element method) [38]. Significant improvements can be obtained with this approach by using efficient algorithms for generating the mesh at the dielectric boundaries [39,41]. FAMBE is one program that is available to compute the electrostatic potential using this method [42]. Finally, a different (but physically equivalent) approach to incorporate the influence of a polar solvent, in which the solvent is modeled by a discrete lattice of dipoles that reorient under the influence of applied electric fields, has been proposed and developed by Warshel and coworkers [43].

It is generally assumed that the dielectric constant is uniform everywhere except in the vicinity of the solute/solvent boundary. If all the solute degrees of freedom are treated explicitly and the influence of induced electronic polarization is neglected, the positiondependent dielectric constant ε(r) varies sharply from 1, in the interior of the solute, to εv in the bulk solvent region outside the solute. Such a form for ε(r) follows rigorously from an analysis based on a statistical mechanical integral equation under the assumption that there are only short-range direct correlations in the solvent [44]. To estimate the electrostatic contribution to the solvation free energy, the reaction field Φrf used in Eq. (17) is obtained as the electrostatic potential calculated from Eq. (18) with the nonuniform dielectric constant ε(r), minus the electrostatic potential calculated with a uniform dielectric constant of 1.

Results obtained using macroscopic continuum electrostatics for biomolecular solutes depend sensitively on atomic partial charges assigned to the nuclei and the location of the dielectric boundary between the solute and the solvent. The dielectric boundary can be constructed on the basis of the molecular surface [34,35] or the solvent-accessible surface (constructed as a surface formed by overlapping spheres) [27]. The parametrization of an accurate continuum electrostatic model thus requires the development of optimal sets of atomic radii for the solutes of interest. Various parametrization schemes aimed at reproducing the solvation free energy of a collection of molecules have been suggested [27,28]. From a fundamental point of view, the dielectric boundary is closely related to the nearest density peak in the solute–solvent distribution function [45]. As a consequence, the optimal radius of an atom is not a property of that atom alone but is an effective empirical parameter that depends on its charge, on its neighbors in the solute, and also on the nature of the molecules forming the bulk solvent. In contrast to the radii, the partial charges of the solute are generally taken from one of the standard biomolecular force fields without modification and are not considered as free parameters.

Continuum electrostatic approaches based on the Poisson equation have been used to address a wide variety of problems in biology. One particularly useful application is in the determination of the protonation state of titratable groups in proteins [46]. For