Dielectric Screening Treatment of Electrostatic Solvation - The Journal

11226

J. Phys. Chem. B 1997, 101, 11226-11236

Dielectric Screening Treatment of Electrostatic Solvation Rui Luo,†,‡ John Moult,† and Michael K. Gilson*,†,§ Center for AdVanced Research in Biotechnology, UniVersity of Maryland Biotechnology Institute, 9600 Gudelsky DriVe, RockVille, Maryland 20850, Department of Chemistry and Biochemistry, UniVersity of Maryland at College Park, College Park, Maryland 20742, and National Institute of Standards and Technology, Gaithersburg, Maryland 20899 ReceiVed: July 29, 1997; In Final Form: September 23, 1997X

It has been shown previously that the solvation energies of small molecules can be computed to useful accuracy with continuum models of the solvent that separate the solvation energy into electrostatic and nonpolar parts. The electrostatic part of the solvation energy can be computed with detailed solutions of the Poisson equation, but such calculations are time-consuming. This paper examines the reliability of an approximate method for computing electrostatic solvation energies within the continuum model. The central approximation is that, although the electrostatic field is weakened by the high-dielectric solvent, the shape of the field is not influenced by the solvent. This “dielectric screening” approximation appears to work well when used with a nonlinear, lattice-based model of the solvent. Here, it is developed in the context of linear, continuum electrostatics. The optimal screening function is found to depend upon the distribution of charge within the solute. The present dielectric screening method performs well in comparisons with finite-difference solutions of the Poisson equation and with experimental data, for a range of molecular systems.

1. Introduction Accurately representing the influence of solvent upon molecules in solution is an important and long-standing goal of theoretical chemistry. Ideally, one would like to provide a detailed description of the solvent through explicit simulation of a large number of solvent molecules. This approach is frequently used in molecular dynamics and Monte Carlo simulations of solutions. In many applications, however, the solute is the focus of interest, and the detailed properties of the solvent are not of central importance. In such cases, a simplified representation of the solvent can be employed. This paper concerns a class of simplified solvent models that separates the free energy of solvation into two contributions.1-3 The first contribution is the work of forming in the solvent a nonpolar cavity with the shape of the solute and having van der Waals interactions with the solvent. The second is the work of moving the atomic charges, within a low dielectric body having the shape of the molecule, into the solvent cavity. The first work term is typically estimated from the surface area of the cavity.4 The second work term may be estimated by a classical electrostatic model. The electrostatic model represents the solute as a dielectric body whose shape is defined by atomic coordinates and radii.5-7 The solute contains a set of charges that produce an electrostatic field in and around the solute. The solvent is represented as a dielectric continuum surrounding the solute. The electrostatic fields in such a system may be computed by solving the Poisson equation.8-10 A dissolved electrolyte may be accommodated by use of the PoissonBoltzmann equation instead of the Poisson equation.11 However, the present paper is limited to the case where no dissolved ions are present, i.e., to the Poisson equation. Continuum solvent models of this type provide a well-defined physical picture and may be parametrized to yield solvation free energies * Corresponding author. Phone: (301) 738-6217. FAX: (301) 7386255. E-mail: [email protected]. † Center for Advanced Research in Biotechnology. ‡ University of Maryland at College Park. § National Institute of Standards and Technology. X Abstract published in AdVance ACS Abstracts, November 15, 1997.

S1089-5647(97)02483-8 CCC: $14.00

for small molecules that agree well with experiment.2,3 They also yield relatively accurate values for the pKa’s of difunctional acids and bases,12,13 and of proteins (see, e.g., refs 14-22). In some applications of classical electrostatics to molecular systems, only a relatively small number of molecular conformations must be considered. For such applications, the Poisson equation may be solved very accurately with finite-difference, boundary element, or volume element methods (see, e.g., refs 8, 9, 23-27). In other applications, however, it is necessary to solve the Poisson equation for many different molecular conformations. Use of highly accurate numerical methods then becomes intractable, and computationally faster approximations become valuable. The approximations that have been proposed include image-charge methods,28,29 methods based upon the use of coarse finite-difference grids30 and upon the use of coarse boundary-element representations,31 the generalized Born approximations,32,33 and the induced multipole solvation approximation.34,35 The present paper considers the utility of a different approximation for computing electrostatic solvation energies. This “screening function” approximation has proved successful in the context of the Langevin Dipoles (LD) model,36,37 as now discussed. The LD model provides an alternative method of estimating electrostatic solvation energies. It represents the solvent as a fixed cubic lattice of point polarizabilities outside the solute. These polarizabilities respond to electrical fields according to the Langevin equation.36 Like classical electrostatic models, the method gives good results for the solvation free energies of small organic molecules,38 and has also proved useful for computing the pKa’s of ionizable groups in proteins.36 In the full implementation of the LD model, each point polarizability “feels” both the field of the solute and that of all the other point polarizabilities. An iterative calculation provides a selfconsistent set of point polarizations.36 Obtaining unambiguous results also appears to require evaluation of the lattice polarizations for multiple positions of the lattice relative to the solute.37 For these reasons, the full implementation of the LD model can be somewhat time-consuming. © 1997 American Chemical Society

Dielectric Screening Treatment

J. Phys. Chem. B, Vol. 101, No. 51, 1997 11227

However, a faster noniterative form of the LD model has also been implemented. The noniterative LD model is based upon the observation that the self-consistent lattice polarizations stabilize the system by reducing the energy stored in the electrostatic field, i.e., by weakening the field. Therefore, the field that acts on each point polarizability is weaker than the vacuum electrostatic field due to the solute alone. Although the factor by which the field is weakened actually varies from one lattice site to another, the approximation is made that the fields are weakened by a characteristic factor. It is also assumed that the final field has the same direction as the vacuum field. Thus,

Eint(xk) ≈ E0(xk)/dLD(xk)

Theory. We are interested in the electrostatic work associated with step 2 of the solvation process described in the Introduction. In this step, the charges and the low dielectric body of the solute are transferred from vacuum into the prepared van der Waals cavity in the solvent. The dielectric constants of the solute and solvent are respectively i and s. The work associated with this process results entirely from the change in polarizability around the solute and will therefore be termed the polarization free energy Gpol. From classical electrostatics, this quantity may be written as follows:

Gpol )

(1)

Here Eint(xk) is the field acting at the position xk of lattice site k, where the subscript “int” indicates that this is the “internal” field acting on a polarizability, as opposed to the classical Maxwell field;39 E0(xk) is the vacuum field due to solute alone; and dLD is a position-dependent screening factor that relates the two fields within the noniterative LD model. Use of this approximation completely eliminates the need to compute the many pairwise interactions among the point polarizabilities on the lattice. The only interactions treated explicitly are those between solute charges and the polarizabilities on the lattice. The computations are therefore markedly faster than with the full LD model. Nonetheless, the noniterative approach also yields good agreement with experiment for single ions.37 The utility of the screening concept in the LD approach suggests that it might also be useful in a fast implementation of the classical electrostatics model. In classical electrostatics, however, it is natural to define screening in terms of the Maxwell field, rather than the internal field. Thus, a classical screening factor should presumably be chosen to make the following approximation a good one:

E(x) ≈ E0(x)/d(x)

2. Methods

(2)

Here E(x) and E0(x) are respectively the Maxwell field and the vacuum field as a function of position x; and d(x) is the screening factor for the classical model. The screening factor d(x) might appear to play the role of a dielectric constant. In fact, the screening function will be identical to the dielectric constant for a system in which the actual field is proportional to the Coulomb field, such as a centrosymmetric Born ion. More generally, however, the ratio of the actual field to the vacuum field is not a constant either inside or outside the solute.5 In such cases, the screening factor becomes a parameter that must be optimized to provide realistic solvation energies. It is worth noting that, although dLD and d are defined somewhat differently, they may be connected formally by equating the polarization fields within the two models. The present paper describes the implementation and evaluation of an approximation to the classical electrostatic solvation model that uses the approximation of uniform dielectric screening (DS). The Methods section derives the necessary formulas and details the approximations used. It then presents a numerical implementation of the formulas and describes how the solvent dielectric screening factor is parametrized. The Results section analyzes the numerical performance of the model and assesses its agreement with experiment and with detailed finite difference (FD) solutions of the Poisson equation. The Discussion and Conclusions section summarizes the chief results and their implications and proposes future investigations that follow from this work.

