Introducing Charge Hydration Asymmetry into the Generalized Born Model.

The effect of charge hydration asymmetry (CHA)-non-invariance of solvation free energy upon solute charge inversion-is missing from the standard linear response continuum electrostatics. The proposed charge hydration asymmetric-generalized Born (CHA-GB) approximation introduces this effect into the popular generalized Born (GB) model. The CHA is added to the GB equation via an analytical correction that quantifies the specific propensity of CHA of a given water model; the latter is determined by the charge distribution within the water model. Significant variations in CHA seen in explicit water (TIP3P, TIP4P-Ew, and TIP5P-E) free energy calculations on charge-inverted "molecular bracelets" are closely reproduced by CHA-GB, with the accuracy similar to models such as SEA and 3D-RISM that go beyond the linear response. Compared against reference explicit (TIP3P) electrostatic solvation free energies, CHA-GB shows about a 40% improvement in accuracy over the canonical GB, tested on a diverse set of 248 rigid small neutral molecules (root mean square error, rmse = 0.88 kcal/mol for CHA-GB vs 1.24 kcal/mol for GB) and 48 conformations of amino acid analogs (rmse = 0.81 kcal/mol vs 1.26 kcal/mol). CHA-GB employs a novel definition of the dielectric boundary that does not subsume the CHA effects into the intrinsic atomic radii. The strategy leads to finding a new set of intrinsic atomic radii optimized for CHA-GB; these radii show physically meaningful variation with the atom type, in contrast to the radii set optimized for GB. Compared to several popular radii sets used with the original GB model, the new radii set shows better transferability between different classes of molecules.

Introduction

The conceptual simplicity of the implicit solvation (solvent treated as a structureless continuum dielectric medium) framework[1−4] facilitates fast quantitative estimates of solvation effects that are key in biomolecular computation. In this framework, one often approximates the solvation free energy—free energy change to transfer a solute from gas phase to aqueous phase—as a sum of its electrostatic (polar) and non-polar parts, that is, ΔGsolv = ΔGpol + ΔGnp. Among these, ΔGpol is by far the largest contribution in the majority of molecules; ΔGnp has not been the computational bottleneck so far, perhaps due to the simplistic nature of the popular approximations most commonly used in this context.[5] The defacto “workhorse” for estimation of ΔGpol in most practical computations is the linear response continuum treatment of electrostatics based on numerical solutions of the Poisson (PE) or Poisson–Boltzmann (PB) equation,[3] or the generalized Born (GB) model,[4] which is an approximation to the PE. Because of its closed analytical form, the GB is arguably the fastest among these implicit solvent models. Some noteworthy applications of GB include protein folding,[6−15] “large scale” motions in biomolecules,[16−18] analysis of nucleic acid conformational energetics,[19,20] binding between proteins and nucleic acids,[21−23] modeling of peptides in membrane environment,[24−28] and many others.[29−32]

A number of approximations of various severity separate the linear response continuum PE/PB or GB models from the more accurate explicit solvent treatment.[33] One such serious approximation is the neglect of charge hydration asymmetry (CHA)[34−41]—the characteristic dependence of solvation free energy, ΔGsolv, on the sign of the charges in a solute. The effect is clearly missing from the linear response continuum model, in which ΔGsolv ∼ q2/R for a single spherical charge q of radius R. A straightforward example of CHA is the large difference between experimental ΔGsolv for oppositely charged ions of similar ionic radius; for example, K+ and F–[42,43] differ roughly by 50 kcal/mol, which is about 50% of their individual ΔGsolv values.

