4404
J. Phys. Chem. B 1998, 102, 4404-4410
Parameter Dependence in Continuum Electrostatic Calculations: A Study Using Protein Salt Bridges Zachary S. Hendsch, Charles V. Sindelar, and Bruce Tidor* Department of Chemistry, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139-4307 ReceiVed: September 3, 1997; In Final Form: February 4, 1998
Continuum electrostatic calculations were carried out for a data set of 21 salt bridges using three sets of parameters (CHARMM, PARSE, and a version of OPLS with slightly exaggerated atomic radii) and three representations of the dielectric boundary (approximate, exact, and smoothed versions of the analytic molecular surface). Comparisons of the overall results and a decomposition into terms due to desolvation, the direct salt bridge, and other interactions with the surrounding protein were used to analyze the effect of differences in parameter sets. It is concluded that (1) while numerical differences among the results can be significant, the underlying tension between desolvation penalty incurred and favorable interactions recovered upon folding is similar in all computations; (2) the largest variation occurs for the desolvation penalty, indicating the importance of fitting this property in parameter development; (3) the approximation of the molecular surface produces the largest desolvation penalties, whereas the smoothed boundary produces the smallest; (4) the exaggerated-radii OPLS implementation used here resulted in underestimates of the desolvation penalty and overestimates of the bridge interaction, relative to the other parameter sets, for 15 of 21 salt bridges; and (5) although the origins and parameter values are different for the CHARMM and PARSE parameter sets, they produce rather similar energetics, due in part to compensating effects.
Introduction The use of continuum electrostatic calculations is becoming a standard tool to analyze properties of proteins and other biological macromolecules (for reviews see refs 1-4). Studies reported in the literature differ from one another in many details of how these methods are applied, including the use of different parameter sets (partial atomic charges and atomic radii), values for the dielectric constant for the biomolecule interior (with 2, 3, and 4 being the most common), definitions of the molecular surface, treatment of the dielectric and ionic-strength boundary between interior and exterior of the biomolecule, boundary conditions (generally at large distance from the molecule), and resolution of the grid used for numerical solution of the Poisson-Boltzmann equation or its linearized form. Because of the increased popularity of these methods and the large number of potential variations in their application, it is important to perform systematic studies of how results depend on these variables. Here we investigate the dependence of the stability effect of salt-bridging side chains on variations in the atomic parameters and the treatment of the dielectric change at the molecular boundary. Specifically, we report the results of investigations using three sets of parameters (atomic radii and partial atomic charges) in common use (CHARMM,5 PARSE,6 and OPLS;7 the OPLS parameters used here purposely have slightly exaggerated radii) and three representations of the dielectric boundary (the DELPHI version 3.0 approximation to the analytic molecular surface, the exact analytic molecular surface,8,9 and the exact analytic molecular surface with a modified smoothing algorithm applied10,11). We tested the effect of changing these parameters on the calculated electrostatic effect of salt bridges on protein stability using a data set from our earlier work.12 Because we find negligible differences between the solutions to the full nonlinear Poisson-Boltzmann equation and its linearized form, we
employed the latter and exploited the linearity to compute three additive contributions to the overall stability effect of each salt bridge: (1) ∆∆Gsolv is the electrostatic penalty for desolvating the two charged side chains in protein folding, (2) ∆∆Gbridge is the favorable electrostatic interaction between bridging side chains in the folded state, and (3) ∆∆Gprotein is the (generally favorable) electrostatic interaction between members of a salt bridge and polar and charged groups in their environment.12 An important aspect of the work here is that we have compared not only the overall electrostatic effect, ∆∆Gtotal, but also the individual contributions to understand better the effects of differences among these parameter sets and boundary treatments. Methods All calculations were performed using a modified version of the DELPHI computer program1,13,14 as described previously (with an internal dielectric of 4 and an ionic strength of 145 mM) except that ∆∆Gsolv was calculated as the difference in total energy here, as opposed to reaction-field energy in our previous work, between folded- and unfolded-state models.12 The linearized Poisson-Boltzmann equation was solved here, though very similar results were obtained for test calculations using the full nonlinear form. The uncertainty due to the granularity of the grid was estimated as 2 times the standard error for 10 translations on the finite-differences grid.15 Other details of the calculations were identical to our previous work unless stated otherwise.12 CHARMM parameters from Karplus and co-workers were from the PARAM19 parameter set.5 Rmin/2 was taken as the atomic radius for each atom, including hydrogen. OPLS parameters were from Jorgensen et al.7 The atomic radii were taken as 2-5/6σ, which is equivalent to Rmin/2; hydrogen atoms had an atomic radius of 1.25 Å. PARSE parameters from Honig and co-workers were adapted from the literature with the following
S1089-5647(97)02866-6 CCC: $15.00 © 1998 American Chemical Society Published on Web 05/07/1998
Continuum Electrostatic Calculations adjustments and additions.6 Charges of +0.125 at Cβ and -0.125 at Cγ were used for Phe, Tyr, Trp, and neutral and charged forms of His. This choice was based on the charge distribution for toluene. Charges of 0.000 were used for the aromatic carbons at Cδ2 and C2 in Trp. CHARMM charges were used for Pro and for disulfide-bonded Cys. Charges of 0.275 were used for Cδ2 and C1 of the charged His side chain. Charges for the N- and C-terminus were transferred from Lys and Glu side chains. It should be noted that the CHARMM and OPLS parameter sets were developed to produce agreement between experiment and structural, energetic, and thermodynamic properties calculated from molecular dynamics and Monte Carlo simulations, but the PARSE parameters were fit to solvation free energies of small molecules calculated with continuum electrostatic methods. Thus, the PARSE parameter set might be expected to perform best for continuum electrostatic calculations. The OPLS parameter set was designed with a hydrogen atom van der Waals radius of 0.00 Å. Here we inflated the value to 1.25 Å to purposefully create a parameter set that would underestimate desolvation penalties, particularly for positively charged side chains. Mohen et al. have shown that the OPLS parameter set with either atomic radii of 2-1σ and hydrogen radii of 1.25 Å or with atomic radii of 2-5/6σ and hydrogen radii of 0.00 Å produced reasonably good agreement with experimentally determined solvation energies for small polar and charged molecules, but that use of the larger value for both hydrogen and non-hydrogen radii tended to underestimate the magnitude of aqueous solvation energies.11 Thus, the current work attempts to understand how differences in parameters are propagated into differences in the types of energies calculated in continuum calculations, rather than as an evaluation of individual parameter sets, particularly because we intentionally chose a nonoptimal set of OPLS-based parameters. In their development, the CHARMM and OPLS parameter sets used a dielectric constant of 1 and PARSE used a dielectric constant of 2. Here an interior dielectric constant of 4 was used to account for the larger-scale atomic and dipolar fluctuations expected for a protein compared to a small molecule. The exterior dielectric was 80. The method for generating the dielectric boundary surface in DELPHI version 3.0 was an approximation generated by assigning grid lines whose midpoint fell within the solvent accessible surface the interior dielectric, finding a mesh of points that lie on this solvent accessible surface, and then reducing the solvent accessible surface to the molecular surface by changing the dielectric constant assigned to grid-line midpoints within a one-water-molecule radius of a “mesh point” back to the dielectric constant of the exterior.13 This algorithm fails when the number of mesh points is insufficient to define the reentrant surface of the molecule. The surface generation method has been improved in more recent versions of the program.6 We have written our own program for determining the intersections of the analytic molecular surface (using the framework of Connolly9,16) with the grid lines used in the finitedifference Poisson-Boltzmann calculation and employed it here. In standard continuum electrostatic implementations, the midpoints of grid lines interior to the dielectric boundary are assigned one value for the dielectric and those exterior are assigned another; there are no intermediate values. This method was employed for both the DELPHI v3.0 and the analytic molecular surface here. In additional calculations, the modified smoothing algorithm of McCammon, Pettit, and co-workers was used with the analytic molecular surface.10,11 In this method, a grid line that intersects the dielectric boundary is assigned a
J. Phys. Chem. B, Vol. 102, No. 22, 1998 4405 TABLE 1: Parameter Dependence of Salt-Bridge Stabilitya surface representation parameter set
DELPHI
v3.0
analytic
with modified smoothing
CHARMM PARSE OPLS
∆∆Gtotal (Total Contribution) 4.78 (0.21) 2.86 (0.15) 3.14 (0.13) 2.61 (0.11) 0.84 (0.08) 0.39 (0.05)
1.93 (0.05) 1.84 (0.04) -0.23 (0.05)
CHARMM PARSE OPLS
∆∆Gsolv (Desolvation Penalty) 13.02 (0.23) 10.88 (0.17) 11.29 (0.14) 10.33 (0.11) 10.50 (0.10) 9.70 (0.08)
10.35 (0.05) 9.96 (0.03) 9.48 (0.02)
CHARMM PARSE OPLS
∆∆Gbridge (Bridge Term) -5.05 (0.05) -4.93 (0.03) -4.88 (0.04) -4.61 (0.03) -6.27 (0.05) -6.02 (0.04)
-5.22 (0.03) -4.87 (0.03) -6.33 (0.04)
CHARMM PARSE OPLS
∆∆Gprotein (Environment Term) -3.18 (0.03) -3.10 (0.02) -3.26 (0.03) -3.11 (0.02) -3.39 (0.03) -3.29 (0.03)
-3.19 (0.02) -3.25 (0.02) -3.38 (0.02)
a All values in kcal/mol. Positive ∆∆G values indicate terms unfavorable to salt-bridge formation. The average error for each value, estimated as twice the standard deviation of the mean, is listed in parentheses.
dielectric value between the interior and exterior value determined by the equation
)
extint aint - (1 - a)ext
(1)
where ext is the dielectric constant for the solvent, int is the dielectric constant for the solute, and a is the fraction of the grid line corresponding to solvent.11 Results Overall Averages. The overall averaged numerical results of this study are presented in Table 1. Positive ∆∆G values indicate terms unfavorable for salt-bridge formation. ∆∆Gtotal shows that CHARMM and PARSE parameters produce similar overall results (more so using the analytic surface representation, with or without modified smoothing), and results with OPLS are quite different and underestimate the electrostatic penalty of forming salt bridges relative to the other two. Also apparent is the trend that the DELPHI v3.0 approximate surface produces the most destabilizing salt-bridge electrostatics, the analytic surface reduces this by 20-50%, and addition of modified smoothing reduces the average value further still, for all three parameter sets. The overall result is that salt bridges are calculated to be electrostatically destabilizing by 2-5 kcal/mol for either the CHARMM or PARSE parameters over the range of dielectric boundaries used, whereas the exaggerated OPLS parameters find little or no electrostatic penalty (or benefit) over the same range. It should be noted that the comparison of overall average results, while useful because it allows general similarities and differences to be detected, tends to mask the individual variability seen between parameter sets (see below). Whereas the difference between CHARMM and PARSE results using the analytic molecular surface is 0.25 kcal/mol when the overall averages are compared, the rms difference between the individual results is 1.1 kcal/mol. To investigate further the source of differences among these parameter sets and surface representations, Table 1 also presents the three contributions ∆∆Gsolv, ∆∆Gbridge, and ∆∆Gprotein. Interestingly, all nine variations give very similar values for
4406 J. Phys. Chem. B, Vol. 102, No. 22, 1998
Hendsch et al.
Figure 2. Comparison of the atomic charge parameters for the different parameter sets used in this study. (a) CHARMM vs PARSE parameters, (b) CHARMM vs OPLS parameters, and (c) PARSE vs OPLS parameters. Figure 1. Comparison of the radius parameters, in Å, for the different parameter sets used in this study. (a) CHARMM vs PARSE parameters, (b) CHARMM vs OPLS parameters, and (c) PARSE vs OPLS parameters.
∆∆Gprotein. Since many of the bridges are largely buried, one would expect this term to be more sensitive to the magnitudes of charges than to the sizes of the atomic radii or the definition of the molecular boundary. Thus, it is not surprising that this term is relatively insensitive to changes in the treatment of the molecular boundary. CHARMM and PARSE give similar values for ∆∆Gbridge that are less favorable than those given by OPLS. This result could be consistent with the OPLS radii and/or charges being of somewhat larger magnitude than for the other parameter sets. Increased charge magnitudes obviously increase the interaction between oppositely charged groups; likewise, increased radii might be expected to decrease the effective dielectric constant for the bridges. Finally, the solvation penalties reveal the greatest discrepancy among parameters, with results using CHARMM being larger than PARSE, which are larger than OPLS, and with the approximate DELPHI v3.0 surface being larger than the analytic molecular surface, which is larger than when modified smoothing was applied.
Comparison of Parameter Values. A comparison of the atomic radii from the three parameter sets is given in Figure 1. This shows that (a) the CHARMM radii are uniformly larger than the corresponding PARSE values except for two atom types, which are both hydrogen, that have larger PARSE values; (b) the hydrogen radii are larger in the OPLS set than in CHARMM, but the non-hydrogen radii are roughly similar between these sets (though the smaller atoms are somewhat larger in OPLS and the larger atoms tend to be larger in CHARMM, resulting in a larger overall range for the CHARMM non-hydrogen radii; one anomaly is the Arg Cζ, which has a radius of 2.1, 1.7, and 1.263 Å in CHARMM, PARSE, and OPLS parameter sets, respectively); and (c) the OPLS radii are uniformly larger than the PARSE radii, except for the single exception of Arg Cζ, which lies below the line in Figure 1c. It should be noted that the CHARMM parameter set used contains only strongly polar hydrogen atoms, the OPLS parameter set also includes the thiol hydrogen of Cys, and the PARSE parameter set further includes aromatic hydrogen atoms. A similar comparison of the charge sets is given in Figure 2. This shows that (a) while the CHARMM and PARSE charge sets are clearly distinct, there does not appear to be an obvious
Continuum Electrostatic Calculations
Figure 3. Comparison of ∆∆Gtotal calculated using different parameter sets with the analytic surface. All values are in kcal/mol. (a) CHARMM vs PARSE parameters, (b) CHARMM vs OPLS parameters, and (c) PARSE vs OPLS parameters.
systematic difference between them (i.e., one is not much more polarized than the other); and (b) the OPLS charge set tends to be more polarized than the CHARMM and PARSE sets, with positive charges tending to be more positive and negative charges more negative. Detailed Free Energy Results. Figures 3-6 compare the values of ∆∆Gtotal, ∆∆Gsolv, ∆∆Gbridge, and ∆∆Gprotein calculated for each individual salt bridge using the three different parameter sets with the analytic surface description. There is good general agreement between CHARMM and PARSE for ∆∆Gtotal, with little apparent systematic difference between the two. Nevertheless, the rms difference between values calculated with the two methods is 1.1 kcal/mol, and the largest discrepancy is 2.6 kcal/ mol. The OPLS parameters used here give values for ∆∆Gtotal that are systematically smaller than those from CHARMM or PARSE, with discrepancies of up to 7.2 kcal/mol. The parameter dependence for the cost of desolvating the side chains (∆∆Gsolv) is examined in Figure 4. As expected, the OPLS parameters tend to give the smallest ∆∆Gsolv of the three parameter sets, though not uniformly so. Indeed, 6 of
J. Phys. Chem. B, Vol. 102, No. 22, 1998 4407
Figure 4. Comparison of ∆∆Gsolv calculated using different parameter sets with the analytic surface. All values are in kcal/mol. (a) CHARMM vs PARSE parameters (b) CHARMM vs OPLS parameters, and (c) PARSE vs OPLS parameters.
the 21 salt bridges actually have a larger desolvation penalty with the OPLS parameters than with the CHARMM; thus, while the OPLS parameters used here produce a smaller desolvation cost on average, this may not be the case for any particular example. Interestingly, the OPLS parameters give a larger desolvation penalty than do the CHARMM for each negatively charged side chain, but the reverse is true for all but one positively charged side chain. This is due to larger charges in OPLS that are offset by larger hydrogen radii for the positively charged side chains, but not the negatively charged ones (which have no polar hydrogens). The OPLS parameters also result in a ∆∆Gbridge term that is approximately 1 kcal/mol more stabilizing than the other two parameter sets (using the analytic surface; see Table 1 and Figure 5). This difference appears relatively systematic (see Figure 5). The cause of the difference in ∆∆Gbridge was investigated by repeating the calculation with a hybrid parameter set consisting of CHARMM charges and OPLS radii. This gave an average ∆∆Gbridge of 5.45 kcal/mol, which is about halfway
4408 J. Phys. Chem. B, Vol. 102, No. 22, 1998
Hendsch et al.
Figure 5. Comparison of ∆∆Gbridge calculated using different parameter sets with the analytic surface. All values are in kcal/mol. (a) CHARMM vs PARSE parameters, (b) CHARMM vs OPLS parameters, and (c) PARSE vs OPLS parameters.
Figure 6. Comparison of ∆∆Gprotein calculated using different parameter sets with the analytic surface. All values are in kcal/mol. (a) CHARMM vs PARSE parameters, (b) CHARMM vs OPLS parameters, and (c) PARSE vs OPLS parameters.
between the averages calculated with the full CHARMM and OPLS parameters. This indicates that the difference in ∆∆Gbridge is due both to the larger OPLS radii (which decreases the effective dielectric between the salt-bridging side chains) and to the more polarized OPLS charges. The average ∆∆Gprotein for the 21 salt bridges shows very little dependence on parameters (see Table 1). A comparison of ∆∆Gprotein for each salt bridge with the analytic surface (see Figure 6) shows very little variation between the CHARMM and OPLS parameters (except for one outlier) and moderate variation when comparing PARSE with the other two. It is surprising that the OPLS and CHARMM results are more similar than the CHARMM and PARSE results, since the latter’s averages are nearly identical, but Figure 6b shows that OPLS is systematically more negative than CHARMM, which has a cumulative effect on the difference of the average. The one outlier between the CHARMM and OPLS parameters, a salt bridge between Glu35 and Arg10 from the N-terminal domain of 434 repressor, was examined in more detail. ∆∆Gprotein is more favorable by 2.3 kcal/mol with the OPLS
parameters, due mainly to an increase in the interactions between Glu35 and two surface residues that bury it, Lys7 and Gln17 (whose favorable interaction with Glu35 improves by 1.8 and 0.9 kcal/mol, respectively). To test whether this difference is due to the increased radii, which would lessen the solvent screening of these interactions, or to the larger partial charges in the OPLS parameter set, the calculation was redone with a hybrid parameter set consisting of CHARMM charges and OPLS radii and a second hybrid set with OPLS charges and CHARMM radii. For Lys7 the largest effect on the magnitude of the interaction is from the difference in atomic radii. Using CHARMM charges with the OPLS radii improves the interaction between Glu35 and Lys7 by 1.4 kcal/mol over using all CHARMM parameters, whereas using OPLS charges with CHARMM radii results in an increase of only 0.3 kcal/mol. The improvement in the interaction between Gln17 and Glu35 is due in about equal parts to the change in charges and the change in radii. Both CHARMM charges with OPLS radii and OPLS charges with CHARMM radii increase the interaction between Gln17 and Glu35 by 0.4 kcal/mol.
Continuum Electrostatic Calculations
Figure 7. Comparison of ∆∆Gtotal and its components calculated using the analytic surface and the smoothed analytic surface. All values are in kcal/mol and were calculated using the CHARMM parameters. (a) ∆∆Gtotal, (b) ∆∆Gsolv, (c) ∆∆Gbridge, and (d) ∆∆Gprotein.
Interestingly, the modified smoothing algorithm applied to the analytic surface produces a ∆∆Gtotal that is less destabilizing (see Table 1 and Figure 7). This is caused by a decrease in ∆∆Gsolv and an increase in ∆∆Gbridge, in parallel with what would be expected with increased atomic radii. Discussion Here we have carried out continuum electrostatic calculations to study the effect of 21 salt bridges on protein stability using three different parameter sets and three representations of the dielectric boundary. One of the parameter sets was an interpretation of the OPLS model intentionally chosen to underestimate the desolvation penalty. The range of results for the average over 21 salt bridges is -0.23 to 4.78 kcal/mol. The PARSE and CHARMM sets alone (neglecting the OPLS set with purposefully inflated radii) covered a range of 1.84-4.78 kcal/ mol, and when the approximate surface representation is ignored, the range narrows to 1.84-2.86 kcal/mol, with the difference arising more from whether modified smoothing is applied to the boundary than from differences in the radius and charge sets. Interestingly, this is true even though significant differences do exist between CHARMM and PARSE radius and charge
J. Phys. Chem. B, Vol. 102, No. 22, 1998 4409 sets, with the PARSE radii smaller for all atoms but hydrogen. This is particularly remarkable given the rather different origins and purposes of these two sets. Furthermore, we note that our analytic surface representation and the new DELPHI surface representation give very similar results (data not shown). Importantly, compensating changes occur between parameter sets that moderate differences in the computed energetics. For example, the two largest differences in the ∆∆Gbridge term (which represents a single side chain-side chain interaction) between the CHARMM and PARSE parameter sets are 4.32 kcal/ mol for Arg10-Glu35 in 434 repressor (-10.56 vs -14.88 kcal/ mol) and 2.52 kcal/mol for Asp10-Arg148 in T4 lysozyme (-8.80 vs -11.32 kcal/mol), with the stronger interaction computed using PARSE parameters in both cases. In each case compensation by differences in ∆∆Gsolv and ∆∆Gprotein terms results in rather similar computed differences between the parameter sets for the total effect of the salt bridge; the differences in ∆∆Gtotal are 0.18 kcal/mol in 434 repressor and 0.44 kcal/mol in T4 lysozyme. This is the type of effect that one would expect to see if one parameter set overestimated interaction strengths with solvent and with other protein groups relative to another parameter set. While such parameter uncertainties can clearly affect some of the computed values, they do not, in general, affect the overall interpretation one would draw from the computations. Here for instance, the conclusion that salt bridges do not provide substantial electrostatic stabilization for proteins due to a large desolvation penalty that is rarely fully recovered in favorable electrostatic interactions in the folded state, is common to all nine sets of calculations performed. Moreover, the use of exaggerated radii for the OPLS computations further strengthens these conclusions, since this parameter set would tend to underestimate the desolvation penalty and overestimate the favorable electrostatic interactions on average, yet the net effect was still not found to be strongly stabilizing. One further observation regarding the OPLS parameter set is that the charges tend to be more polarized than the CHARMM and PARSE parameter sets. The largest variation in the results presented here is due to variation in computations of the desolvation penalty. Thus, this critical property should be given high priority in the development of parameters for continuum electrostatic models. That the PARSE parameter set was fit to solvation energies suggests that it should be an especially accurate parameter set for a variety of folding and binding studies.6 In summary, parallel calculations with a variety of parameter sets and surface descriptions using a data set of 21 salt bridges have documented a relatively narrow range of variability in the results and elucidated a number of important trends. Moreover, the use of a thermodynamic decomposition based on the linearized Poisson-Boltzmann equation has allowed overall patterns in the computed energetics to be traced to differences in the parameter sets. This information is expected to be useful not only in the development and application of parameter sets but also in molecular design, where the effect of functional group substitutions (sulfhydryl for hydroxyl, for example) could be estimated through the appropriate changes of partial atomic charge, geometry, and radii. (In the PARSE parameter set S/H (O/H) charges (e) are -0.29/+0.29 (-0.49/+0.49) and radii (Å) are 1.85/1.00 (1.40/1.00).6) The decomposition involves three terms: desolvation, direct salt-bridging interaction, and interactions with other protein groups. While effects of changing charge parameters are as anticipated (larger charge magnitudes lead to larger desolvation penalties and stronger interactions), an unanticipated result is that larger radii not only
4410 J. Phys. Chem. B, Vol. 102, No. 22, 1998 decrease desolvation penalties but also enhance interactions through enlarging the low dielectric region and thereby leading to lowered effective dielectrics. Because the effect of altering charge tends to cancel between desolvation and interaction terms, it may not lead to as strong effects as changes to atomic radii, whose effect on desolvation reinforces that of interaction. This could be exploited in design studies through the selective replacement (enlargement) of groups on the surface involved either directly in electrostatic interactions or indirectly, by burying other interacting polar or charged groups. Acknowledgment. We thank Barry Honig for making the computer program available and Barry Honig and Richard A. Friesner for helpful discussions. This work was supported in part by the National Institutes of Health (GM47678 and GM55758).
DELPHI
References and Notes (1) Sharp, K. A.; Honig, B. Annu. ReV. Biophys. Biophys. Chem. 1990, 19, 301.
Hendsch et al. (2) Davis, M. E.; McCammon, J. A. Chem. ReV. 1990, 90, 509. (3) Honig, B.; Sharp, K.; Yang, A.-S. J. Phys. Chem. 1993, 97, 1101. (4) Honig, B.; Nicholls, A. Science (Washington, D.C.) 1995, 268, 1144. (5) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (6) Sitkoff, D.; Sharp, K. A.; Honig, B. J. Phys. Chem. 1994, 98, 1978. (7) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657. (8) Richards, F. M. Annu. ReV. Biophys. Bioeng. 1977, 6, 151. (9) Connolly, M. L. J. Appl. Crystallogr. 1983, 16, 548. (10) Davis, M. E.; McCammon, J. A. J. Comput. Chem. 1991, 7, 909. (11) Mohan, V.; Davis, M. E.; McCammon, J. A.; Pettitt, B. M. J. Phys. Chem. 1992, 96, 6428. (12) Hendsch, Z. S.; Tidor, B. Protein Sci. 1994, 3, 211. (13) Gilson, M. K.; Sharp, K. A.; Honig, B. H. J. Comput. Chem. 1988, 9, 327. (14) Gilson, M. K.; Honig, B. Proteins: Struct., Funct., Genet. 1988, 4, 7. (15) Taylor, J. R. An Introduction to Error Analysis; University Science Books: Mill Valley, CA, 1982. (16) Connolly, M. L. Science (Washington, D.C.) 1983, 221, 709.