∫

∫

1 1 EF‚DFdV EI‚DIdV 8π 8π

(3)

Here the subscripts I and F indicate the initial and final states of the process; E and D are the Maxwell field and the displacement fields, respectively; and the integrals extend over all of space. Subtraction of the vacuum electrostatic energy, 1/ ∫E ‚D dV, from both terms in eq 3 leaves G 8π 0 0 pol unchanged and permits the use of the following identity:39

-

∫

∫

∫

1 1 1 P‚E0 dV ) E‚D dV E0‚D0 dV 2 8π 8π

(4)

where P is the classical polarization vector. The result is that

Gpol ) -

∫

1 (PF - PI)‚E0 dV 2

(5)

Note that, although E0 is infinite at the point charges used to model the electrostatic field of the solute, the integral in eq 5 is finite. Briefly, this is because the divergence of E0 is nonzero only inside the solute, and the divergence of (PF - PI) is always zero inside the solute. Green’s theorem can be used to show that the integral of the dot-product of these fields does not diverge.40 Equation 5 expresses as a volume integral the electrostatic solvation energy that we wish to compute. Evaluating this integral accurately requires solving the Poisson equation for the field and, thus for the polarization, throughout space. However, eq 2 permits the field to be estimated from the vacuum field, which is readily computed, by use of the screening factor d. Thus,

EI ≈ E0/dI EF ≈ E0/dF

(6)

where the subscripts of d indicate the screening factors appropriate to the initial and final state of the process of electrostatic solvation. The relationship between polarization and field is

P)

( - 1) E 4π

(7)

where is the dielectric constant. Both quantities depend, in general, upon position. Equations 6 and 7 permit eq 5 to be written approximately as

Gpol ≈ -

(

)

1 1 F - 1 I - 1 2 E0 dV 2 4π dF dI

∫

(8)

Although not written explicitly, the values of and d will assume different values for the solute and solvent regions.

11228 J. Phys. Chem. B, Vol. 101, No. 51, 1997

Luo et al.

This expression may be simplified as follows. First, note that F ) I ) i within the solute. Also, it is convenient to assume that the screening factor with the solute is not affected by the transfer from vacuum to solvent. A correction for this approximation is described below. Therefore, the integrand of eq 8 becomes zero within the solute, so the integration may be restricted to the solvent region. Outside the solute, I ) 1 and F ) s. The result is that the second term in the parentheses in eq 8 is zero outside the solute. Equation 8 may therefore be rewritten as follows:

∫solventE20 dV

1 Gpol ≈ - R 2 R≡

(9)

s - 1 4πdF

Born Onsager Gbulk pol ≈ Gpol + Gpol

The quantity R may be viewed as the effective polarizability of the solvent. This approximate expression directly relates the vacuum field in the region occupied by solvent to the polarization free energy which is the quantity we wish to compute. The effective polarizability R will be treated as a fitting parameter that will be adjusted to optimize agreement with accurate numerical solutions of Poisson’s equation for systems of interest. This is equivalent to and somewhat more convenient than adjusting the screening factor. Moreover, R will be adjusted further to correct for the approximation that the screening function does not change within the solute. The polarization energy cannot be evaluated directly with eq 9 in its present form, because the formula requires integrating over an infinite region of space. It is therefore convenient to separate the integral into two parts. The first, Gnear pol , is a shortrange contribution that will be evaluated numerically with a spherical grid around the solute. The second, Gbulk pol , is a longrange contribution from the bulk of the solvent surrounding the spherical grid. The bulk contribution will be evaluated with an analytic approximation. A similar separation of terms is used in LD methods.37 The separation permits eq 9 to be rewritten as bulk Gpol ≈ Gnear pol + Gpol

(10)

∫sphereE20 dV

1 Gnear pol ≡ - R 2

where the integral is the same as that in eq 9, except that it is now restricted to a spherical region around the solute. Computational Methodology. Numerical EValuation of the Polarization Energy. This section describes the evaluation of the polarization energy as approximated in eq 10. The numerical evaluation of Gnear pol with a grid is described and then analytic approximations to Gbulk pol are presented. Integration Grid. The grid is a cubic lattice truncated by a sphere that is centered at the center of geometry of the solute. The grid spacings examined here range from 0.1 to 0.5 Å. The radius of the sphere a is the greatest atom-to-center distance of the solute, plus a “buffer” distance typically between 3 and 9 Å. The choice of grid spacing and buffer distance is discussed below. The integral in eq 10 is evaluated as the following sum over the grid vertexes k:

Gnear pol ≈ RV

∑k ckE20,k

depending upon whether grid point k is in the solute, on the solute-solvent boundary, or in the solvent. Grid points within 1/32 of a grid unit from the van der Waals surface of the solute are considered boundary points. Other grid points are unambiguously inside or outside. The solute-solvent boundary thus is defined by the van der Waals surface of the solute. Evaluation of the Contribution of Bulk Solvent. One of two methods is used to calculate Gpol bulk. The choice depends on the size of the solute. For small molecules, the bulk contribution may be estimated as the sum of the two leading terms of a multipole expansion of the bulk water polarization energy.37 The first term has the form of the Born solvation energy of an ion in a spherical cavity. The second term has the form of the Onsager solvation energy of a dipole in a spherical cavity. Thus

(11)

Here V is the volume of one grid box, E0,k is the vacuum field at grid point k, and ck assumes values of 0, 0.5, or 1,

(12)

where

GBorn pol ≡ -

( )

(13)

µ2 s - 1 a3(2s + 1)

(14)

Q2 1 12a s

GOnsager ≡pol

Here Q is the net charge of the solute, a is the radius of the spherical grid, and µ is the total dipole moment with respect to the center of the sphere. The approximation in eqs 12-14 will be termed the Born/Onsager method. It is accurate only when the quadrupole of the solute is small, or when the solvent buffer around the solute is thick. This makes it unwieldy for large molecules, such as proteins. Therefore, Friedman’s image charge method28 is used for large molecules. This method is accurate even when the solute occupies most of the spherical grid.28 The performance of the Born/Onsager and Friedman methods is examined below. EffectiVe Polarizability R. This section has so far described numerical procedures for computing Gpol, but it has not yet specified the value or values of R to be used in the integral expression for Gnear pol . Obtaining good agreement between the present DS approximation and full solutions of the Poisson equation for real molecules will require using appropriate values of this parameter. It can be shown that the optimal value of R depends upon the charge distribution of the solute. This can be appreciated by considering two simple solutes. First, the field that a Born ion generates in the solvent is

E ) E0/s

(15)

Comparison with eq 2 shows that the optimal screening factor in this case is simply s. From eq 9, this corresponds to the following effective polarizability:

Rion ≡

s - 1 4πs

(16)

When this expression is used with eq 9 to compute the solvation energy of a Born ion, the agreement with the analytic Born energy is limited only by numerical precision. This is because the screening function involves no approximation, and because an integral in the solvent region captures the entire electrostatic energy. That is, the value of the integrand P ‚ E0 within the solute is not affected by solvating the ion. A different screening factor is appropriate for the case of an ideal dipole in a spherical cavity. It is readily shown that the field such a solute generates in the solvent is39


E)


3E0 (2s + i)

(17)

By inspection, the screening factor in this case is (2s + i)/3. This expression can be combined with eq 9 to write a preliminary expression for the effective polarizability of the solvent:

R′dipole ≡

3 s - 1 4π 2s - i

(18)

With this formula, the integral of -1/2RE20 dV in the solvent region would agree with the analytic result for this region to within numerical precision. However, the quantity P ‚ E0 within a dipolar solute is influenced by solvation, in contrast with the case of the Born ion. Therefore, Gpol computed with R′dipole as an integral in the solvent region would not agree perfectly with the analytic formula for the solvation energy of a dipolar solute. An expression for R that corrects for the internal integral can be derived by referring to the analytic formula for transfer of an ideal dipole of magnitude µ in a spherical dielectric body of radius a and internal dielectric constant i from vacuum to a solvent with dielectric constant s:

3(s - 1) µ2 Gdipole ) 3 (2 + i)(2s + i) a

(19)

This expression can be derived as a generalization of the case of a dipolar solute with i ) 1.39 Following eq 9, it is now required that