Within the continuum framework, experimental solvation energies of single-charge spherical ions could be reproduced accurately by ad hoc adjustments made to the ion radii.[44,45] These adjustments essentially amount to redefining the dielectric boundary, which (to an extent) implicitly accounts for the missing CHA by making the boundary depend on the ion charge. While for single ions the procedure works well, its generalization to the multi-atom case did not yet lead to a single consensus[46] set of atomic radii; a set that works well for a certain class of molecules does not necessarily perform equally well for a different class of structures (see, for example, Table 1). We argue that to a significant degree the problem is the missing CHA; we show that once the effect is explicitly added into even the simplest of the continuum models, the GB, a more universal set of atomic radii can be developed. Several existing “beyond PE” semi-empirical approaches such as the semi-explicit assembly (SEA)[51] (assembling hydration shells of solutes using the precomputed properties of hydration shells of charged spheres in explicit solvent) or integral equations based reference interaction site models (RISM, 3D-RISM)[52−55] already implicitly account for CHA, along with many other effects missing from the linear response continuum framework. However, these methods are not comparable to the GB model in computational efficiency and conceptual simplicity, which gives us an additional motivation to attempt to introduce CHA directly into the GB model.

Table 1

Root Mean Square Error (rmse, in kcal/mol) of the GB ΔGpol Relative to the Explicit Solvent (TIP3P) Referencea

radii set	small molecules	amino acid analogs
Bondi[48] (10)	1.55	1.99
Parse[49] (10)	2.33	7.45
ZAP9[50] (13)	0.82	2.88

The numerical R6 GB calculations[47] (see Methods) are based on three common sets of atomic radii;[48−50] the number of atom types for each radius set is shown in parentheses. The corresponding errors in numerical PE ΔGpol (not shown) are not smaller.

Elsewhere,[56] we presented a derivation of a CHA-aware analogue of the Born equation,[57] which can be regarded as a single spherical charge limit of the generalized Born approximation. A key novelty of our approach lies in a distinct separation between the placement of the solute dielectric boundary and the manner in which one accounts for the asymmetric response of the solvent (water) to the sign of the solute charge. Within the approach, the CHA effects no longer have to be mimicked by a redefinition of the dielectric boundary via ad hoc adjustments to the atomic radii. Instead, the asymmetric response is accounted for by a scaling factor applied to the Born electrostatic solvation energy; the scaling contains the details of the asymmetric response specific to the water model one wants to emulate. Here, we apply the strategy to derive a CHA-aware analogue of the GB equation, which we call the charge hydration asymmetric–generalized Born (CHA–GB) model. By direct comparison with explicit solvation energies, we will show how CHA–GB improves the accuracy of the GB framework. The new model was also compared with several other recent semi-continuum models such as SEA[51] and 3D-RISM.[52−55]

Charge Hydration Asymmetric–Generalized Born: CHA–GB

The “CHA-aware” modified Born formula for a single spherical ion of radius ρ and charge q is given by[56]where εin and εout are dielectric constants of the solute and the solvent, respectively, and ρw = 1.4 Å is the water probe radius. The first part in eq 1 describes the charge-symmetric dielectric response of the original Born equation,[57] accurate in the limit of a continuum “charge-symmetric” solvent. (Throughout this article we will loosely use the word “symmetric” for models devoid of charge hydration asymmetry.) One such model is the mean spherical approximation[58] (MSA) for which the Born model is exact; the dielectric boundary is shifted up by Rs = 0.52 Å relative to the surface of the spherical solute of radius ρ. Later in this work, we will extend this dielectric boundary definition to solutes. Once the charge-symmetric solvation is accounted for in the “Symmetric Born” part of eq 1, the charge hydration asymmetry (CHA) correction enters via a multiplicative scaling factor η. In η, the CHA of a water model being emulated is controlled by the model dependent parameter ROH = Q̃/p, where Q̃ and p are respectively the primitive (non-traceless) quadrupole and dipole moments of the n-point water model (Table 2).

Table 2

ROH = Q̃/p for the Three Water Models Used in This Worka

Q̃ = ∑qz2 and p = qz, where z is the azimuthal-symmetry coordinate of charge q with respect to the molecule center. A geometric interpretation of ROH for the TIP3P water model is shown in the schematic on the right.

We now seek an analogue of the CHA-aware Born equation, eq 1, for molecules. Without the CHA, such analogue is well known—the canonical (charge-symmetric) generalized Born (GB) equationwhere q is the atomic partial charge of atom i(j). The empirical function fGB iswhere r is the distance between the ith and jth atom, and R and R are their effective Born radii. The effective Born radius of an atom characterizes its degree of burial inside the solute. Many efficient methods exist for computing this quantity (see ref (33) for a recent review).

The generalized Born approximation can be viewed as an interpolation between the two limiting cases: case (a) r → ∞, that is, charged atoms are treated as separate distant (non-interacting) charged spheres and case (b) r → 0, that is, all charges merge to form a single sphere of charge ∑q and effective Born radius R. The empirical correction, exp(−r2/(4RR)) in fGB, eq 3, originally introduced by Still et al.,[59] performs the interpolation in a way that at least partially accounts for non-spherical shapes of realistic molecules.[4] Our strategy is to utilize the same interpolation to derive the CHA-aware analogue of eq 2 by matching the asymptotes, (a) and (b), to the CHA-aware Born formula (eq 1). In what follows, we consider these two limiting cases in detail.

Case (a): From eq 2, the ΔGpol = (−1/2)(1/εin – 1/εout)∑(q2/R), that is, the sum of self-energies of every atom in the molecule. Comparing with the symmetric Born part of eq 1, the effective Born radius, R, is seen as equivalent to the dielectric boundary radius, ρ + Rs. To add CHA correction in this limit, we scale each GB self-energy term by the corresponding η = {1+sgn[q]ROH/(R – Rs + ρw)}−1 from eq 1. This implies ΔGpol → (−1/2)(1/εin – 1/εout)∑(q2/R)η as r → ∞.

Case (b): In this limit, all the charges are merged together to form a single sphere of effective Born radius, R, with the total charge of ∑q. Using charge-symmetric eq 2, we can write ΔGpol = (−1/2)(1/εin – 1/εout)(∑q)2/R. To add the CHA according to eq 1, the overall ΔGpol must now be scaled with the CHA factor that depends upon the net charge of the molecule, that is, η = {1 + sgn[Σq]ROH/(R – Rs + ρw)}−1. This implies ΔGpol → (−1/2)(1/εin – 1/εout){∑(q2)/R}η as r → 0.

In short, for case (a), each self-term of eq 2 has to be scaled by its respective CHA scaling factor, whereas for case (b), a global CHA-scaling is needed. To consistently satisfy both of these asymptotes while preserving the overall GB functional form set by eqs 2 and 3, we need a single functional form for the CHA scaling factor η that is aware of other charges in the molecule. To this end, we make η a function of ∑qe–τ and scale the effective Born radii accordinglywhere τ is a positive constant that controls the effective range of the neighboring charges (j) affecting the CHA of atom (i). For realistic molecules, the r’s are of the same order of R’s, and hence, it can be argued that τ is of the order of 1. Indeed, the value of τ that we obtain by systematic fitting against explicit solvent reference ΔGpol’s (see Methods) is close to 1. Using eq 4, the charge hydration asymmetric–generalized Born (CHA–GB) approximation can be written aswhere fCHA–GB = (r2 + R̃R̃ exp (−r2/(4RR)))1/2. In line with our philosophy of separating the CHA effects from the purely geometric dielectric boundary, we do not use the CHA-corrected effective radii in the exp (−r2/(4RR)) factor due to its “geometry correction” role in the GB model.

Charge-Symmetric Dielectric Boundary