1 Gpol ) - Rdipole 2 )-

∫a∞E20 dV ) Gdipole

(20)

4πµ2 Rdipole 3a3

so

Rdipole )

s - 1 3 3 4π 2 + i 2s - i

(21)

This is identical with the preliminary expression for R′dipole, except for an additional factor of 3/(2 + i). Given that the external integral is exact when R′dipole is used, it is evident that this new factor corrects for the change in the internal polarization of the solute when it is transferred from vacuum to solvent. Accordingly, it disappears when i ) 1, i.e., when there is no change in the internal polarization upon solvation. For the physically reasonable value i ) 2, the correction is 0.75. The present study uses the corrected value Rdipole. Comparison of the screening factors d given above for an ion and a dipole shows that a dipolar field is screened more than a ionic field, as noted previously.5 This result may be explained intuitively in terms of boundary charges. Charges inside a solute induce boundary charges at the solute-solvent surface. The boundary charges are of opposite sign to the inducing charges. Therefore, they weaken or “screen” the field of the source charges. A dipolar solute induces both positive and negative boundary charges, and the partial cancellation of these charges leads to reduced screening of the field. Because the value of R appropriate for a Born ion differs significantly from that for an Onsager dipole, it is reasonable to suppose that the optimal screening factor for a real molecule will depend upon its charge distribution. It seems likely that

the screening factor of F for a solute with a single charge is the maximal possible value, because the boundary charges induced by the charge are unopposed. It also seems likely that the screening factor for an ideal dipole is the minimal value. This is because an ideal dipole consists of two infinitely close charges of opposite sign, an arrangement that presumably results in the maximal cancellation of induced boundary charges. These conjectures were tested empirically by determining the screening factors that yielded optimal agreement between Gpol computed for small molecules with full numerical solutions of the Poisson equation and by the present DS method. As anticipated, the optimal screening factors consistently fall between the optimal values for a single ion and a dipole. These results motivated the development of an empirical expression for the effective polarizability R, whose value ranges between Rion and Rdipole, and depends upon the asymmetry of the charge distribution of the solute. An index of charge asymmetry γ is defined as follows:

{

|Q+/Q-|, γ ≡ |Q-/Q+|, ∞

|Q+| g |Q-| |Q+| < |Q-| Q+ ) 0 or Q- ) 0

(22)

Here, |Q+| and |Q-| are respectively the absolute values of the total positive and negative partial atomic charges of the solute atoms. Then the effective polarizability is given by

(

)(

)

1 1 1 1 1 ) + R Rion Rdipole Rion k (γ - 1) + k (γ - 1)1/2 + 1 1 2 (23) where k1 and k2 are constants whose values must be adjusted to optimize the agreement of the present DS model with the results of detailed solutions of the Poisson equation. This adjustment is described below. Equation 23 yields R ) Rdipole when the net charge is zero, i.e., when γ ) 1. It yields R ) Rion when only positive or only negative charge is present, i.e., when γ f ∞. Computational Parameters. Except as otherwise specified, all calculations use a solute dielectric constant of 1.0 and a solvent dielectric constant of 78.5, the thickness of the spherical buffer region surrounding the solute is set to 6 Å, and the Born/ Onsager method is used to calculate Gpol bulk. The solutesolvent boundary is defined by the van der Waals surface of the solute. FD solutions of the Poisson equation are obtained with the program UHBD.41 For small molecules, including the dipeptides (see below), the FD grids are dimensioned to 50 × 50 × 50, and the grid spacings are roughly 0.25 Å or less. The FD calculations for the protein fragments use grids of 75 × 75 × 75, with grid spacings of roughly 0.35 Å. For the proteins, the FD grids are 150 × 150 × 150 with grid spacings of roughly 0.35 Å. All FD solutions use dielectric boundary smoothing to improve accuracy.42 The boundary potentials of the grid are set with Coulomb’s law, with the dielectric constant of the solvent. The ionic strength is zero. Further details of specific calculations are included in Results. 3. Results This section examines the sensitivity of the DS method to the spacing of the grid used for the numerical integration, and to the thickness of the “buffer” region that defines the size of the spherical grid around the solute (see above). It is demonstrated that the method precisely reproduces analytic results for


Luo et al.

TABLE 1: Sensitivity of Gpol to Thickness of Solvent Buffera Born/Onsager

image

molecule

3Å

6Å

9Å

3Å

6Å

9Å

ser thr lys cys met asp glu asn gln arg phe tyr trp his aspglutyrlys+ arg+ his+

-40.52(-0.67) -39.26(-0.49) -37.11(-0.29) -18.48(-0.40) -20.82(-0.72) -50.83(-0.35) -50.51(-0.27) -69.81(-1.69) -68.90(-1.28) -74.94(-0.24) -16.03(0.00) -47.40(-0.22) -45.31(-0.06) -70.85(-1.39) -349.68(-107.99) -345.82(-100.67) -343.53(-91.13) -313.59(-81.88) -319.47(-74.07) -297.01(-89.69)

-40.53(-0.20) -39.30(-0.16) -37.14(-0.11) -18.49(-0.12) -20.87(-0.25) -50.94(-0.12) -50.61(-0.10) -69.87(-0.57) -68.96(-0.46) -74.99(-0.11) -16.08(0.00) -47.52(-0.08) -45.43(-0.03) -70.97(-0.50) -350.55(-73.12) -346.26(-69.47) -345.64(-64.51) -314.80(-59.99) -322.70(-55.92) -300.07(-64.36)

-40.54(-0.09) -39.30(-0.07) -37.15(-0.06) -18.49(-0.05) -20.87(-0.12) -50.95(-0.05) -50.62(-0.05) -69.88(-0.25) -68.96(-0.22) -75.00(-0.05) -16.09(0.00) -47.54(-0.04) -45.45(-0.01) -70.99(-0.23) -351.05(-55.32) -346.63(-53.13) -346.90(-50.07) -315.58(-47.41) -324.63(-44.92) -301.71(-50.18)

-40.53(-0.69) -39.31(-0.54) -37.15(-0.33) -18.49(-0.41) -20.88(-0.79) -50.97(-0.50) -50.64(-0.40) -69.88(-1.76) -68.97(-1.35) -75.00(-0.31) -16.10(-0.07) -47.56(-0.38) -45.48(-0.24) -71.00(-1.53) -348.42(-106.73) -344.65(-99.50) -342.92(-90.53) -312.64(-80.92) -318.66(-73.26) -296.19(-88.88)

-40.53(-0.20) -39.30(-0.17) -37.15(-0.12) -18.49(-0.13) -20.88(-0.26) -50.96(-0.14) -50.63(-0.12) -69.88(-0.57) -68.97(-0.48) -75.00(-0.12) -16.09(-0.02) -47.55(-0.11) -45.47(-0.06) -70.99(-0.53) -349.65(-72.21) -345.41(-68.61) -344.93(-63.80) -314.06(-59.26) -322.03(-55.25) -299.32(-63.61)

-40.54(-0.09) -39.30(-0.08) -37.15(-0.06) -18.49(-0.05) -20.88(-0.12) -50.96(-0.06) -50.63(-0.05) -69.88(-0.25) -68.97(-0.22) -75.00(-0.06) -16.09(0.00) -47.55(-0.05) -45.46(-0.03) -70.99(-0.24) -350.36(-54.63) -345.96(-52.47) -346.31(-49.47) -314.99(-46.82) -324.08(-44.36) -301.10(-49.57)

near bulk bulk a Sensitivity of G pol ) Gpol + Gpol to thickness of solvent buffer for two different approximations. Values in parentheses are Gpol . Energies are in kJ/mol. (Same in all the following tables.) Compounds are analogs of the specified amino acid side chains, and begin with the Cβ atom.