The dielectric boundary, which separates the low dielectric solute from the high dielectric solvent, is a key concept in the standard continuum solvation approach. Just like in the original GB model, the boundary is needed in our CHA–GB formalism for the computation of the effective Born radii via an integral over the solute surface (see Methods). The dielectric boundary definition we propose, Figure 1, is a generalization of the single spherical ion case, eq 1, designed to account for the charge-symmetric part of the solvent dielectric response. The idea is to make the boundary manifestly invariant to inversion of charges within the solute, as is the case in the classical macroscopic continuum electrostatics framework[60] on which the GB or PE are based. Charge hydration asymmetry of a more realistic solvent has to be added as an external correction to electrostatic solvation energy, for example, as done above for the case of the proposed CHA–GB formalism. A sharp dielectric boundary is defined by a set of atomic radii and the rolling probe radius. In principle, one could find a set of optimal charge-symmetric atomic radii ρ by an optimization against solvation energies computed within a charge-symmetric explicit water model such as BNS[61] or an even simpler model in which the appropriate point dipole is placed at a single vdW center.[56] In practice, we will first add CHA via eqs 4 and 5, and then iteratively optimize the radii to obtain best fit against a more realistic charge-asymmetric explicit solvent (TIP3P) ΔGpol for a representative set of realistic molecules.

Figure 1

Charge-symmetric dielectric boundary used in CHA–GB for multi-atom solutes. The schematic inset shows the key construction steps. Each atomic radius (ρ, red circle) is increased (purple area) by the same correction Rs = 0.52 Å. The dielectric surface (blue line) is outlined by rolling a probe of radius ρw – Rs. The smaller than the standard ρw = 1.4 Å probe ensures invariance of the solvent accessible surface around the solute.

Methods

Molecule Sets

Four molecule sets[38,62−64] were used to test the performance of the proposed model.

Set 1: This is a polar solute series[38] of several “bracelet”-like regular planar polygonal (triangle to octagon) molecules made with a single aromatic carbon “ca” atom type (Figure 2). These charge-inverted bracelets were specifically designed to study solvent-dependent CHA effects.[38] Two charge configurations for each molecular geometry were considered, the N and P bracelets, that is, one atom in the molecule charged −1e (N bracelet) or +1e (P bracelet), whereas the rest of the atoms were equally charged such that the overall molecule is neutral. The explicit ΔGpol in TIP3P and TIP4P-Ew were obtained using standard thermodynamic integration (TI)[65] in Amber 12,[66] whereas the TIP5P-E ΔGpol were deduced from ref (38) (Supporting Information).

Figure 2

Polar solvation free energies of the charged-inverted “bracelets”. Atomic structures of the N-bracelets are shown on the top horizontal axis, while the P-bracelets are shown at the bottom. Explicit solvent energies for three different water models ((a)TIP4P-Ew, (b)TIP3P, and (c)TIP5P-E) are denoted by red dots for N-bracelets and blue dots for P-bracelets. The corresponding CHA–GB energies are shown by the black triangles; the ROH parameter that controls the propensity for CHA in CHA–GB, eq 5, is set to the value appropriate for the given explicit water model (Table 2).

Set 2: A set of 248 rigid diverse neutral small molecules was selected from a larger set of 504 molecules commonly used to test solvent models;[47,51,62,67] the corresponding explicit water (TIP3P) ΔGpol values were taken from ref (62). The rationale for focusing on the rigid molecules is as follows. It is well known that solvation free energies are sensitive to the solute conformations;[67] poor conformational sampling may lead to errors in estimated solvation energies relative to the experiment. However, our primary goal here is to improve the implicit solvent estimates of the polar part of solvation energy, and so we believe that a comparison with the single conformation polar solvation free energy in the explicit solvent is the best strategy. That is, we would like to minimize the number of other factors that affect the ΔGpol estimates. In order to compare with the explicit ΔGpol, we therefore chose the most rigid 248 molecules that showed small conformational variability during the molecular dynamics simulations performed on the same 504 molecule set as in ref (67) (see the Supporting Information for details).

Set 3: This is a set of 48 structures of single amino acid analogs of the form N-acetyl-X-N′-methylamide (X refers to one of the 20 amino acids). Among these, 40 structures were taken directly from ref (63), to which we added eight additional molecules by neutralizing the charged structures: ASP, LYS, GLU, and ARG (Supporting Information).

Set 4: A set of 19 neutral small protein structures with average sizes of 30 amino acid residues were randomly selected from a larger set used in ref (64) (Supporting Information).

Parameter Optimization: Realistic Molecule Sets

For CHA–GB, 10 parameters (ρ of nine atom types, Table 3, and τ = 1.47) were optimized against the explicit (TIP3P) solvation free energy with a training set comprised of 124 molecules from the rigid small molecule set and 24 molecules from the set of amino acid analogs. These molecules in the training set were sampled to approximate an equal representation of the span of solvation free energies as well as of the atom types of the entire set (Supporting Information). We used the Nelder–Mead simplex algorithm[68] with the objective function that weighs equally RMS deviations from the reference ΔGpol of the small molecules and the amino acid analogs. Convergence to the optimum set of parameters (Table 3) was ensured within 100 independent runs with random initial guess parameters (within a physical range) (Supporting Information). For all realistic molecules, we use the analytical linearized Poisson–Boltzmann (ALPB) model,[69] which is a correction to GB, eq 2, that introduces a physically correct dependence on the dielectric constants while keeping the computational efficiency of the original. In ALPB, ΔGpol is approximated aswhere β = εin/εout, α = 0.571412, and A is the electrostatic size of the molecule, essentially the overall size of the structure computed analytically.[69] We add the CHA correction to ALPB by replacing fGB by fCHA–GB similar to eq 5. To avoid uncertainties associated with the effective Born radii estimation, we use numerical “R6” radii,[70] which were shown[47,71] to deliver ΔGpol estimates essentially as accurate as those based on the perfect[72] PE-based radii. Specifically, the effective radii are computed viawhere r is the position of the infinitesimal surface element dS that encloses the charge-symmetric dielectric boundary ∂V (Figure 1), and r is the coordinate of atom i. Unless otherwise stated, for canonical GB, we use the Lee–Richards molecular surface[73] obtained with the standard water probe radius ρw = 1.4 Å. We use a numerical approximation of eq 7, representing the dielectric boundary as a triangulated mesh with density of 6 vertices/Å2 obtained using MSMS[74] package. For proteins structures, following ref (71), we add a uniform correction of 0.028 Å–1 to the R–1 in eq 7.

Table 3

Optimized Intrinsic Atomic Radii Sets for CHA–GB and GB

	radii set (Å)
	C	H	N	O	S	F	Cl	Br	I
CHA–GB	1.56	0.47	1.59	1.37	1.88	1.44	1.84	1.92	2.29
GB	1.76	1.29	1.46	1.50	2.04	1.16	1.25	2.04	1.72

Parameter Optimization: N and P Bracelets

The intrinsic radius of the aromatic carbon, “ca” atom type, was obtained by fitting against the TIP4P-Ew ΔGpol – optimal ρ(“ca”) = 1.37 Å. The same value of τ = 1.47 as above was used in eq 4.

3D-RISM

We computed single point 3D-RISM ΔGpol (TIP3P) using the implementation[75] in AMBER-12[66] with the protocol from http://dansindhikara.com/Tutorials/. These estimates were adjusted[76] by two parameters, a1 and a2, fit against the explicit solvent ΔGpolwhere ΔGpol3D-RISM/GF is the 3D-RISM ΔGpol using the Kovalenko–Hirata closure[77] assuming Gaussian fluctuation of the solvent, ΔGpolcorr is the corrected polar solvation energy, V is the partial molar volume, and ρ = 0.0333A–3 is the solvent number density. The optimized values of a1 and a2 are provided in the Supporting Information.

Performance of the New Model