a Born ion and a spherical dipolar solute. The method is then shown to yield electrostatic solvation energies of small molecules, protein fragments, and proteins that agree well with those from FD solutions of the Poisson equation, as well as with experimental data. Good agreement also is obtained for the variation of the electrostatic solvation energy of alanine dipeptide with conformation. Finally, the method is applied to the computation of pKa shifts for small dicarboxylic acids. Numerical Precision of the Method. SensitiVity to Grid Spacing. The present method involves a numerical evaluation of the approximate formula for Gpol in eq 10. Obtaining a precise value requires that the grid used in the numerical integration not be too coarse. The necessary grid spacing is arrived at here by examining the variation of computed values of Gpol as a function of the grid spacing for a Born ion of 1e charge, and a spherical dipolar solute of dipole moment 1 e Å. The cavity radius is 2 Å in both cases, and a 6 Å buffer thickness is used initially; therefore, the outer radius of the integration grid is 8 Å. The value of Gpol is computed with grid spacings ranging from 0.05 to 1.0 Å. The results for the ion and the dipole are plotted in Figure 1a and 1b, respectively. The plots suggest that the calculations are converged for grid spacings e0.1 Å and that the error remains small for grid spacings of 0.2-0.3 Å. These results do not change significantly when the origin of the grid is shifted relative to the solute, or when the thickness of the buffer region is changed to 3 or 9 Å. Grid spacings of 0.2 Å or less are used for all calculations presented in the remainder of the paper, except as otherwise noted. SensitiVity to Thickness of the Buffer. The numerical results of the present method can be influenced not only by the grid spacing, but also by the thickness of the buffer region around the solute. However, for the simple ion and dipole cases examined above, the thickness of the buffer region is not important because the bulk solvent contribution, Gbulk pol , is given accurately and precisely by the Born/Onsager treatment. (See the following section.) Tests of the influence of the buffer thickness must use more complex solutes in order to be instructive. Therefore, the present analysis considers analogues of both neutral and ionized amino acid side chains (Table 1),

for which both the Born/Onsager method and the Friedman image method are approximate. The PARSE parameter set2 is used for the partial charges and radii. The effective polarizability of the solvent R is computed according to eq 23, with optimized values of k1 and k2 obtained from the parametrization described below. Table 1 presents the values of Gpol for buffer thicknesses of 3, 6, and 9 Å, computed with both the Born/ Onsager approximation and the Friedman image approximation. As shown, the results depend only weakly upon the thickness of the buffer. Upon the basis of this analysis, the calculations of Gbulk pol in the remainder of this paper use a 6 Å buffer and the Born/Onsager approximation, except as otherwise noted. Accuracy for Ideal Spherical Solutes. It is expected that the present method will yield highly accurate values of Gpol for a Born ion and a spherical dipolar solvent. This is because the effective polarizabilities for these cases, Rion and Rdipole, are without approximation, as is the Born/Onsager treatment of Gbulk pol . It is confirmed that this is indeed the case. Comparisons between analytic values of Gpol with the values from the present numerical method, for various charges and radii of an ideal Born ion and similar comparisons for a spherical, dipolar solute indicate that the numerical and analytic results agree to high precision, as expected: the largest absolute error in any of these comparisons is only 0.01 kJ/mol, and the largest fractional error is only 0.01%. This result supports the validity of the theory underlying the present method. Application to Molecular Systems.This section of Results examines the performance of the DS method for molecules. The first subsection describes the parametrization of eq 23 for the effective polarizability. Subsequent sections compare the DS method with accurate FD solutions of the Poisson equation and, where possible, with experiment. Parametrization of the EffectiVe Polarizability. The final step in formulating the DS method is to establish values for the two adjustable parameters k1 and k2 in eq 23. These are adjusted to optimize the agreement with FD solutions of the Poisson equation for the seven small, organic ions listed in Table 2. The resulting values of k1 and k2 are 0.5833 and 1.1666, respectively. With these values, the correlation coefficient for



TABLE 2: Parametrization of Effective Polarizability ra molecule

γ

GDS pol

GFD pol

∆Gpol

relative error

asp glutyrcyslys+ arg+ his+

11.000 11.000 2.600 ∞ 4.125 1.571 2.427

-350.53 -346.27 -345.64 -326.21 -314.80 -322.70 -300.08

-361.95 -362.07 -355.76 -338.58 -314.80 -319.85 -289.97

11.41 15.80 10.12 12.37 0.00 2.84 10.12

3.15% 4.36% 2.84% 3.65% 0.00% 0.89% 3.49%

-

a Parametrization of effective polarizability R. γ: charge asymmetry, defined by eq 22. GDS pol: electrostatic solvation free energy from DS method. GFD pol: electrostatic solvation free energy from FD method. FD FD ∆Gpol: |GDS pol - Gpol|. relative error: ∆Gpol/Gpol. See Table 1 for further details.

Figure 2. (a) Correlation and linear regression of DS electrostatic solvation energies with those of FD method for neutral molecules. (b) Correlation and linear regression of the DS solvation free energies with experimental results for the same molecules, but excluding ethanediol and ethanediamine (see text).

Figure 1. (a) Sensitivity of numerical integration to grid spacing for a Born ion. Solid line: DS results. Dashed line: analytical result. (b) Sensitivity of numerical integration to grid spacing for an ideal dipolar solute. Solid line: DS results. Dashed line: analytical result.

the DS and FD results is 0.99. The average unsigned error for the seven molecules is 8.95 kJ/mol, and the relative errors are consistently less than 4.5%. Table 2 provides further details of the comparison. Subsequent calculations in this paper use the optimized values of k1 and k2. SolVation Energies of Small Molecules. The DS method is used to compute the electrostatic solvation energies Gpol of the 116 small molecules used to train and test the validity of the PARSE parameter set.2 This parameter set was optimized to reproduce experimental solvation energies when used with the Poisson-Boltzmann electrostatics model. The present calculations use the PARSE radii and charges and also use the same molecular conformations used in developing and testing the parameter set.2 As dictated by the PARSE model, the internal dielectric constant of the solutes is set to 2. Note that this

parameter influences the results of the DS method through its effect upon Rdipole (eq 23), and thus upon R, the effective polarizability. To allow comparison with experimental solvation energies, the electrostatic solvation energies obtained here are added to the nonpolar cavity terms given previously for these molecules.43 The DS method performs well for both neutral and ionized molecules, as now shown. The results for neutral are detailed in Table 3 and are summarized by the scatter plots of Figure 2. The agreement between DS and FD results is examined in Figure 2a, which presents a scatter plot of the electrostatic energies computed with the present method versus the FD results. A linear regression fit of these data yields a correlation coefficient of 0.99, a slope of 0.98, and y-intercept of -1.73 kJ/mol. The mean unsigned deviation between the DS and FD results is 1.45 kJ/mol. The agreement of the DS results with experiment is examined in Figure 2b, which presents a scatter plot of the DS energies incremented by the appropriate cavity terms, versus the measured solvation energies. A linear regression fit of these data yields a correlation coefficient of 0.97, a slope of 0.89, and y-intercept of -2.41 kJ/mol. The mean unsigned error for these neutral molecules is 1.87 kJ/mol. The results for ionized molecules are detailed in Table 4, and are summarized in Figure 3. The agreement between DS and FD results is examined in Figure 3a, which presents a scatter plot of the electrostatic energies computed with the present method versus the FD results. A linear regression fit of these data yields a correlation coefficient of 0.99, a slope of 1.00, and y-intercept of -6.15 kJ/mol. The mean unsigned deviation


Luo et al.

TABLE 3: Calculated and Measured Solvation Free Energies for Neutral Moleculesa molecule

GDS pol

GFD pol

Gcav

DS Gsolv

exp Gsolv

∆Gerror,DS

molecule

methanol ethanol propanol butanol isopropyl alcohol 2-butanol 3-methyl-1-butanol ammonia methylamine ethylamine propylamine N-butylamine dimethylamine diethylamine methylthiol ethylthiol methylethyl sulfide dimethyl sulfide diethyl sulfide acetone 2-butanone 2-pentanone 3-pentanone 3-methyl-2-butanone 2,4-dimethyl-3-pentanone acetaldehyde proprionaldehyde acetic acid proprionic acid butyric acid acetamide proprionamide N-methylacetamide N-propylguanidine benzene toluene ethylbenzene pyridine 4-methylpyridine 2-methylpyridine phenol p-cresole 2-methylphenol methylindole methylimidazole pentanol 4-methyl-2-pentanol 2-pentanol 2-methyl-2-butanol 2,3-dimethylisobutanol 3-pentanol hexanol 3-hexanol 2-methyl-3-pentanol heptanol 4-heptanol 2-methyl-2-propanol

-29.89 -28.72 -29.01 -28.76 -26.46 -26.54 -28.67 -23.62 -28.21 -27.30 -27.17 -27.17 -27.55 -24.58 -13.63 -13.38 -15.38 -16.47 -14.21 -25.79 -25.21 -24.70 -24.70 -24.58 -23.41 -23.78 -22.95 -37.03 -37.41 -37.03 -51.71 -51.16 -53.67 -55.43 -12.04 -11.83 -11.62 -28.72 -29.47 -26.38 -35.32 -35.20 -34.03 -33.40 -52.50 -28.97 -26.67 -26.63 -25.37 -25.25 -26.50 -28.84 -26.29 -26.25 -29.34 -26.08 -26.00

-31.18 -30.26 -30.18 -30.26 -28.21 -28.88 -30.30 -24.49 -29.05 -28.42 -28.51 -28.47 -28.97 -26.88 -14.21 -14.04 -15.93 -16.93 -15.01 -26.13 -25.46 -25.29 -24.79 -24.79 -23.70 -24.16 -23.41 -37.29 -36.87 -36.99 -51.12 -50.75 -54.88 -59.44 -13.96 -13.84 -13.71 -30.81 -31.68 -28.55 -38.04 -37.95 -37.16 -36.57 -53.75 -30.05 -28.55 -28.59 -27.76 -27.55 -28.63 -30.22 -28.63 -28.13 -30.10 -28.34 -28.09

7.40 8.15 8.82 9.49 8.74 9.40 9.99 6.52 7.57 8.23 8.90 9.61 8.40 9.82 7.94 8.61 9.36 8.65 10.03 8.69 9.36 9.95 9.99 9.86 10.83 7.90 8.61 8.23 8.95 9.53 8.32 8.99 9.20 10.28 8.99 9.74 10.41 8.82 9.66 9.66 9.28 10.07 9.99 10.66 9.15 10.16 10.53 10.07 9.91 10.37 10.07 10.83 10.74 10.58 11.49 11.41 9.28

-22.49 -20.57 -20.19 -19.27 -17.72 -17.14 -18.68 -17.10 -20.65 -19.06 -18.27 -17.56 -19.14 -14.76 -5.68 -4.77 -6.02 -7.82 -4.18 -17.10 -15.84 -14.76 -14.71 -14.71 -12.58 -15.88 -14.34 -28.80 -28.47 -27.50 -43.39 -42.18 -44.48 -45.14 -3.05 -2.09 -1.21 -19.90 -19.81 -16.72 -26.04 -25.12 -24.03 -22.74 -43.35 -18.81 -16.13 -16.55 -15.47 -14.88 -16.43 -18.02 -15.55 -15.67 -17.85 -14.67 -16.72

-21.23 -20.48 -20.19 -19.73 -19.90 -19.14 -18.48 -18.02 -19.10 -18.81 -18.35 -18.31 -17.93 -17.01 -5.18 -5.43 -6.23 -6.44 -5.98 -16.09 -15.22 -14.76 -14.25 -13.54 -11.45 -14.63 -14.38 -28.01 -27.04 -26.58 -40.63 -39.38 -42.13 -45.65 -3.64 -3.18 -3.34 -19.65 -20.61 -19.35 -27.67 -25.62 -24.54 -24.70 -42.84 -18.68 -15.59 -18.35 -18.52 -16.39 -18.18 -18.22 -17.05 -16.26 -17.72 -16.76 -18.85

1.25 0.08 0.00 0.46 2.17 2.01 0.21 0.92 1.55 0.25 0.08 0.75 1.21 2.26 0.50 0.67 0.21 1.38 1.80 1.00 0.63 0.00 0.46 1.17 1.13 1.25 0.04 0.79 1.42 0.92 2.76 2.80 2.34 0.50 0.59 1.09 2.13 0.25 0.79 2.63 1.63 0.50 0.50 1.96 0.50 0.13 0.54 1.80 3.05 1.50 1.76 0.21 1.50 0.59 0.13 2.09 2.13

2-methyl-2-pentanol octanol pentaneamine hexaneamine N,N-dipropylamine N,N-dibutylamine 4-methyl-2-pentanone 4-heptanone 5-nonanone 2-hexanone 2-heptanone 2-octanone 2-nonanone butanal pentanal hexanal heptanal octanal nonanal N-methylformamide cyclopentanol cyclohexanol cycloheptanol pyrrolidine piperazine pipyridine propylbenzene o-xylene m-xylene p-xylene naphthalene anthracene phenanthracene 1,2,4-trimethylbenzene methylethylbenzene butylbenzene methylpropylbenzene dimethyethylbenzene dimethylpropylbenzene dimethylethylphenol 3-methylpyridine 2-ethylpyridine 3-ethylpyridine 4-ethylpyridine 2,4-dimethylpyridine 3,4-dimethylpyridine 3,5-dimethylpyridine 2,3-dimethylpyridine 2,5-dimethylpyridine 2,6-dimethylpyridine dimethylethylpyridine 2-methylpyrazine 1,2-ethanediamine 1,2-ethanediol av dev av error

GDS pol

GFD pol

-25.54 -29.18 -27.84 -27.63 -24.66 -25.21 -24.91 -24.54 -24.45 -25.21 -25.41 -25.25 -25.29 -23.53 -23.49 -23.45 -23.53 -23.41 -23.45 -34.74 -27.88 -27.88 -27.17 -26.88 -49.11 -26.38 -11.66 -11.87 -11.66 -11.75 -17.31 -22.07 -22.15 -11.58 -11.41 -11.62 -11.41 -11.33 -11.29 -41.97 -28.59 -26.21 -28.17 -28.63 -25.92 -28.01 -27.88 -26.08 -25.92 -24.03 -27.88 -41.47 -55.47 -52.75

-27.80 -30.10 -28.47 -28.47 -27.13 -27.46 -25.12 -24.62 -24.58 -25.25 -25.16 -25.21 -25.29 -23.37 -23.37 -23.37 -23.32 -23.45 -23.58 -35.20 -29.09 -29.05 -28.93 -28.55 -51.41 -28.09 -13.63 -13.88 -13.79 -13.71 -19.69 -25.16 -25.21 -13.79 -13.54 -13.63 -13.54 -13.42 -13.38 -37.37 -30.05 -28.21 -29.76 -30.22 -28.13 -29.68 -29.34 -27.92 -27.84 -26.42 -29.76 -43.47 -56.81 -53.96 1.45

Gcav

DS Gsolv

exp Gsolv

∆Gerror,DS

10.53 12.16 10.28 10.95 11.16 12.50 10.49 11.24 12.58 10.70 11.37 12.04 12.71 9.28 9.95 10.62 11.29 11.95 12.62 8.40 9.40 9.91 10.37 9.03 9.40 9.53 11.08 10.37 10.53 10.53 10.45 11.87 11.79 11.12 10.95 11.75 11.54 11.33 11.87 11.66 9.61 10.28 10.28 10.24 10.41 10.20 10.37 10.03 10.41 10.41 11.20 9.49 8.69 8.32

-15.01 -17.01 -17.56 -16.68 -13.50 -12.71 -14.42 -13.29 -11.87 -14.50 -14.04 -13.21 -12.58 -14.25 -13.54 -12.83 -12.25 -11.45 -10.83 -26.33 -18.48 -17.97 -16.80 -17.85 -39.71 -16.85 -0.59 -1.50 -1.13 -1.21 -6.86 -10.20 -10.37 -0.46 -0.46 0.13 0.13 0.00 0.59 -30.30 -18.98 -15.93 -17.89 -18.39 -15.51 -17.81 -17.51 -16.05 -15.51 -13.63 -16.68 -31.98 -46.77 -44.43

-16.43 -17.10 -17.14 -16.85 -15.30 -13.92 -12.79 -12.25 -11.16 -13.75 -12.71 -12.04 -10.41 -13.29 -12.67 -11.75 -11.16 -9.57 -8.69 -41.80 -22.95 -22.86 -22.95 -22.91 -30.85 -21.36 -2.22 -3.76 -3.47 -3.39 -9.99 -17.68 -16.51 -3.59 -1.25 -1.67 -1.88 -1.84 -0.75 -24.75 -19.94 -18.10 -19.23 -18.27 -20.31 -21.82 -20.23 -20.19 -19.73 -19.23 -18.64 -23.07 -31.77 -32.02