CHA–GB closely reproduces the explicit polar solvation energy (ΔGpol) (Figure 2) of the charge inverted N and P “bracelet” molecules in TIP3P (rmse = 0.67 kcal/mol), TIP4P-Ew (rmse = 1.02 kcal/mol), and TIP5P-E (rmse = 0.50 kcal/mol) water models emulated in CHA–GB by using the appropriate value of ROH from Table 2. Note that using identical (for all water models) dielectric boundaries for the N and P bracelets of same size properly captures the propensity for CHA of each of the three explicit water models tested. It is interesting that within CHA–GB, the CHA contributions from both self-terms and cross terms in eq 5 are significant. For example, for TIP3P, the ratio of magnitude of CHA in cross terms to that in the self-terms, ⟨|ΔΔGpol(cross)/ΔΔGpol(self)|⟩ = 0.6. These contributions have opposite effects on the overall CHA; the CHA in self-terms is partially negated by the CHA in cross terms. As expected, the canonical GB does not capture the CHA (Figure 2a); the best fit to explicit solvent ΔGpol results in very different values of ρ(“ca”) for different water models as the readjustment of the dielectric boundary compensates for the missing CHA. The inability of the Poisson treatment to capture CHA of the same “bracelet” molecules was noted earlier.[38] In contrast, our fairly simplistic way of introducing the CHA into the linear response continuum already captures the effect at a level similar to considerably more complex “beyond Poisson” models such as 3D-RISM[76,78] and SEA[51] (Figure 3). A noteworthy observation is that the atomic radii optimized for CHA–GB have physically meaningful variations with the atom type (Table 3) unlike the ones optimal for GB that have ρ(I) < ρ(Br) and ρ(H) > ρ(Cl). It is also reassuring that the optimal hydrogen radius ρ(H) = 0.47 Å in the CHA–GB model is such that the realistic distance[79] of about 2 Å observed between the hydrogen and its electronegative partner in the hydrogen bond can now be reproduced. Compared against the reference (TIP3P) ΔGpol, CHA–GB shows approximately 40% improvement in accuracy (root mean square error, rmse) over the canonical GB (Table 4) for the 248 small molecules and 48 amino acid analogs. The largest errors and the average errors are also reduced in the CHA–GB model. For small proteins, CHA–GB also shows an improvement in accuracy; note that this type of structure was not used to train the radii sets. Arguably, the most impressive achievement of CHA–GB is that it provides equally good accuracy of better than 1 kcal/mol simultaneously for both the small molecules and amino acid analogs. This is in contrast to GB. For example, GB with a ZAP9 set (13 atom types) specifically optimized for small molecules can yield similar good accuracy (ΔGpol rmse = 0.82 kcal/mol) for the 248 small molecules, but then the corresponding errors for the 48 amino acid analogs become large (ΔGpol rmse = 2.88 kcal/mol) (Table 1). We have checked directly that the uniformity of the accuracy improvement or equivalently the improved transferability of the new radii set is indeed the consequence of introducing the CHA into the GB model. Namely, even if the new dielectric boundary of Figure 1 is used in the charge-symmetric GB, the model still cannot reach the 1 kcal/mol rmse accuracy of ΔGpol simultaneously for both the small molecules and amino acid analogs(Supporting Information). Our overall conclusions do not change, and the accuracy of CHA–GB remains essentially the same if we extend the molecule set to include the flexible molecules from the original set of 504 neutral molecules (Supporting Information). A comparison of the performance of CHA–GB and 3D-RISM[77] is presented in Table 5. Note that the same training set and optimization technique (see Methods and Supporting Information) was used to obtain best fit parameters for both models. A systematic comparison of CHA–GB with other “beyond PE” models is out of the scope here. However, we note that the published performance of the SEA[51] model relative to the explicit solvent on the same set of 504 small molecules (same partial atomic charges) (ΔG rmse = 0.81 kcal/mol) is comparable to that of CHA–GB (rmse = 0.89 kcal/mol).

Figure 3

Asymmetry in the polar part of solvation free energies, ΔΔGpol = |ΔGpol(N) – ΔGpol(P)| for N/P bracelets using different methods. For the SEA model, the data are adapted from ref (51).

Table 4

Accuracy of ΔGpol Estimated by CHA–GB and GB Based on Their Respective Optimal Atomic Radii (Table 3)a

	small molecules		amino acid analogs		proteins
	GB	CHA–GB	GB	CHA–GB	GB	CHA–GB