1.42 0.08 0.42 0.17 1.80 1.21 1.63 1.04 0.71 0.75 1.34 1.17 2.17 0.96 0.88 1.09 1.09 1.88 2.13 15.47 4.47 4.89 6.14 5.06 8.86 4.51 1.63 2.26 2.34 2.17 3.13 7.48 6.14 3.13 0.79 1.80 2.01 1.84 1.34 5.56 0.96 2.17 1.34 0.13 4.81 4.01 2.72 4.14 4.22 5.60 1.96 8.90 15.01 12.41 1.87

a GDS: electrostatic solvation free energy from DS method. GFD: electrostatic solvation free energy from FD method. G : computed cavity cav pol pol DS exp exp component of solvation free energy (ref 43). GDS solv: Gpol + Gcav. Gsolv: Measured solvation free energy (refs 43, 55, 56, and 57). ∆Gerror,DS: Gsolv DS DS FD Gsolv. av dev.: ∑|Gpol - Gpol|/N. av err.: ∑|∆Gerror,DS|/N. See Table 1 for other details.

between the DS and FD results is 5.74 kJ/mol. The agreement of DS with experiment is examined in Figure 3b which presents a scatter plot of the DS energies incremented by the appropriate cavity terms, versus the measured solvation energies. A linear regression fit of these data yields a correlation coefficient of 1.00, a slope of 1.02, and y-intercept of 5.07 kJ/mol. The mean unsigned error for these ionized molecules is 2.22 kJ/mol. This is greater than the error for the neutral molecules, but the fractional error is smaller. As summarized in the previous two paragraphs, the electrostatic solvation energies computed with the DS method agree well with FD results. Accordingly, the solvation energies obtained by adding cavity-formation contributions to the

TABLE 4: Calculated and Measured Solvation Free Energies for Ionic Moleculesa molecule ammonium cysaspglutyrlys+ his+ arg+ av dev av error a

GDS pol -348.11 -325.96 -341.17 -336.99 -319.48 -297.70 -277.47 -286.87

GFD pol

Gcav

GDS solv

Gexp solv

∆Gerror,DS

-347.57 7.15 -340.96 -340.80 0.17 -330.85 7.86 -318.10 -320.98 2.88 -348.82 8.11 -333.06 -337.12 4.05 -345.60 8.82 -328.17 -330.72 2.55 -328.80 10.03 -309.45 -313.54 4.10 -302.26 9.86 -287.83 -289.42 1.59 -280.44 9.24 -268.23 -268.06 0.17 -294.27 10.28 -276.59 nab na 5.73 2.22

See previous tables for details. b na ) not applicable.



Figure 3. Same as Figure 2, but for ionized molecules.

electrostatic solvation energies agree well with experiment. However, the computations show particularly large errors relative to experiment for the last two molecules listed in Table 3, 1,2-ethanediol and 1,2-ethanediamine. This is despite the good agreement between the FD and the DS results for these molecules. It has been proposed previously that the discrepancy between the PARSE results for these molecules and the measured solvation energies results from a change in conformation upon transfer from vacuum to solution.43 Corrections for the influence of conformational change upon the solvation energies of ethanediol and ethanediamine have been estimated at 16.72 and 9.20 kJ/mol, respectively.43 When corrected by these amounts, the solvation energies computed with the DS method agree better with experiment: the errors of ethanediol and ethanediamine fall to 1.71 and 5.73 kJ/mol, respectively. It is of interest to place the performance of the DS method into the context of previously published methods for computing

the solvation energies of small molecules. This is done in Table 5, which shows that the DS method with PARSE parameters compares favorably with six other continuum or grid-based solvation models. The comparison must be interpreted cautiously, because the various models use several different sets of atomic parameters. However, it would appear that the DS method is a competitive method from the standpoint of accuracy. SolVation Energies of Protein Fragments and Proteins. Experimentally measured solvation energies are not available for proteins and protein fragments. Rather than comparing with experiment, therefore, this section examines the agreement of the DS model with FD solutions of the Poisson equation for these entities. The protein fragments examined are several short peptides extracted from dihydrofolate reductase and from Streptomyces griseus proteinase A. The atomic coordinates are drawn from the Protein Data Bank44 files 3DFR45 and 2SGA.46 The peptides are excised from loop regions of the proteins by cutting the terminal peptide bonds, and no change is made to the resulting termini. Polar hydrogens are added with the program BHYDR2 based on the algorithm in ref 47. Atomic charges and radii are those of the PARSE parameter set.2 As shown in Table 6, the agreement with FD calculations is reasonably good, especially for the charged fragments. Even for the neutral fragments, however, the relative errors are within 10%. The results of a similar comparison for several small proteins are reported in Table 7. Here, polar hydrogens are added with the program Quanta 4.11,48 with partial charges from the CHARMm 22 parameter sets.49 The atomic radius of each atom is set to the Rmin Lennard-Jones parameter of the appropriate atom type, with the exception of the hydrogen atoms, whose radii are set to 1.2 Å. To save computer time, the spacing of the integration grid is increased slightly to 0.3 Å, and the diameter of the grid is set to the largest solute dimension plus a 3 Å buffer. The polarization energy of the bulk solvent Gbulk pol is computed with the Friedman image-charge method. As shown in Table 7, the DS results are within about 5% of the FD results, and the DS method requires less computer time. Electrostatic SolVation Energy as a Function of Conformation. The studies presented above demonstrate that the DS method works well for a variety of different molecules. This section examines the ability of the DS method to yield the solvation energies of two flexible molecules, alanine dipeptide and aspartate dipeptide, as a function of conformation. These molecules differ only in the side-chain attached to the central R-carbon atom. The parameters and starting coordinates of both molecules were created with CHEMNOTE in QUANTA 4.11,48 with the CHARMm 22 parameter set.49 Electrostatic solvation energies are computed as the two main chain torsion angles vary between 0 and 360° in 10° steps. The resulting

TABLE 5: Solvation Free Energies Computed with Various Methodsa molecule

DS

GB

IMPS

FDPB

PDLD

QM/GB

QM/FD

expt

methanol ethanol 2-propanol acetone acetic acid acetamide benzene toluene phenol av error

-21.32 -20.48 -20.06 -17.14 -28.84 -43.47 -2.93 -2.09 -25.92 0.92

-25.92 -21.74 -17.97 -13.38 -27.17 -44.31 -4.18 -0.42 -26.33 2.22

-19.23 -7.94 -7.11 -12.54 na -39.71 na na na 6.44

-22.57 -20.90 -19.23 -16.30 -27.59 -40.96 -4.60 -3.34 -26.75 0.54

-24.66 -21.32 nab na -28.84 -46.40 -2.93 -5.02 -25.50 2.22

-22.15 -19.23 -15.05 -19.65 -36.78 na -1.25 -0.42 -22.57 3.72

-21.32 -18.81 -18.81 -18.39 -32.19 -44.31 na -4.60 -26.33 1.92

-21.32 -20.48 -20.06 -16.30 -28.01 -40.55 -3.76 -3.34 -27.59

a DS: Current method with PARSE parameters. GB: Generalized Born method (ref 58). IMPS: Inducible multiple solvation model (refs 59, and 60) with OPLS parameters. FDPB: Finite difference Poisson-Boltzmann with PARSE parameters (ref 43). PDLD: Protein dipoles Langevin dipoles model (ref 38). QM/GB: Quantum mechanical solute with generalized Born solvation (ref 61). QM/FD: Quantum mechanical solute with FDPB solvation model (ref 62). Exptl: experimental data (refs 55, 56, and 57). av. error: average unsigned error. b na ) not applicable.


Luo et al.

TABLE 6: Electrostatic Solvation Energies for Protein Loopsa crystal structure

position

sequence

charge status

GDS pol

3DFR 3DFR 3DFR 3DFR 3DFR 3DFR 2SGA 2SGA

20-23 27-30 51-54 64-67 120-123 136-139 113-116 132-136

PWHL LHYF ARPL HQED GSFE FTKV GRVY VGQAV

neutral neutral charged charged charged charged charged charged

-163.35 -196.88 -414.66 -982.51 -475.39 -541.10 -478.23 -200.43

GFD pol

∆Gpol

-164.48 1.13 -219.41 22.57 -410.68 3.93 -937.66 44.89 -465.74 9.66 -525.72 15.42 -475.64 2.59 -196.21 4.22

a Position: sequence numbering of residues in original protein. Sequence: Amino acid sequence of peptide. Charge status: “neutral” unless ionized group(s) present. See previous tables for other details.

TABLE 7: Electrostatic Solvation Energies and CPU Usage of Small Proteinsa crystal structure

no. of atoms

GDS pol (CPU)

GFD pol (CPU)

1BPI 1PGA 1CRN

577 535 396

-6065 (563.69) -4493 (416.58) -1593 (237.92)

-6094 (1308.52) -4598 (1306.45) -1668 (1266.17)

a Electrostatic solvation energies of small proteins of known crystal structure. CPU: approximate computer time (in seconds) for each calculation is measured on a single processor of an SGI Onyx InfiniteReality with 4 194-MHz R10000 processors. See text and other tables for details.

Figure 5. Variation of electrostatic solvation energy with torsion angles of aspartate dipeptide for FD method (a) and DS method (b).

Figure 4. Variation of electrostatic solvation energy with torsion angles of alanine dipeptide for FD method (a) and DS method (b).

conformations are used as is, i.e., they are not reoptimized for each set of dihedral angles. Figures 4 and 5 compare the DS and FD results as functions of the torsion angles. The shape of the energy-surfaces generated by the two methods are very similar. Figures 6a,b compare the energy surfaces further through scatter plots of the DS versus the FD results. The correlation coefficients are 0.99 for both molecules.

Figure 6. Correlation of DS and FD results for alanine dipeptide (a) and aspartate dipeptide (b) as a function of conformation (see Figure 5).



TABLE 8: pKa’s of Dicarboxylic Acidsa DS

exptl

compound

pKa1

pKa2

pKa1

pKa2

succinic acid glutaric acid adipic acid

2.85 3.56 4.14

4.55 4.71 5.23

4.21 4.35 4.42

5.64 5.40 5.41

a

Computed and measured (ref 63) pKa’s of dicarboxylic acids.

TABLE 9: pKa’s of Dicarboxylic Acidsa DS

exptl

compound

pKaint

∆pKa

pKaint

∆pKa

succinic acid glutaric acid adipic acid

3.15 3.86 4.44

1.70 1.15 1.09

4.51 4.64 4.72

1.43 1.09 0.99

a Computed and measured pKa’s of dicarboxylic acids, expressed as intrinsic pKa and ∆pKa.

pKa Shifts in Small Dicarboxylic Acids. As an additional test, the DS method is used to compute the known pKa values of several dicarboxylic acids. The present treatment of the energetics of ionization is that used previously.12 Thus, a chemically similar “model” compound of known pKa is matched to each acidic group. In each case, the carboxylic acid with the appropriate number of carbons is used, i.e., propionic (pKa 4.87), butyric (pKa 4.82), and pentanoic adipic (pKa 4.84) acids for succinic, glutaric, and adipic acids, respectively. Symmetry numbers are accounted for.12 Ionization is modeled as a change in the partial charges of the atoms, with explicit addition or loss of a proton from the ionization site.12 The internal dielectric constant is set to 2 and that of the solvent to 78. The OPLS parameter set is used.50 Coulombic energy is calculated by pairwise summation, and the solvation energy is computed by the DS method with a grid spacing of 0.1 Å. The Friedman image-charge method is used to compute the energy contribution of the bulk of the solvent. Table 8 compares measured pKa values with those computed with the current method, and Table 9 presents the same results in the form of the intrinsic pKa values, pKaint, and interaction energies ∆pKa ) pKa1 - pKa2.12 The computations yield very good results for ∆pKa: the errors are less than 0.3 unit in all three cases. However, the computations overestimate the difference between pKaint and the model pKa’s. Previously published FD calculations for these molecules yielded somewhat better results for the intrinsic pKa’s.12 However, the present errors are well-behaved in the sense that they diminish smoothly with the length of the alkane chain separating the two acid groups. 4. Discussion and Conclusions The DS method yields electrostatic solvation energies for molecular systems that agree well with experiment and with the results of detailed solutions of the Poisson equation. This is despite the rather drastic approximations used in the method, and the modest use of parametrization. The chief approximation in the DS method is that the field generated by the solute is proportional to the Coulombic field that the solute would generate in a uniform dielectric medium. The constant of proportionality depends only upon the distribution of the charges. With this approximation, the shape of the boundary does not influence the field within the solvent. This contrasts with full solutions of the Poisson equation, for which the field within the solvent does depend upon the shape of the solute-solvent

boundary. Nonetheless, the shape of the dielectric boundary between solute and solvent does strongly influence the final results of the DS method, because the boundary determines the region where the effective polarizability of the medium changes, i.e., the volume of integration in eq 9. The success of the method suggests that this is the most important effect of the molecular shape upon electrostatic solvation energies. Still, there are applications for which the influence of molecular shape upon field must be included in order to reproduce experimental data.9 The present method resembles the noniterative LD models,36 which likewise use a screening factor in the solvent region and neglect changes of the field energy within the solute. However, the present method differs from the LD models in at least two important respects. First, the DS method assumes that the dielectric response is linear, whereas the LD model allows for dielectric saturation according to the Langevin function. Second, the DS method uses a closely spaced integration grid in order to model a dielectric continuum. In contrast, the LD method typically uses a grid that is coarse on the atomic size scale in order to model the granularity of the solvent. The accuracy of the DS method suggests that nonlinearity and granularity are not critical for the success of the LD models. It is worth noting that the DS method is quite different from the screened Coulomb potentials (SCP) that have been discussed previously,36,51,52,53 despite the similarity in name. The SCP models approximate the interaction energy of a pair of solvated or partly solvated charges with a function of distance that differs from the r-1 dependence of Coulomb’s law. The DS method offers no such distance-dependent screening function. In particular, it would be incorrect to use the screening factor of the DS method to compute charge-charge interactions. Rather, the interaction energy of two charges is the sum of their Coulombic interaction and a solvent-dependent volume integral that must be evaluated numerically. Such interactions are implicit in the reported ∆pKa values of dicarboxylic acids (see Results). The field energy (FE) method,54 like the present DS method, relies upon the assumption that the dielectric inhomogeneity of the system does not distort the electrostatic field. However, the DS method differs from the FE method in a number of respects. For example, the FE method treats self-energies and charge-charge interaction energies separately, but the DS method does not make this distinction. More interesting, perhaps, is that the DS method is based upon an integral of the energy density over the solvent region, but the FE method uses an integral over the interior of the solute. For small solutes, the FE method is likely to require less computer time. However, for larger solutes, such as proteins, the difference in computer time is likely to be small. This is because DS uses an analytic approximation for the contribution of bulk solvent, so the region occupied by the integration grid is about the same size as the region within the protein. An apparent drawback of integrating over the molecular interior is that the energy density varies sharply near solute charges. This can make it difficult to achieve accurate results and there is, in fact, relatively little information on the accuracy of the FE method.54 The DS model has only two adjustable parameters. These are optimized to reproduce FD results for seven small, organic molecules that are represented with the PARSE parameter set. However, the DS method yields good results with several different parameter sets, including PARSE, CHARMm22, and OPLS. Thus, the parametrization appears to be insensitive to the choice of atomic parameters. This supports the underlying validity of the method.