rmse	1.24	0.88	1.26	0.81	10.24	8.95
⟨err⟩	–0.53	–0.37	0.24	0.09	–6.12	–2.80
⟨|err|⟩	0.93	0.63	0.90	0.64	7.79	7.63
r2	0.86	0.93	0.998	0.999	0.99	0.99
%(|err| > 2kBT)	30.6	14.9	25.0	16.7	84.2	89.5
rmse worst 5%	3.20	2.55	4.14	2.11	25.53	16.20

The accuracy is assessed relative to explicit solvent (TIP3P) ΔGpol (kcal/mol): root mean square error (rmse), mean error (⟨err⟩), mean absolute error ⟨|err|⟩, r2 correlation, percentage of molecules with absolute error >2kBT, and rmse of the 5% molecules with largest |err|.

Table 5

Error (rmse relative to explicit solvent, in kcal/mol) in ΔGpol Estimates Using CHA–GB and 3D-RISM

	small molecules	amino acid analogs	proteins
CHA–GB	0.88	0.81	8.95
3D-RISM	0.50	5.28	18.36

Because many disparate approximations separate reality from the linear response continuum framework,[80] it is noteworthy that adding only CHA to it (here via GB) yields accuracy in ΔGpol estimates similar to or even better than those of some more complex approaches that implicitly incorporate not only CHA but also other explicit solvent effects missing from the linear response continuum. At the same time, the reported computational cost of such methods, for example, SEA[51] and especially of 3D-RISM,[75] is much higher than that of the canonical GB,[47] while CHA–GB is only about ∼1.5 folds slower than GB.

By adding the non-polar contribution to the polar part (Supporting Information), we can compare the total solvation energy against the experiment. The solvation free energy predicted by CHA–GB (rmse = 1.22 kcal/mol) against the experiment is about 20% more accurate than that of the GB (rmse = 1.45 kcal/mol) and is comparable to the accuracy of calculations in explicit (TIP3P) water[62] (rmse = 1.26 kcal/mo) for the same 504 small molecules. A standalone code to perform the GB and CHA–GB calculations described here can be downloaded from http://people.cs.vt.edu/~onufriev/software.php.

Conclusion

The work has two main results. First, we have demonstrated how charge hydration asymmetry (CHA)—so far missing from the conceptual framework of the linear response continuum electrostatics—can be incorporated into what is arguably the simplest commonly used approximation based on the framework, namely, the generalized Born model. The resulting CHA–GB model is as conceptually simple and should be nearly as computationally efficient as the original. At the same time, it provides a noticeable improvement in the accuracy of electrostatic solvation energy estimates relative to explicit solvent reference. What is perhaps even more noteworthy is that the introduction of CHA leads to a uniform improvement in the performance of the model on different classes of molecules, including charged compounds. Note that standard continuum electrostatic calculations based on commonly used sets of intrinsic atomic radii do not show the uniform accuracy across a range of charge states and structures, which had led to the development of specialized radii sets. These atomic radii often exhibit unphysical trends in the atom sizes, for example, radius of fluorine may become smaller than that of hydrogen. This is where the second main result of the work is: the development of a new dielectric boundary definition to which the CHA effects are “external”. In the proposed formalism, the dielectric boundary is manifestly charge-symmetric, while the CHA effects are added to the computed solvation energy as a correction. Within this framework, atomic radii optimized for solvation energy calculations against explicit solvent reference show physically meaningful variation with the atom type. As is the case with any model, there are limitations of the proposed approach. The CHA scaling has so far been incorporated by matching the two obvious limiting cases of the original generalized Born equation. While these two limits are physically well grounded, the intermediate region is treated simply as a mathematical interpolation. Another limitation of our current model is that the CHA is solely determined by the sign of the solute charge, that is, independent of its magnitude. While this is a good approximation for ions (integer charge), a more sophisticated approach may be needed to treat fractional atomic charges.[41] One should also not forget that even at the level of continuum solvent various CHA-unrelated microscopic solvent effects[81] are not considered by the proposed CHA–GB model. Nevertheless, it is noteworthy that the accuracy of the basic continuum electrostatics model can be improved noticeably by adding only CHA to it.