11236 J. Phys. Chem. B, Vol. 101, No. 51, 1997 Although the DS method is relatively accurate, it is not computationally efficient. Therefore, in its present form, it is not likely to be more useful than existing alternatives. However, what has been learned from this analysis can serve as a basis for the development of more efficient approximations. Acknowledgment. We thank Drs. D. Sitkoff, K. A. Sharp, and B. Honig for providing coordinates of the molecules used in training and testing the PARSE parameters and M. Potter for help with the coordinates and parameters used in the pKa calculations. This work was supported in part by NIH Grant GM41034 to J.M. and NIH Grant GM54053 to M.K.G. and by NIST and the ATP program of NIST. Certain commercial products are identified in this paper in order to specify the methods accurately. Such identification does not imply recommendation or endorsement by the National Institute of Standards and Technology, nor does it imply that the products identified are necessarily the best available for the purpose. References and Notes (1) Honig, B.; Sharp, K.; Yang, A.-S. J. Phys. Chem. 1993, 97, 11011109. (2) Sitkoff, D.; Sharp, K. A.; Honig, B. J. Phys. Chem. 1994, 98, 19781988. (3) Schmidt, A. B.; Fine, R. M. Mol. Simul. 1994, 13, 347-365. (4) Sharp, K. A.; Nicholls, A.; Friedman, R.; Honig, B. Biochemistry 1991, 30, 9686-9697. (5) Gilson, M. K.; Rashin, A. A.; Fine, R.; Honig, B. J. Mol. Biol. 1985, 183, 503-516. (6) Gilson, M. K.; Honig, B. Proteins: Struct., Funct., Genet. 1988, 3, 32-52. (7) Gilson, M. K.; Honig, B. Proteins: Struct., Funct., Genet. 1988, 4, 7-18. (8) Warwicker, J.; Watson, H. C. J. Mol. Biol. 1982, 157, 671-679. (9) Klapper, I.; Hagstrom, R.; Fine, R.; Sharp, K.; Honig, B. Proteins: Struct., Funct., Genet. 1986, 1, 47-79. (10) Gilson, M. K.; Sharp, K. A.; Honig, B. H. J. Comput. Chem. 1988, 9, 327-335. (11) McQuarrie, D. A. Statistical Mechanics; Harper & Row: New York, 1973. (12) Potter, M. J.; Gilson, M. K.; McCammon, J. A. J. Am. Chem. Soc. 1994, 116, 10298-10299. (13) Rajasekaran, E.; Jayaram, B.; Honig, B. J. Am. Chem. Soc. 1994, 116, 8238-8240. (14) Bashford, D.; Karplus, M. Biochemistry 1990, 29, 10219-10225. (15) Langsetmo, K.; Fuchs, J. A.; Woodward, C. Biochemistry 1991, 30, 7603-7609. (16) Oberoi, H.; Allewell, N. M. Biophys. J. 1993, 65, 48-55. (17) Yang, A.-S.; Gunner, M. R.; Sampogna, R.; Sharp, K.; Honig, B. Proteins: Struct., Funct., Genet. 1993, 15, 252-265. (18) Antosiewicz, J.; McCammon, J. A.; Gilson, M. K. J. Mol. Biol. 1994, 238, 415-436. (19) Antosiewicz, J.; McCammon, J. A.; Gilson, M. K. Biochemistry 1996, 35, 7819-7833. (20) Demchuk, E.; Wade, R. C. J. Phys. Chem. 1996, 100, 1737317387. (21) Alexov, E. G.; Gunner, M. R. Biophys. J. 1997, 74, 2075-2093. (22) Gilson, M. K. In Computer Simulations of Biomolecular Systems; Escom: Amsterdam, 1997. (23) Zauhar, R. J.; Morgan, R. S. J. Mol. Biol. 1985, 186, 815-820. (24) Juffer, A. H.; Botta, E. F. F.; van Keulen, B. A. M.; van der Ploeg, A.; Berendsen, H. J. C. J. Comput. Phys. 1991, 97, 144-171. (25) You, T. J.; Harvey, S. C. J. Comput. Chem. 1993, 14, 484-501.

Luo et al. (26) Nicholls, A.; Honig, B. J. Comput. Chem. 1991, 12, 435-445. (27) Holst, M.; Kozack, R. E.; Saied, F.; Subramaniam, S. Proteins: Struct., Funct., Genet. 1994, 18, 231-245. (28) Friedman, H. L. Mol. Phys. 1975, 29, 1533-1543. (29) Abagyan, R.; Totrov, M. J. Mol. Biol. 1994, 235, 983-1001. (30) Gilson, M. K.; Davis, M. E.; Luty, B. A.; McCammon, J. A. J. Phys. Chem. 1993, 97, 3591-3600. (31) Horvath, D.; Lippens, G.; van Belle, D. J. Chem. Phys. 1996, 105, 4197-4210. (32) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. Chem. Phys. Lett. 1995, 246, 122-129. (33) Qiu, D.; Shenkin, P. S.; Hollinger, F. P.; Still, W. C. J. Phys. Chem. 1997. (34) Davis, M. E. J. Chem. Phys. 1994, 100, 5149-5158. (35) David, L.; Field, M. J. Chem. Phys. Lett. 1995, 245, 371-376. (36) Warshel, A.; Russell, S. T. Quantum ReV. Biophys. 1984, 17, 285422. (37) Lee, F. S.; Chu, Z. T.; Warshel, A. J. Comput. Chem. 1993, 14, 161-185. (38) Warshel, A.; Chu, Z. T. In Cramer, C. J., Truhlar, D. G., Eds., Structure and ReactiVity in Aqueous Solution. Characterization of Chemical and Biological Systems; ACS Symposium Series; American Chemical Society: Washington, D. C., 1994; pp 71-94. (39) Bo¨ttcher, C. J. F.; van Belle, O. C.; Bordewijk, P.; Rip, A. Theory of Electric Polarization, 2nd ed.; Elsevier Scientific Publishing Company: Amsterdam, 1973; Vol. 1. (40) Arfken, G. Mathematical Methods for Physicists, 3rd ed.; Academic Press, Inc.: San Diego, CA, 1985. (41) Davis, M. E.; Madura, J. D.; Luty, B. A.; McCammon, J. A. Comput. Phys. Commun. 1991, 62, 187-197. (42) Davis, M. E.; McCammon, J. A. J. Comput. Chem. 1991, 12, 909912. (43) Sitkoff, D.; Sharp, K. A.; Honig, B. J. Phys. Chem. 1994, 98, 19781988. (44) Bernstein, F. C.; Koetzle, T. F.; Williams, T. F.; Meyer, G. J. B., Jr.; Brice, M. D.; Rodgers, J. R.; Kennard, O.; Shimanouchi, T.; Tasumi, M. J. Mol. Biol. 1977, 112, 535-542. (45) Bolin, J. T.; Filman, D. J.; Matthews, D. A.; Hamlin, R. C.; Kraut, J. J. Biol. Chem. 1982, 257, 13650. (46) Moult, J.; Sussman, F.; James, M. N. G. J. Mol. Biol. 1985, 182, 555. (47) Brunger, A. T.; Karplus, M. Proteins: Struct., Funct., Genet. 1988, 4, 148-156. (48) QUANTA; Molecular Simulations, Inc.: Waltham, MA. (49) CHARMm, Version 22; Molecular Simulations, Inc.: Waltham, MA, 1992. (50) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657-1666. (51) Hingerty, B. E.; Ritchie, R. H.; Ferrel, T. L.; Turner, J. E. Biopolymers 1985, 24, 427. (52) Lavery, R.; Sklenar, H.; Zakrzewska, K.; Pullman, B. J. Biomol. Struct. Dyn. 1986, 3, 989. (53) Mehler, E. L.; Eichele, G. Biochemistry 1984, 23, 3887-3891. (54) Schaefer, M.; Froemmel, C. J. Mol. Biol. 1990, 216, 1045-1066. (55) Cabani, S.; Gianni, P.; Mollica, V.; Lepori, L. J. Solution Chem. 1978, 7, 721-770. (56) Wolfenden, R.; Andersson, L.; Cullis, P. M.; Southgate, C. C. Biochemistry 1981, 20, 849. (57) Wolfenden, R. Biochemistry 1978, 17, 199. (58) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. J. Am. Chem. Soc. 1990, 112, 6127-6129. (59) Davis, M. E. J. Chem. Phys. 1994, 100, 5149-5159. (60) David, L.; Field, M. J. Chem. Phys. Lett. 1990, 90, 509-521. (61) Cramer, C. J.; Truhlar, D. G. J. Am. Chem. Soc. 1991, 113, 83058311. (62) Tannor, D. J.; Marten, B.; Murphy, R.; Friesner, R. A.; Sitkoff, D.; Nicholls, A.; Ringnalda, M.; Goddard, W. A. I.; Honig, B. J. Am. Chem. Soc. 1994, 116, 11875-11882. (63) Martell, A. E.; Smith, R. M. Critical Stability Constants; Plenum Press, Inc.: New York, 1974; Vol. 1-6.

Dielectric Screening Treatment of Electrostatic Solvation - The Journal

Recommend Documents