Extending the Solvation-Layer Interface Condition Continum

Article pubs.acs.org/JCTC

Extending the Solvation-Layer Interface Condition Continum Electrostatic Model to a Linearized Poisson−Boltzmann Solvent Amirhossein Molavi Tabrizi,† Spencer Goossens,† Ali Mehdizadeh Rahimi,† Christopher D. Cooper,‡ Matthew G. Knepley,¶ and Jaydeep P. Bardhan*,† †

Department of Mechanical and Industrial Engineering, Northeastern University, Boston, Massachusetts 02115, United States Departamento de Ingeniería Mecánica and Centro Científico Tecnológico de Valparaíso (CCTVal), Universidad Técnica Federico Santa María, Valparaiso, Chile ¶ Department of Computational and Applied Mathematics, Rice University, Houston, Texas 77005, United States ‡

ABSTRACT: We extend the linearized Poisson−Boltzmann (LPB) continuum electrostatic model for molecular solvation to address charge-hydration asymmetry. Our new solvationlayer interface condition (SLIC)/LPB corrects for first-shell response by perturbing the traditional continuum-theory interface conditions at the protein−solvent and the Sternlayer interfaces. We also present a GPU-accelerated treecode implementation capable of simulating large proteins, and our results demonstrate that the new model exhibits significant accuracy improvements over traditional LPB models, while reducing the number of fitting parameters from dozens (atomic radii) to just five parameters, which have physical meanings related to first-shell water behavior at an uncharged interface. In particular, atom radii in the SLIC model are not optimized but uniformly scaled from their Lennard-Jones radii. Compared to explicit-solvent free-energy calculations of individual atoms in small molecules, SLIC/LPB is significantly more accurate than standard parametrizations (RMS error 0.55 kcal/mol for SLIC, compared to RMS error of 3.05 kcal/mol for standard LPB). On parametrizing the electrostatic model with a simple nonpolar component for total molecular solvation free energies, our model predicts octanol/water transfer free energies with an RMS error 1.07 kcal/mol. A more detailed assessment illustrates that standard continuum electrostatic models reproduce total charging free energies via a compensation of significant errors in atomic self-energies; this finding offers a window into improving the accuracy of Generalized-Born theories and other coarse-grained models. Most remarkably, the SLIC model also reproduces positive charging free energies for atoms in hydrophobic groups, whereas standard PB models are unable to generate positive charging free energies regardless of the parametrized radii. The GPU-accelerated solver is freely available online, as is a MATLAB implementation. asymmetry.14−18 Hydration asymmetry is readily seen for the case of a Born ion: an anion and cation of equal valence and radii have different solvation free energies and entropies. However, most continuum models predict solvation free energies that do not depend on the sign of the charge, and fail to predict entropies entirely.19 To account for noncontinuum response of waters in the hydration shell, several groups have developed implicit-solvent theories to account for such asymmetries in solvation free energies.20−25 The Born equation, which relates an ion’s radius and its electrostatic solvation free energy, motivates one of the prevalent approaches to incorporate charge-sign asymmetry: adjusting atom radii to match solvation free energies. For Born ions, radii can always be parametrized to fit exactly, because each radius is fit independently to a single free energy. The necessary adjustments can be rationalized on physical grounds,26 and

1. INTRODUCTION Many implicit-solvent models treat the electrostatic component of a molecular solute’s interactions with the surrounding solvent using macroscopic dielectric theories and the Poisson or Poisson−Boltzmann equation.1−3 As a result of the significant approximations required to apply semimacroscopic models at atomistic or near-atomistic length scales,1,4 these simple but fast models can exhibit unpredictable deviations from experimental results5 or reference free-energy calculations from explicitsolvent calculations.6 Many efforts have been made to address the shortcomings of the Poisson−Boltzmann equation,7−9 for example; in addition, several theories have been proposed that use nonlinear dielectric response to model dielectric saturation at high field strength.10,11 Popular Poisson-based continuum theories also cannot treat hydrogen bonds between solvent and solute, motivating the development of a variety of corrections.12,13 This paper addresses a different long-standing deficiency in popular continuum theories, the problem of charge-hydration © XXXX American Chemical Society

Received: August 22, 2016 Published: April 5, 2017 A

DOI: 10.1021/acs.jctc.6b00832 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

to emphasize that our model does not address the longstanding controversies regarding the “correct” choice of a solute−solvent interface, or regarding the density of water near hydrophilic or hydrophobic surface patches (for discussion, see ref 33). In a similar vein, the charge density on the nominal boundary should be regarded as an effective charge density, in the sense that it does not model the local charge density that one might calculate by integrating a thin shell at the solute− solvent interface. The effective charge density is instead intended to reproduce the correct reaction potential f ield inside the solute. In other words, the nonlinear boundary condition leads to a nonlinear boundary-integral equation (BIE), whose solution (an effective surface charge) correctly reproduces the solute reaction potential. In reality, of course, the solvent responds with a volumetric polarization charge density, and this volume charge density exhibits numerous complicated behaviors and oscillations. Inside the solute, however, the resulting reaction-potential must still satisfy the Laplace equation, and for any boundary that contains the solute charges, one can find a surface charge density that reproduces the proper reaction potential.29 We proposed the modified interface condition considering the molecular surface (solvent-excluded surface) to be the nominal dielectric boundary.32 Modified interface conditions, or at least reparameterized versions, will be needed for alternative surface definitions or solvent models. For example, we recently developed a version of SLIC using the mean spherical approximation (MSA) model34,35 as a basis to predict ion solvation thermodynamics in different solvents.36,37 This SLIC/MSA model did not address charge hydration asymmetry, and instead showed that solvation thermodynamics can be modeled accurately using a temperature-dependent SLIC model, that is, a temperature-dependent perturbation to the usual (macroscopic) dielectric interface condition. We note that a wide variety of implicit-solvent models are presently under development:38 the spectrum includes liquidstate integral equations,39−41 free-energy functional and variational methods,42−45 advanced Generalized-Born (GB) models,13,28 Fennell and Dill’s Semi-Explicit Assembly method,46 the morphometric theory of Roth et al.,47,48 local molecular field (LMF) approaches,49−51 inhomogeneous solvation theory.52,53 We also note that Pomogaeva and Chipman have introduced an electrostatic continuum model that adds terms related to the maximum and minimum electric field on the boundary.54 The more sophisticated models frequently treat the nonpolar and electrostatic components of solvation in a unified way, and the charge models themselves have a large impact on computed free energies.55 The simple model that we propose here reproduces the electrostatic component of the solvation free energy as defined by Roux and Simonson,1 who expressed the solvation free energy as a sum of two terms: first, the free energy associated with creating the uncharged molecular cavity in the solvent, and second the free energy of creating the solute’s charge distribution in the cavity. The latter quantity is also the charging free energy, if one performs MD free-energy calculations; however, it should be noted that this free energy includes changes in nonpolar interactions that occur as part of the charging process. Remsing and Weeks have used LMF recently to assess different thermodynamic cycles associated with these steps, which offers complementary views into the roles of steric effects and the charging process.49 The following section introduces the standard, macroscopic dielectric model as well as boundary-integral approaches to

Latimer et al. presented an early, compelling argument based on solvent sterics, namely, that positively charged water hydrogens can approach a negative ion more closely than can the negatively charged water oxygens approach a positive one.14 On the other hand, many studies, including our own, have demonstrated that buried atoms can experience significant asymmetric solvation even if adjusting their radii would not affect the dielectric boundary.27−29 A more detailed, or at least nuanced, model is therefore needed. The solute−solvent interface potential field16,29,30 also contributes to asymmetries in hydration free energies. Using large-scale molecular dynamics (MD) simulations of a fully uncharged protein surrounded by explicit solvent, Cerutti, Baker, and McCammon illustrated that this field exhibits complicated behavior near the solute−solvent interface, but is nearly constant in the protein interior.31 This field, which arises from a water structure at and near the interface, is fundamentally noncontinuum in nature: modeling an uncharged protein using standard continuum theory would result in the electrostatic potential being identically zero. To emphasize that the nonzero field exists independently of the solute charge distribution, we termed this the static potential f ield φstatic. By investigating hydration asymmetry (i.e., ΔGcharging(q = +1) − ΔGcharging(q = −1)) as a function of a charge’s distance from the center of a sphere of radius 5 Å, we found that the static potential dominates asymmetry for buried charges, whereas for solvent-exposed charges the solvent-steric effect contributes slightly more than the static potential.29 The rapid decay of the solvent-steric effect led us to question whether it would be possible to formulate a charge-sign-asymmetric model in which the dielectric boundary did not depend on the solute charge. Analysis of the Latimer, Mukhopadhyay, and Purisima asymmetries led us to propose a solvation-layer interface condition (SLIC) that captures the effect of changing the boundary.32 The present work extends the earlier SLIC model to account for the static potential as well as the solvent-steric component of asymmetry. We also extend SLIC to address ionic solutions via the linear Poisson−Boltzmann equation. We previously demonstrated that in nonionic solutions, SLIC is more accurate than standard Poisson models for calculating the charging free energies of titratable amino acids.32 These results suggested that SLIC might be able, in principle, to provide more accurate calculations of protein pKa values and pKa shifts; however, biologically relevant calculations must also model the effects of mobile salt ions, which our earlier work did not address. Our model therefore needed to be extended at least to the linearized Poisson−Boltzmann (LPB) equation. This extension is the subject of the present work. The new model incorporates both of the dominant mechanisms underlying charge-hydration asymmetry: the static potential via a constant-field approximation, and the steric-size effect via the SLIC. Calculating the electrostatic potential in the solvent region required a modification to SLIC, to account for the fact that in the solvent, the reaction field in the original SLIC model does not identically satisfy Gauss’s law (being an effective charge approach, this was not especially surprising). The new model satisfies Gauss’s law by construction, and results show a significant improvement over the original SLIC model, as well as comparable salt dependence. Throughout this paper we will refer to the (chargeindependent) solute−solvent interface as the nominal dielectric boundary, or the effective boundary. This language is intended B


Article

Journal of Chemical Theory and Computation solving it. We also detail the mechanisms underlying asymmetry and the motivation for the solvation-layer interface condition (SLIC). Section 3 presents the new model, which includes the static-potential contribution to asymmetry as well as a second charge layer at the Stern layer, to ensure that outside the solvation layer, the potential obeys the usual macroscopic Gauss’s law. This section also presents a boundary-integral formulation for the new model, and details on numerical implementation. Calculations on Born ions and amino acids are presented in section 4, as are comparisons to detailed MD simulations and a brief investigation of the model’s implications for protein titration. Section 5 discusses some of the simplifying assumptions and modeling approaches in light of the results. Section 6 concludes the paper by summarizing the work and highlighting questions for further research.

where κ is the inverse Debye screening length. In the solute, region III, the solute’s charge distribution ρ(r) consists of Nc discrete point charges, the ith of which is located at ri and has value qi. In region III, the potential satisfies the Poisson equation

2. THEORY 2.1. Continuum-Electrostatic Models with Linear Boundary Conditions. Figure 1 is a schematic of the problem

The surface a models the dielectric boundary separating the solute from the ion-exclusion layer. The surface normal is denoted as n̂(r) where the hat denotes unit length, and defined to point outward into region II. Similarly, the Stern surface b separates the ion-exclusion layer and the electrolyte, and its normal points outward into region I. Three boundary conditions complete the specification of the problem. First, the potential is assumed to decay sufficiently quickly as |r| → ∞. Second, the potential is continuous across the interfaces:

∇2 φIII(r) = −

(2)

Between the electrolyte solution and the solute is the Stern layer, region II, which is a high-dielectric region without mobile ions. In this region the potential satisfies the Laplace equation

∇2 φII(r) = 0

(3)

δ→ 0

lim[φII(ra + δn(̂ ra)) − φIII(ra − δn(̂ ra))] = 0

(4)

lim[φI(rb + δn(̂ rb)) − φII(rb − δn(̂ rb))] = 0

(5)

δ→ 0

Figure 1. Schematic of a continuum-electrostatic model for a single solute (region labeled III) in an infinite solvent (region labeled I), where the solvent is modeled as a dilute electrolyte obeying the linear Poisson−Boltzmann equation. In the ion-exclusion layer (labeled II), which is bounded by the dielectric boundary a and the Stern surface b, the potential obeys a Laplace equation with the solvent dielectric constant. In the new SLIC model, the solvation-layer interface condition is enforced at the dielectric boundary a and the Gauss-Law correction is enforced at the Stern boundary b, regardless of whether the potential in region I is modeled using the Laplace equation or the linear Poisson−Boltzmann equation.

The present work focuses on the third boundary condition, which relates the normal components of the electric displacement fields. For a linear dielectric medium, the displacement field is D(r) = ϵ(r)E(r)

(6)

= −ϵ(r)∇φ(r)

(7)

and in standard continuum-electrostatic theory, the normal displacement field is continuous across the boundary; that is,

we consider in this paper. The dilute electrolyte solvent extends to infinity and is labeled region I. The thin solvent region immediately surrounding the solute, called the ion-exclusion or Stern layer, is inaccessible to dissolved ions and labeled region II. The region corresponding to the interior of the solute is labeled region III, and the only fixed charges in the system are located in III. Ion-exclusion layers are usually defined to have a thickness corresponding to the radius of the mobile ions in the solution, and its overall effect is to limit ionic screening because the actual charges of finite-size ions cannot reach the actual protein surface.56 In the present work, all Stern layers are 2 Å thick. The dielectric constant in each region is denoted by the symbol ϵ with the appropriate subscript; for example, ϵIII is the solute dielectric constant, and usually ϵII = ϵI. In the solvent (region I) the potential satisfies the linearized Poisson−Boltzmann equation ∇2 φI(r) = κ 2φI(r)

ρ(r) ϵIII

lim n(̂ ra)· (DII(ra + δn(̂ ra)) − DIII(ra − δn(̂ ra))) = 0

δ→ 0

(8)

so ϵI

∂φI

ϵII

∂n

(rb) = ϵII

∂φII ∂n

∂φII ∂n

(ra) = ϵIII

(rb)

∂φIII ∂n

(ra)

(9)

(10)

Using Green’s theorem and well-known derivations,57,58 we may obtain a coupled system of four boundary-integral equations (BIEs) for the potential and normal derivative on the dielectric boundary and the Stern surface:

(1) C


Article

Journal of Chemical Theory and Computation ⎡⎛ 1 ⎞ ⎢⎝⎜ I + KIII,aa⎠⎟ − G III,aa 2 ⎢ ⎢⎛1 ϵ ⎞ ⎢ ⎜ I − KII,aa⎟ + G II,aa II ⎠ ϵ III ⎢⎝2 ⎢ ϵ ⎢ −K + G II,ba II II,ba ⎢ ϵ III ⎢ ⎢ ⎢ ⎣

+ KII,ab ⎛1 ⎞ I + KII,bb⎟ ⎝2 ⎠

⎜

⎛1 ⎞ I − KI,bb⎟ ⎝2 ⎠

⎜

Solvent finite-size effects, which represent the second major contribution to charge-hydration asymmetry, can be seen in the second derivative of the charging free energy curve around q = 0. The second derivative changes from one value to another, in a very narrow window around q = 0, which leads to the two different expressions in eq 12. For solvent-exposed charges, the difference between L(+) and L(−) leads to substantial energetic differences: for a unit charge with sodium’s Lennard-Jones parameters, the difference is 60 kcal/mol; for a unit charge buried in a sphere of radius 5 Å, the difference is still 20 kcal/ mol.29 As understood for decades,14 the asymmetry for solventaccessible charges arises in large part because a positive charge is predominantly solvated by larger water oxygens, which cannot approach the charge as closely as water hydrogens can approach a negative charge. However, the difference becomes relatively small for charges more than approximately 3 Å from the interface,29 motivating the functional form of other asymmetric hydration models.28 In previous work60 we showed that eq 12 successfully models charge-hydration asymmetry for single ions as well as for charges in a larger spherical solute by combining the distinct physical phenomena discussed above. That is, for a single point charge, solvent finite-size effects can be modeled with the signdependent quadratic term, and the linear term associated with φstatic results from water structure independent of the solute charge distribution. The latter dominates for buried charges, and the former is increasingly important as a function of proximity to the solute−solvent interface. Models for chargehydration asymmetry that do not distinguish the two mechanisms cannot be expected to properly resolve the length scales over which solvent-sterics affect asymmetric solvation. In the next section, we discuss modifying the usual macroscopic dielectric interface condition so that the resulting model exhibits the desired sign-dependent quadratic for a single charge. 2.3. Solvation-Layer Interface Condition (SLIC): A Nonlinear Boundary Condition at the Solute−Solvent Interface. We introduce SLIC by way of considering the standard Poisson model, that is, in nonionic solution and using the familiar macroscopic dielectric interface condition (MDIC) of eq 10. The reaction potential can be computed as the Coulomb potential that arises from the induced surface charge at the dielectric interface σMDIC. This surface charge satisfies the BIE known as the polarizable continuum model (PCM)61−66

⎤ ⎥⎡ φIII(ra) ⎤ ⎥ ⎥⎢ ⎥⎢ ∂φ (ra) ⎥ III ⎥ − G II,ab ⎥⎢ ⎥⎢ ∂n(ra) ⎥ ⎥ ⎥⎢ ⎢ ⎥ − G II,bb ⎥⎢ φII(rb) ⎥ ⎥ ⎢ ⎥ ∂φ (r ) ⎥ ⎥⎢ II b ⎥ ϵI + G I,bb ⎥⎢⎣ ∂n(rb) ⎥⎦ ϵ II ⎦

⎡ qG ⎤ ⎢∑ i III ⎥ ⎥ ⎢ 0 ⎥ =⎢ ⎥ ⎢ 0 ⎥ ⎢ ⎦⎥ ⎣⎢ 0

(11)

In eq 11 GII,ab denotes the single-layer potential operator associated with the free-space Green’s function for region II, which maps from a distribution on the boundary b to the potential on the boundary a. Similarly, KII,ab denotes the corresponding double-layer potential operator, φIII(ra) denotes the potential on a (the limit taken from the interior, that is, region III). Full details of this derivation are provided in the Appendix. 2.2. Dominant Physical Mechanisms Underlying Hydration Asymmetry. We previously identified the relative roles of two physical mechanisms which dominate chargehydration asymmetry:29 the static potential field in the solute,16,30 and the influence of the solvent molecules’ finite size and internal charge distribution.14,18,23,28 In this section, we describe the mathematical structure of our near-linear-response continuum model, which we call a piecewise-affine response model, that includes both of these mechanisms. Explicit-solvent MD free-energy perturbation (FEP) calculations29 of monatomic ions and of single charges in a spherical solute showed that for −1e ≤ q ≤ + 1e, charging free energies were well captured by a model of the form

ΔGcharging

⎧ 1 (−) 2 L q + φstaticq , q ≤ 0 ⎪ ⎪2 =⎨ ⎪ 1 L(+)q2 + φ q , q ≥ 0 ⎪ static ⎩2

(12)

Here φstatic is the static potential field at the charge location when q = 0, and L(−) and L(+) represent the proportionality coefficients for solvent response given a negative or positive charge, respectively. The term affine refers to the fact that the electrostatic potential is not a pure linear function of the charge, but has a constant offset, which leads to the linear term in eq 12. The static potential field is the first major component of charge-hydration asymmetry. This field is related to the liquid− vapor interface potential,30,31,59 and its existence is evident from charging-free energy calculations demonstrating that the derivative of the free energy with respect to the charge q is not zero around q = 0 (as predicted by linear-response theory). The slope at q = 0 is in fact the potential induced by the complex solvent structure around even a fully uncharged solute. Cerutti et al. have characterized this potential for a complex solute, a protein.31 They found that this field was nearly constant, which is expected because the field must satisfy the Laplace equation, because it arises exclusively from solvent charges, which are outside the solute.32

⎛ ϵ̂ ⎞ σ MDIC 1 ⎜ I + K ′⎟ = ⎝2 ⎠ ϵIII ϵIII

NC

⎛ ∂G ⎞ ⎟q ∂n ⎠ i

∑ ⎝⎜− i=1

(13)

where ϵ̂ =

ϵII + ϵIII ϵIII − ϵII

(14)

For sufficiently smooth surfaces, the operator K′ always has an eigenvalue −1/2 whose corresponding eigenfunction is a constant.67,68 Now we consider the spherically symmetric case of a Born ion, a point charge q centrally located in a sphere of radius Rnom and dielectric constant ϵIII. We then have a constant surfacecharge density σ on the sphere’s surface, and the electrostatic solvation free energy can be written in terms of the electrostatic interaction energy between the point charge and surface charge distribution: D


Article

Journal of Chemical Theory and Computation 1 1 qR nomσ 2 ϵIII

ΔG =

ratios clearly depend on the local electric field, with a primary dependence on the field’s sign, with a much smaller and more complicated dependence on its magnitude. The field-signdependent ratio observed in the figure motivated capturing asymmetric response using a step function or an approximation to one, such as a function based on a hyperbolic tangent such as

(15)

In the classical dielectric model, eq 13 simplifies because the 1 surface charge distribution is constant, so 2 I + K ′ σ = 0 and

(

)

thus ⎛ ϵIII ⎞ MDIC ∂φ− = ⎟σ ⎜1 + ϵII − ϵIII ⎠ ∂n ⎝

h(En) = α tanh(βEn − γ ) + μ

(17)

In eq 17, α, β, γ, and μ are model parameters with physical interpretations: α models the difference in response between inward and outward pointing fields; β captures the change in electric field strength required to transition between the two response regimes; and γ and μ center the response function. In particular, γ/β is the electric field at which the response is halfway between the two limits, and μ = α tanh(−γ) ensures that the response function matches bulk results in a suitable limiting case of weak electric fields. Said differently, this definition of μ enforces that h(0) = 0, so that in the limit of a vanishing electric field, the model reverts to standard (chargesign symmetric) response. Optimizing the remaining parameters to match the ratios from explicit-solvent MD illustrates that eq 17 seems capable of capturing much of the charge-signdependent response. Here we are using MD free-energy calculations essentially to “computationally measure” ΔG, computing two charging free energies for each of the CHARMM monovalent and divalent ions: one simulation charges the ion from q = 0 to +1e, and the other charges the ion from q = 0 to −1e. That is, to provide a suitably thorough coverage of the relationship between Coulomb field and effective surface charge, we are defining atoms with realistic sizes though not always realizable charges. The narrow transition and near-step-like behavior suggested modifying the PCM integral equation to solvation-layer steric effects using the modified, nonlinear BIE

(16)

−

where ∂φ is the potential just inside the boundary. One ∂n recovers the usual Born ion free-energy by substituting σMDIC prescribed by eq 16. Furthermore, it can be seen that in the standard PB model, the relationship between the ion charge’s ∂φ direct Coulomb field E = ∂n and σ/E is constant. On the other hand, eq 15 makes clear that ΔG and Rnom, when taken together as reference quantities, can be used to establish an ef fective surface charge density σ (on the boundary Rnom) that by definition reproduces the correct ΔG. Using explicit-solvent FEP calculations in the CHARMM force field for charging monatomic ions, and the standard CHARMM radii (Rmin/2) as continuum model radii, we find that the MDderived surface charges follow a surprisingly simple relationship to the electric field (Figure 2). In contrast to the standard PB model where the ratio is a constant (red line), the MD-derived

⎛ ϵ̂ ⎞ σ 1 ⎜ I + K ′ + h(E )⎟ = n ⎝2 ⎠ ϵIII ϵIII

⎛ ∂G ⎞ ⎟q ∂n ⎠ i

∑ ⎜⎝− i

(18)

where the field-dependent operator h(En) is defined so that h(En)σ = h(En(r))σ(r)

(19)

that is, the nonlinear operator represents a pointwise multiplication (local scaling), and En(r) = EnCoul(r) + Enreac(r −)

(20)

N

q ⎛ ∂G ⎞ =∑ qi⎜ − ⎟ + ( −K ′)σ ⎝ ∂n ⎠

Figure 2. Motivating the functional form of the SLIC model’s modified dielectric interface condition, using explicit-solvent MD simulations as a reference and a comparison to the standard macroscopic dielectric interface condition. In this plot, the horizontal axis represents the bare Coulombic electric field at the dielectric boundaries of monovalent Born ions (including ions with opposite sign to their usual). The vertical axis is the ratio of this electric field to the “exact” surface charge that would give the ion’s electrostatic solvation free energy (as calculated using explicit-solvent free energy perturbation32). In the standard macroscopic dielectric model (red line), the ratio is independent of the electric field. Explicit-solvent MD simulations (blue circles) indicate a strong dependence on whether the field points in or out. The proposed functional form for the modified dielectric interface condition can reproduce the correct behavior (black squares). Note that for illustration purposes, here the x axis is the normal derivative of the Coulomb potential due to the atomic charge, not the total electric field as is actually used in the SLIC model; hence, the parameters used to generate this plot are for illustration purposes only (α = 0.3, β = −60, γ = 0.37, and μ = −0.107).

i=1

(21)

Note that the normal electric-field operator is −K′ because the electric field is the negative gradient of the potential; also, the r− in eq 20 is a reminder that we are considering the limit of the electric field as we approach the dielectric boundary from the inside. The narrow transition region suggests that solvation energies should be relatively insensitive to the precise functional form of the transition between the limits. Nevertheless, the specific functional form warrants further investigation because some choices may offer superior theoretical properties, for example, allowing derivation of variational principles, or existence and uniqueness proofs. Earlier work established that the SLIC Poisson model of eq 18 works well for modeling charge-sign hydration asymmetry.32 However, some types of analysis, and numerical calculations E


Article

Journal of Chemical Theory and Computation based on volumetric techniques (finite-difference and finiteelement methods69−72) require an explicit form for the new nonlinear boundary condition in terms of the displacement fields on both sides of the interface. Such an explicit representation is most easily obtained by deriving the standard PCM eq 13 and then introducing the nonlinear term h(En). Equation 13 can be derived by considering an equivalent problem to the mixed-dielectric Poisson problem, in which ρ(r) and the boundary a remain as in the original problem, but the dielectric constant is ϵIII everywhere (making the system a homogeneous dielectric) and we introduce a surface charge σ(r) on a, which is as yet unspecified. The total potential in the homogeneous medium, which we denote by φ̃ (r), is the sum of the Coulomb potentials induced by the solute charge distribution ρ(r) and by the surface charge distribution σ(r): Nc

φ̃(r) =

1 [∑ q G(r; ri) + ϵIII i = 1 i

∫a σ(r′)G(r; r′) dr′]

(22)

f (En) =

ϵIII − h(En) ϵII − ϵIII

(30)

1 ∂φ̃− σ SLIC = 1 + f ∂n ϵIII

(31)

∂φ̃− ∂φ̃+ 1 ∂φ̃− − = ∂n ∂n 1 + f ∂n

(32)

⎛ ∂φ̃+ 1 ⎞ ∂φ̃− = ⎜1 − ⎟ 1 + f ⎠ ∂n ∂n ⎝

(33)

(24)

(1 + f )

(25)

∂φ̃ ∂n

∂φ̃+ ∂φ̃− =f ∂n ∂n

(34)

Therefore, if the perturbation h(En) = 0, we recover the standard macroscopic dielectric interface condition. 2.4. Field Renormalization for Exterior Problems. As demonstrated in previous SLIC studies,32,36 the modified interface condition improves accuracy for the potential in the solute interior by computing an effective surface charge at the dielectric interface. However, just as an image point charge for a planar half-space does not produce the correct potential in the half-space where it resides,73 one cannot assume that SLIC correctly predicts the potential outside the solute. One clear failing is that outside the solute−solvent dielectric boundary, the total displacement field does not satisfy the macroscopic dielectric form of Gauss’s law, which says that for a solute with total fixed charge qTOTAL = ∑qi, one expects that if a surface Γ properly contains the solute region,

Note that for a given σ(r), eq 22 represents a complete solution for the homogeneous-medium problem, giving rise to a particular interior potential. Equation 25 specifies the discontinuity in the electric field across the interface. The standard PCM solves the original mixed-dielectric problem, that is, ensures that the interior potentials φ(r−) and φ̃ (r−) are equal, by specifying that σ(r) must reproduce the displacement-field condition of the original problem (eq 10). By identifying φ̃ (r−) with φIII(r) and φ̃ (r+) with φII(r), we can use eq 7, the simple scalar relationship between the local electric field and the local displacement field, to eliminate

(29)

so that h(En) = 0 recovers the standard macroscopic dielectric interface condition (MDIC). Substituting eq 29 into eq 25 gives the explicit form of the nonlinear displacement-field boundary condition, which is the desired replacement for the MDIC (eq 10):

(23)

+

(28)

where

Last, Gauss’s law implies that in the homogeneous medium, the normal component of the electric field is discontinuous as r passes through the surface: σ(r) ∂φ̃ − ∂φ̃ + = (r ) − (r ) ϵIII ∂n ∂n

(27)

∂φ̃− σ = ϵIII ∂n

(1 + f (En))

and the scaled identity operator arises from assuming that the boundary is sufficiently smooth (has a tangent plane defined everywhere). In operator notation, eq 23 is written ⎛1 ⎞ ⎤ ∂φ(̃ r) 1 ⎡ ∂G ′ ⎟σ ⎥ = + ⎜ I + KIII,aa ⎢⎣∑ qi ⎝ ⎠ ⎦ 2 ∂n(r) ϵIII ∂n

ϵ ∂φ̃− ∂φ̃− σ − III = ∂n ϵII ∂n ϵIII

Substituting eq 23 into eq 28 gives eq 13. In other words, the boundary conditions define the self-consistency condition that leads directly to the BIE. As we have seen, the Born-ion problem suggests a different, nonlinear BIE (eq 18), and working backward led to the new boundary condition. Noting that the left-hand side of eq 28 expresses the effect of the dielectric constants in a form reminiscent of a perturbation, we proposed a nonlinear analogue of eq 28, namely

N ∂φ(̃ r) 1 ⎡⎢ c ∂G(r; ri) ⎛ 1 = + ⎜ σ(r) + ∑ qi ∂n(r) ϵIII ⎢⎣ i = 1 ∂n(r) ⎝2

∫̵

(26)

⎛ ϵIII ⎞ σ ∂φ̃− = ⎟ ⎜1 + ϵII − ϵIII ⎠ ϵIII ∂n ⎝

To compute the normal derivative of the potential on the boundary, we take the gradient of eq 22 for the electric field and let r approach the surface a from inside region III. The normal component of the field is

⎞⎤ ∂G(r; r′) σ(r′) dr′⎟⎥ a ∂n(r) ⎠⎥⎦

∂φ̃− ∂φ̃+ σ = − ϵIII ∂n ∂n

∫Γ

q total ∂φ(r+) = − dA ∂n(r+) ϵ

(35)

(assuming that the dielectric constant is ϵ everywhere on Γ).

: F


Article

Journal of Chemical Theory and Computation Notice that the integral of eq 35 is independent of both the ion radius and the surface Γ; however, in real water and in the SLIC model, this condition is not necessarily satisfied. In real water, the field in the first few solvation shells exhibits complicated oscillations;31 the failure mode for the SLIC model is readily apparent for the Born ion. For the standard PB model, the ion radius Rnat is chosen such that the Poisson-model surface charge density σnat reproduces the correct solvation free energy; furthermore, by definition, the potential field outside the ion satisfies eq 35. By contrast, the SLIC radius is a nominal or effective radius, Reff, and leads to an SLIC effective surface charge density σeff. For the surface charge densities in the two models to generate the same solvation free energy, regardless of whether it is taken f rom experiment or explicit-solvent MD, the four quantities must satisfy σeff =

φIII(ra) = φII(ra) f (En(r −a ))

∂n

⎞ ∂φ (rb) dA ⎟ I ⎠ ∂n

∂φII(ra) ∂n

(38) (39)

∫b

∂φII(rb) ∂n

⎞ ∂φ (rb) dA ⎟ I ⎠ ∂n

(40)

Because the potential in each region still satisfies a linear PDE, we may still invoke Green’s theorem as in the standard electrostatic model. Applying the boundary conditions above, and defining

This relation implies that the total SLIC surface charge varies as a function of the ratio of natural and nominal ion radii; in turn, ∂φ this means that the integral ∫ ∂n dA depends on the two Γ radiiin contradiction of eq 35. We remedy this with a correction term that ensures that the exterior potential obeys Gauss’s law outside the solvation shell. Resolving Angstrom-scale and sub-Angstrom scale details of solvent structure in the hydration shell will require much more sophisticated models; the purpose of SLIC is to find a maximally simple continuum theory that reproduces charging free energies and related properties. Therefore, to preserve long-range macromolecular interactions, we impose the correction at the Stern boundary, which is approximately one water layer away from the dielectric boundary; we have used a Stern layer width of 2 Å as commonly employed. The correction replaces the standard Stern-boundary displacement-field condition (eq 9) with ∂φII(rb)

∂n

= (1 + f (En(r −a )))

⎛ q total ⎞ ∂φ (rb) ⎛ ⎜⎜ − ⎟⎟ II =⎜ ϵ ∂n ⎝ ⎝ I ⎠

d1 = −

∫b

∂φIII(ra)

φII(rb) = φI(rb)

R nat σnat R eff

⎛ q total ⎞ ∂φ (rb) ⎛ ⎜⎜ − ⎟⎟ II =⎜ ⎝ ⎝ ϵ+ ⎠ ∂n

(37)

d2 =

q total ϵII

∫b

(41)

∂φII(rb) ∂n

dA

(42)

we obtain ⎤ ⎡⎛ 1 ⎞ − G III,aa ⎥⎡ φIII(ra) ⎤ ⎢⎜⎝ I + KIII,aa⎟⎠ ⎥ ⎥⎢ ⎢ 2 ⎥⎢ ∂φ (ra) ⎥ ⎢⎛1 f ⎞ III ⎥ ⎢ ⎜ I − KII,aa⎟ + G II,aa + KII,ab − G II,ab ⎥⎢ ⎠ 1+f ⎥⎢ ∂n(ra) ⎥ ⎢⎝2 ⎢ ⎥ ⎥ ⎢ ⎢ φ (r ) ⎥ f ⎛1 ⎞ ⎥ ⎢ −K ⎜ I + K ⎟ + G II,ba − G II,bb ⎢ II b ⎥ II,ba II,bb⎠ ⎥ ⎢ 1 + f ⎝2 ⎥⎢ ∂φ (r ) ⎥ ⎢ ⎢ d1 ⎥⎢ II b ⎥ ⎛1 ⎞ ⎜ I − K ⎟ ⎥⎢⎣ ∂n(rb) ⎥⎦ ⎢ I,bb⎠ + G I,bb ⎝2 d2 ⎦ ⎣ ⎡ qG ⎤ ⎢∑ i III ⎥ ⎥ ⎢ 0 ⎥ =⎢ ⎥ ⎢ 0 ⎥ ⎢ ⎥⎦ ⎢⎣ 0

(43)

Note that in the standard dielectric model, f/(1 + f) does not depend on r or the electric field at r, but is merely a constant ratio of permittivities, i.e. a constant times the identity. For SLIC, however, f/(1 + f) can vary along the boundary, and therefore the scaling no longer commutes with the integral operator. Instead, it must be applied to the normal field before the single-layer potential operators GII,aa and GII,ba. We model the static potential as simply as possible, assuming that it takes a constant value throughout the solute.29,31 The net contribution to the solvation energy is then equal to the solute net charge times the static potential. The value of the constant depends on the specific water model to which SLIC parameters are fit.16,74 In this work, we use φstatic = 10 kcal/mol for all calculations parametrized against explicit-solvent MD calculations29,32 in TIP3P water.75 After solving eq 43 selfconsistently by Picard iteration,32 we can obtain the induced reaction field in the solute via Green’s theorem (eq 58). Including the static potential contribution as well, the electrostatic charging free energy is

(36)

Note that this follows the usual convention in which the dielectric constant is the same on both sides of the Stern boundary. Integrating both sides then shows that as intended, the solution satisfies Gauss’s law just outside the surface. The choice of an appropriate surface is empirical, and we use the Stern boundary to limit model complexity (because using a different width shell for the renormalization would require having a total of three surfaces). Stern layers are often used in LPB calculations to prevent unphysically large ion concentrations, and therefore represent a well-understood surface at which to impose this Gauss-Law correction. Note also that the normalization condition cannot be imposed at the dielectric interface itself. The Born ion problem demonstrates why: it has only a single unknown (σ) but two equations, the SLIC and the Gauss-Law correction.

3. A SLIC MODEL WITH IONIC SCREENING We now present a boundary-integral approach for the improved SLIC model that is suitable for either the mixed-dielectric Poisson problem or an LPB problem with a Stern layer. The governing PDEs in each region are as before in the symmetric model; however, the boundary conditions for the new SLIC model are

ΔGcharging =

= G

1 2

1 T induced q φ + q Tφstatic 2

∑ qiφinduced(ri) + φstatic ∑ qi

(44)

(45)


Article


Newton method, may be required for more complicated nonlinear boundary conditions. 3.3. Implementation on an Existing Fast Solver. By design, the SLIC+LPB formulation is very similar to the original Yoon-Lenhoff formulation shown in eq 11, which makes it easy to port to an existing code base. To illustrate its incorporation into existing PB solvers and test the model using atomistically detailed protein models, we implemented it into the GPU-accelerated open-source code PyGBe.86 PyGBe numerically solves the Yoon-Lenhoff formulation using flat triangular panels with piecewise-constant basis functions and centroid collocation. It is well-known that dense BEM solvers require O(N2) memory and O(N2) time for the matrix-vector product calculated at every GMRES iteration. We perform this operation in O(N log N) time and memory using a treecode algorithm.87 The treecode computes the integrals with Gauss quadrature rules, and approximates the influence of far-away Gauss nodes on a center of expansion using Taylor series. This acceleration allows us to run problems with one hundredthousand to one million panels in reasonable computational time. We get further speed-up by offloading computationally intensive parts of the treecode algorithm to a graphics processing unit (GPU) card. The solver implementation was validated by comparison against the MATLAB panel BEM implementation.

Section 5 includes a more detailed discussion of the justification for using a linear-response-type expression for reaction potential. 3.1. Computing the Electric Field Inside the Protein− solvent Interface. SLIC requires computing the normal electric field just inside the dielectric boundary. Theoretically, this can be accomplished easily using the gradient of the potential from Green’s theorem and then the usual limit as the field point approaches the surface. However, this approach involves the hypersingular operator W, which adds substantial complexity to model analysis, discretization, and computation.76−79 To avoid using W, we compute the normal electric field using a simple alternative: after solving the appropriate Green’s-theorem problem for the surface potentials and their normal derivatives, we find an effective surface charge on the dielectric boundary that generates the surface potential (more precisely, we solve an interior Dirichlet problem for each solute using a first-kind single-layer formulation80). For example, in Figure 1, we take the full Cauchy data on a (the surface potential to be matched and its normal derivative) from eq 43 and solve G III,aaσ D = G III,aa

⎛ 1 ⎞ − ⎜ − I + KIII,aa⎟φIII ⎝ ⎠ 2 ∂n

∂φIII

(46)

for the effective surface charge distribution σD. We can then calculate the desired normal electric field using eq 21, that is, using the same method that we used previously for the PCM formulation.32 3.2. Numerical Implementation. We tested the new model using two types of boundary-element method (BEM) discretization. The first was the surface-point-charge BEM used in the original work,32 in which the surface distributions were approximated as a set of discrete point sources located at the mesh vertices; each point was assigned a weight equal to onethird of the total area of the triangles associated with the given vertex.32 The second scheme involved standard panel-based BEM using triangular panels, piecewise-constant basis functions, and centroid collocation, a discretization method we have analyzed previously.65,66,81,82 For centroid-collocation BEM using the Yoon-Lenhoff formulation, the nonlinear operator h(En) is easy to add because each matrix entry requires only evaluating the normal electric field at the relevant panel centroid. Regardless of discretization scheme, the linear systems were solved using preconditioned GMRES83 to a tolerance of 10−5. Following our previous work, the preconditioner P for a given problem was defined as P = D−1, where D is a block matrix where each nonzero block is diagonal and its entries are the corresponding entries of the actual BEM matrix A. This preconditioner is fast to compute (P is a sparse matrix and factorization produces zero fill-in), and has been effective in previous work.57,84,85 Following the earlier SLIC study, we solve the nonlinear BEM problem using the Picard iteration, that is, a sequence of linear problems A(i)x(i) = b, with A(i) defined in terms of x(i−1) (and an initial guess of x(0) = 0, although alternatives might accelerate convergence). Tests showed that three iterations were typically required for the solution to converge; reported results indicate the solution after five iterations. More rigorous analysis of convergence, stopping criteria, and the initial guess are an area of current work, along with questions of solution existence and uniqueness. More sophisticated nonlinear solution techniques, for example, the

4. RESULTS Solute−solvent and Stern-layer surfaces for nonspherical solutes were discretized using methods previously described.32,57,65,66,84 The solute−solvent interface was defined to be the solvent-excluded (molecular) surface and the Sternlayer surface was defined to be a solvent-accessible surface in which the solvent probe radius was set to 2.0 Å.88 These surfaces were discretized into planar triangles using MSMS.89 The MD FEP results used as reference data here are taken from our previous work29,32 or using MD calculations following the protocols described. Where reference data are taken from explicit-solvent MD simulations in TIP3P water, the static potential φstatic was taken to be 10 kcal/mol uniformly throughout the solute, which approximates the static potential observed for various water models.16,31 4.1. Simple Model Problems: Spherical Solutes and the Salt-Dependence of Solvation Free Energies. To understand the effects of including the surface potential and charge renormalization surface, we conducted one set of SLIC calculations using the same atomic radii and SLIC parameters as established previously:32 α = 0.5, β = −60, and γ = −0.5 for the SLIC model, φstatic of 10 kcal/mol, and all solute radii set to 0.92 times their MD Rmin/2. Figure 3 plots SLIC and MD charging free energies for the monatomic ions available in the CHARMM force field.90−92 The proposed SLIC model (thick red and blue lines) tracks the MD results (symbols) more accurately than our original SLIC model (thin red and blue lines). The improvement is especially prominent for the smallest ions, where the FEP calculations alternately overestimate and underestimate the predictions of the new model. Future work will investigate whether discrete solvent packing effects, which play important roles in determining solvation energetics, are responsible for the oscillatory differences. Figure 4 plots the change in ion solvation free energy as a function of salt concentration, in both the standard LPB and the new SLIC model. The results show that SLIC+LPB predicts virtually identical salt dependH


Article


due almost entirely to the static potential. Charging free energies for atoms near and at the surface are also well reproduced, which indicates that SLIC with the Gauss-Law correction contributes the appropriate solvent-steric contribution beyond static-potential asymmetry. Overall, the results in Figure 5 illustrate that SLIC+LPB accurately reproduces atomic charging free energies throughout the range of common biomolecular charges, regardless of an atom’s proximity to solvent. For this example, calculations used a point-BEM approach with a sphere radius of 5.067 Å,29 which fit the central q = −1 charge MD result after removing the static-potential contribution (results did not change significantly if one fit the radius using q = +1). The dielectric and Stern surfaces were discretized using surface point charges with densities of 2 vertices/Å2 and 1 vertex/Å2, respectively. Figure 6 plots panel-BEM calculations of charging free energies of protonated and unprotonated forms of titratable amino acids. The blue circles labeled “SLIC extrapolated” represent Richardson estimates of the converged solution, that is, at infinitely fine discretization, assuming linear convergence with respect to the number of boundary elements.86,93 These results agree well in almost all cases, with the largest deviations observed for the protonated forms of arginine and cysteine, and the deprotonated form of tyrosine. The origins of these discrepancies are not yet clear and warrant further investigation. We next address the addition of physiologically relevant concentrations of salt to the SLIC model, by evaluating the saltdependence of shifts in protonation free energies (more precisely, the solvation contribution, that is, without the Coulomb term). Figure 7 contains plots of the SLIC and symmetric models’ predictions of the change in the solvation component of the protonation free energy, between 145 mM salt and no salt:

Figure 3. Charging free energies for Born ions, computed using explicit-solvent MD FEP (symbols),32 the new SLIC model (thick solid red and blue lines), and the original SLIC model (thin dash-dot red and blue lines).32 The black dashed line represents the Bornenergy prediction from the standard (charge-sign-symmetric) Poisson model.

ΔΔΔG prot,salt = (ΔG prot,145mM − ΔGdeprot,145mM) − (ΔG prot,0mM − ΔGdeprot,0mM)

(47)

The symmetric-model results were calculated using our FFTSVD BEM solver.57 The figure shows that the two models predict energetic changes of comparable magnitude, with the SLIC model predicting changes of slightly increased magnitude for four residues (aspartic acid, cysteine, glutamic acid, and lysine). The symmetric and SLIC model predict very similar changes for histidine and tyrosine. The major difference between the models is arginine, for which the SLIC theory predicts salt to increase the solvation component of the protonation free energy. Further investigation will be needed to identify if this discrepancy is chemically or biologically meaningful, because the solvation component is only one aspect of salt-dependent protonation; conformational changes and changes in Coulomb interactions may well play larger roles. Notably, arginine also exhibits the largest deviations for the charged form, so future work will investigate whether detailed parametrization of atomic radii resolves this discrepancy. 4.2. Calculations and Detailed Assessment on Biomolecules. The first study of SLIC parametrized the model using explicit-solvent FEP results on fictitious “rod” and “bracelet” molecules from Mobley et al.,74 and to this point in the paper we have used these parameters. However, to assess the improved SLIC model against more detailed FEP simulations, we reparameterized SLIC using 11 real compounds: six amino-acid side-chain analogues (ethanamide, methyl ethyl sulfide, acetic acid, propanoic acid, 3-methyl-1h-

Figure 4. Salt-dependence of the charging free energy for a Born ion of radius 1 Å, computed using the new SLIC model (red and blue symbols), and using standard LPB theory (black symbols). Results are plotted as functions of the square root of the ionic strength I (M), and relative to their charging free energies in nonionic solvent, that is, when I = 0.

ence for Born ions, with deviations well under 0.1 kcal/mol between the theories. This result is not surprising, because the system’s symmetry, combined with our renormalization condition, cancels the effect of the asymmetry outside the Stern surface, just not at the dielectric boundary itself. We next demonstrate that SLIC+LPB correctly models the absolute and relative contributions to asymmetry due to the complementary surface potential and finite-size effects. As an example, we calculated charging free energies for a single charge in a small sphere (radius approximately 5 Å), located at different distances from the center of the sphere. Figure 5 plots the SLIC predictions as well as MD FEP results obtained in a previous work.29 The SLIC free energies agree extremely well with the results from explicit-solvent MD, regardless of the charge’s proximity to solvent. Buried charges exhibit the expected behavior: the difference in charging free energies is I


Article


Figure 5. Comparison of the new SLIC model and MD FEP charging free energies for a single charge located within a simple model system: a pseudosphere comprising CHARMM90 sodium atoms on a Cartesian lattice of 0.5 Å grid spacing, such that every solute atom is within 4 Å of the sphere center.29 The SLIC results are computed by first fitting the continuum-model sphere radius from the charging free energy for a central charge with q = −1 (left-most red circle: the number x in each legend entry denotes the charge’s location as (x, 0, 0)). The SLIC calculations used the same parameters from our previous work,32 and with φstatic = 10 kcal/mol.

Figure 6. Charging free energies for titratable amino acids, computed using explicit-solvent MD FEP,32 using the new SLIC model, and using the standard symmetric model.94 The SLIC results labeled “extrapolated” were obtained using Richardson extrapolation under the assumption that the BEM calculations converge linearly in the number of boundary elements.

Figure 7. Salt-induced change in the difference between the charging free energies of the protonated and deprotonated states of titratable amino acids, between 0 and 145 mM.

indole, and p-cresol) and five monovalent ions (sodium, potassium, rubidium, cesium, and chloride). Optimizing SLIC to meet the MD results from our calculations for the ions, and those of Mobley et al. for the side-chain analogues,95−97 we obtained α = 0.320464, β = −46.951476, and γ = −1.182070, and again we assumed ϕstatic = 10.0 kcal/mol/e. Using FEP simulations with the same protocol as our earlier work, we computed the charging free energies for each individual atom in phenylalanine with neutral blocking groups added. Figure 8 contains plots of the explicit-solvent MD results, along with continuum-model predictions from standard PB (with Roux radii94) and predictions from the new SLIC model. Three aspects of these results are noteworthy. First, SLIC is clearly much more accurate than standard PB: the RMS error for standard PB is 3.05 kcal/mol, whereas the RMS error for SLIC is 0.55 kcal/mol. Second, the standard PB model is known to obtain total molecular charging free energies

Figure 8. Comparison of explicit-solvent MD charging free energies for every individual atom in phenylalanine with neutral blocking groups added.

J


Article

Journal of Chemical Theory and Computation correctly, because that is how it has been parametrized.94 Our calculations show that the model achieves this overall accuracy by overpredicting some self-energies and underpredicting others, leading to a cancellation of error. The implication is that the reaction fields due to these individual charges are not correct, so solvent influence on intramolecular screening is not as accurately predicted as with SLIC. The relevance to biological applications as well as to coarse-grained theories such as GB remain to be explored. Third, due to the static potential, SLIC is able to accurately reproduce the positive charging free energies observed for atoms that are typically considered as hydrophobic. Standard PB theories, no matter how the radius is parametrized, are unable to reproduce positive charging free energies. To assess SLIC against a diverse set of small molecules we computed electrostatic solvation free energies for the Mobley test set,95 leading to a comparison of 502 free energies; the full set includes 504 compounds but for two molecules, ammonia and N,N-dimethyl-p-methoxybenzamide, our surface meshing software failed to generate a consistent surface. The results are plotted in Figure 9. The comparison shows a strong correlation

vertex per square Angstrom (totalling 21 612 boundary elements) to eight vertices per square Angstrom (120 400 boundary elements). Richardson extrapolation confirms that the implementation converges at the expected linear rate, which is the same rate as the standard PB model (Figure 10). As a

Figure 10. Calculations of the SLIC model using the PyGBe/SLIC GPU-accelerated treecode demonstrate the correct order of convergence (estimated final answer obtained using Richardson extrapolation); dashed lines between symbols are guides for the eye.

Figure 9. Comparison of SLIC continuum-model predictions to the explicit-solvent MD charging free energies computed by Mobley et al.95

model of pH-dependent protein behavior, we calculated how the SLIC model changes predictions of solvation free energy changes in a protein due to side-chain protonation or deprotonation. For consistency with the derivation of the SLIC model as an approximation to the potential of mean force obtained by integrating out water degrees of freedom only,1 we used a dielectric constant of 1 in the protein. We assumed a simple reference charge state relevant for physiological pH; that is, aspartic acids and glutamic acids were deprotonated, and lysines, arginines, and tyrosines were protonated. Figure 11 plots the solvation free energy changes for protonation/ deprotonation, calculated using the standard PB model and the SLIC model. For these calculations the surface discretization had eight vertices per square Angstrom. The results show that

(for a vertex density of 4/Å2, the correlation coefficient is R = 0.9653), with an RMSD 2.07 kcal/mol. It is important to note two details that limit the model’s accuracy in this test. First, the present SLIC results were obtained after a parametrization that included six of these compounds, the side chain analogues noted above, and thus the training set omitted many chemical groups found in this set. Second, the atomic radii used in our calculations were not parametrized but simply set to 0.92Rmin/ 2, where the scale factor 0.92 was determined in previous work using Mobley’s test set of fictitious “rod” and “bracelet” molecules.32 Future work will investigate the dominant sources of error and develop a comprehensive parametrization. For an initial, exploratory investigation of the possible importance of the SLIC model in protein electrostatics, we calculated solvation free energies for bovine pancreatic trypsin inhibitor in multiple protonation states. These calculations involved much larger BEM problems and therefore required use of the modified PyGBe GPU treecode.86 The protein structure was taken from PDB accession code 3BTK98 and prepared in earlier work.99 The discretizations ranged from one

Figure 11. Changes in electrostatic solvation free energy due to changes in protonation state, using the standard PB model with SLIC radii, and the SLIC model with parameters optimized to reproduce explicit solvent calculations on Born ions and side-chain analogues. K


Article


For water, the optimized values were α = 0.626, β = −112.0, γ = −0.93, ϕstatic = −12.1 kcal/mol/e, γA = 0.0037 kcal/mol/Å2, and γB = 2.75 kcal/mol. Our recovery of a static potential from experimental measurements that is opposite in sign to the static potential recovered from semiclassical explicit-solvent models is in agreement with ab initio simulations.105 The values for the electrostatic parameters are close to parameters determined to fit other data, for example, molecular charging free energies, and the nonpolar parameters are comparable to those determined by Honig et al.102 (γA = 0.005 kcal/mol/Å2 and γB = 0.86 kcal/mol). Over a set of 152 neutral small molecules for which the MNSol database has measurements for both water and octanol, the resulting SLIC model for water has an RMS error of 1.78 kcal/mol (mean absolute error of 1.47 kcal/ mol); we note that this set includes the nine neutral compounds in the training set. The calculated results are compared to the experimental measurements in Figure 12. For

SLIC leads to moderately reduced changes in free energy compared to predictions from the standard PB model. These preliminary results are promising because this is the same effect as observed by increasing the protein dielectric constant, namely reducing the magnitudes of computed pKa shifts. Popular approaches for calculating accurate pKa values use high protein dielectric constants to correct for the overpredictions of standard PB models. This empirical approach has been successful, and is motivated by the reasonable assumption that protein flexibility plays an important role in deviations between theory and experiment, but high-dielectric models have not fully corrected for the deviations. We have previously suggested that nonbulk response in the solvation layer contributes additional deviations,100 and the present results motivate future work applying SLIC to protein pKa calculations. For each given residue type, deviations between SLIC and standard PB are remarkably similar. The ASP residues have deviations of 3.7 and 3.4 kcal/mol; GLU residues, 3.9 and 4.7 kcal/mol; TYR residues, −6.0, −6.5, −5.3, and −7.2 kcal/mol; ARG residues, −10.4, −12.0, −10.1, −11.8, and −9.1 kcal/mol; and LYS residues, −8.8, −9.0, −10.8, and −9.4 kcal/mol. 4.3. Calculating Water−Octanol Transfer Free Energies. To test the model on another set of practical problems, we calculated the solvation free energies for a set of 152 small molecules in both water and octanol, and then compared both the computed and experimental solvation free energies and water−octanol transfer free energies. Our reference experimental free energies were taken from the MNSol database.101 To ensure consistency we parametrized both the SLIC electrostatic model as well as a popular but simplistic model for the nonpolar component of the solvation free energy, in which ΔGnp = γAA + γB where A is the solute’s surface area (as given in the MNSol database) and γA and γB are fitting parameters. Note that in contrast to earlier work using this model, parameters for our nonpolar and electrostatic components are determined simultaneously. This difference is important because our electrostatic model is better able to capture the electrostatic charging contribution to solvation free energies, including the unfavorable contributions of hydrophobic groups (Figure 8). In contrast, earlier work has noted that the nonpolar parameters were determined by assuming that the solvation free energies of molecules in the test set (usually alkanes) have zero electrostatic component.102,103 Overall, the complete SLIC implicit-solvent model used here had six free parameters: the four SLIC parameters α, β, γ, ϕstatic, and the two nonpolar parameters γA and γB. These were optimized using nonlinear least-squares to reproduce solvation free energies for three Born ions (sodium, potassium, and chloride) and nine small molecules (acetic acid, ethanol, methanol, p-cresol, propanoic acid, toluene, n-octane, ethylamine, and 1,4-dioxane). These reference compounds satisfied the following desirable properties: (1) cover a substantial representation of the chemical diversity observed in protein side chains; (2) have readily available atom radii, geometries, and charge distributions from MD parameter sets and the Mobley test set;95,96 and (3) have solvation free energies in both water and octanol. The ion solvation free energies in octanol were obtained from Zhao and Abraham.104 The optimizations in each solvent started from the same point and used eight iterations with Matlab’s nonlinear least-squares function, with all optimization parameters set to their default values.

Figure 12. Calculated solvation free energies for 152 neutral compounds in water, compared to experimental measurements available in the MNSol database.101 The RMS error is 1.78 kcal/mol and the mean absolute error is 1.47 kcal/mol.

octanol, the optimized parameters were α = 0.364, β = −108.5, γ = −1.58, ϕstatic = −5.93 kcal/mol/e, γA = −0.0214 kcal/mol/ Å2, and γB = 3.96 kcal/mol. The nonpolar parameters compare reasonably well to those of Park et al.,103 who used the parametrization approach of Honig et al. for octanol and found γA = −0.0328 kcal/mol/Å2 and γB = 1.92 kcal/mol (Park et al. also determined, for water, γA = 0.0086 kcal/mol/Å2 and γB = 1.37 kcal/mol). Over the same test set of 152 neutral small molecules, the RMS error in octanol is 2.31 kcal/mol and the mean absolute error is 1.93 kcal/mol (Figure 13). In Figure 14 we plot the calculated transfer free energies ΔGwater − ΔGoctanol for these 152 compounds, against the experimental transfer free energies. The RMS error is surprisingly small, 1.07 kcal/mol, indicating a favorable cancellation of errors (the mean absolute error is 0.88 kcal/mol). The red asterisks in Figure 14 represent the solutes with AMBER atom type “O”. These have been highlighted because these solutes tend to have overpredicted solvation energies (data not shown), but the transfer energies are well predicted. More detailed analysis is required to understand the origins of the systematic errors. L


Article


atomic charge δq produces a proportional increment to the potential (and the constant of proportionality is set by h). The parameter β controls the width of the transition region between two linear-response regimes, so as β → ∞, we obtain that the 1 reaction component of the charging free energy is 2 ϕreacq, where ϕreac is the potential at the final q. A more detailed explanation of the charging process, and the justification for using the final ϕreac and omitting the transition region of h, may be given as follows. We perform an FEP charging calculation between λ = 0 (uncharged) and λ = 1 (fully charged), so that for a given λ, an atom i whose wild-type charge is qi will have a charge of λqi. In this process, the direct Coulomb electric field at any λ is simply the final electric field (when λ = 1), scaled by λ. We now invoke the Coulomb-field approximation, which assumes that the reaction field is proportional to the direct Coulomb field, but opposite in sign and smaller in magnitude.106,107 As a result, for any finite λ, the total field, which equals the sum of the direct Coulomb field and the reaction field, also takes the sign that it takes at λ = 1. Last, we assume that the SLIC is a step function; this is the limiting case as β → ±∞, which is a reasonable approximation because β has been found to be large. Critically, because the SLIC is a step function, the system responds linearly over the charging process. This can be seen in Figure 2: where SLIC is horizontal, the surface charge varies linearly with the local electric field. As a consequence, the reaction potential also varies linearly with λ over the charging process, except at λ = 0. The charging integral (over λ) therefore takes the usual linearresponse form, and we can write the polarization contribution 1 to the charging free energy as 2 ∑ qiφIII(ri), where φIII(ri) is the reaction potential computed at the ith charge location, given the wild-type charge distribution (i.e., when λ = 1) (for more discussion, see refs 32, 36, and 60.) We have emphasized that SLIC is determined by the normal electric field just inside the dielectric boundary. Standard dielectric theory holds that for a sufficiently smooth surface (one for which each point on the surface has a unique tangent plane), the normal electric field changes discontinuously, with half the change occurring as one moves from the limiting value inside to the surface, and the remaining half as one moves from the surface to the limiting value outside.73 This treatment can be justified from the continuum limit in which representative volume elements of the material have zero size.4 The present work focuses on problems for which this assumption is not justified, and therefore we developed the present model to avoid the issue of “where” the field changes as one moves from one side of the interface to the other. Confirming this physical interpretation, computational tests that replace the interior limit (eq 20) with the actual surface value (eq 24) exhibited substantially poorer results. We have shown that the present single-layer approach suffices for developing LPB-based models involving SLIC. However, alternative BIE formulations for interior Dirichlet problems108 may be preferable for large-scale calculations or for other nonlinear boundary conditions.

Figure 13. Calculated solvation free energies for 152 neutral compounds in octanol, compared to experimental measurements available in the MNSol database.101 The RMS error is 2.31 kcal/mol and the mean absolute error is 1.93 kcal/mol.

Figure 14. Calculated transfer free energies for 152 neutral compounds in octanol, compared to experimental transfer free energies computed from the MNSol database.101 The RMS error is 1.07 kcal/mol and the mean absolute error is 0.88 kcal/mol. The red asterisks represent solutes containing AMBER atom type “O”.

5. DISCUSSION SLIC would seem to complicate the use of simple expressions for the charging free energy. Whereas linear-response theories 1 can express energies via quadratic expressions such as 2 Lq2 for a charge q, in general a nonlinear electrostatic model necessitates that the charging free energy be determined by integrating over a specified charging process. However, in the present work we have used an approximated expression for the solvation free energy, which seems to contradict our premise of a nonlinear model and instead imply linear response. Here we discuss the justification for this approximation. First, studies have shown that FEP charging processes that scale all of the charges equally (the Kirkwood charging process) lead to energies that do change essentially quadratically as a function of the FEP parameter λ.29,94 Second, notice that the nonlinear perturbation h(En) approaches a step function as the parameter β → ∞. For any portions of a charging process in which h is constant, the system responds linearly, that is, an incremental

6. SUMMARY We have extended our solvation-layer interface condition (SLIC) model to address electrolyte solvents using the linearized Poisson−Boltzmann (LPB) equation. The new model includes a second modified boundary condition at the Stern surface, forcing Gauss’s law to be satisfied in the far field. The new model also includes the contribution of the static M


Article

Journal of Chemical Theory and Computation potential φstatic to the charging free energy. The static potential component contributes a component of the charging free energy that is linear in the solute net charge; for neutral compounds, it therefore contributes zero. This static potential is solvent-model dependent, that is, its value should be determined consistently with the other SLIC model parameters, or at the least determined using consistent explicit-solvent simulations, if explicit-solvent results are used as the basis for parametrization.32 For example, the initial model validation results reported here are compared to explicit-solvent simulations in TIP3P water, which is consistent with our prior determination of the static potential in TIP3P water.29 However, for other solvents, including polarizable waters, the value of the field would need to be determined using either model systems that isolate its value, or treat it as an additional fitting parameter while optimizing the model to match solvation free energies. On the other hand, when parametrizing SLIC against experimental solvation free energies, the static potential should be treated as an empirical parameter. Furthermore, the model has been assumed to be a constant, which is consistent with observations of potentials in completely uncharged solutes.29,31,32 Future work should assess the variations in more detail. Overall, the new SLIC+LPB model includes both of the dominant physical mechanisms underlying chargehydration asymmetry. We have used boundary-integral equations (BIEs) and the boundary-element method (BEM) to solve this model numerically, but it could easily be implemented in volumetric solvers. A fast GPU treecode has been implemented to accelerate the model’s calculation and illustrate that hydration-shell response contributes several kcal/ mol to protonation energetics, making it an important factor to consider when modeling protein titration. Our nonlinear model satisfyingly recovers the standard macroscopic continuum model under the appropriate limits. Because the asymmetry magnitude represents the different abilities of water hydrogens and oxygens to screen the local field, future work will explore the interpretation of a “local” solvent dielectric constant that depends on the field and or curvature. For example, in this work we have treated the parameter μ as being fixed by the other model parameters (specifically, μ = −α tanh(−γ)) to ensure that the macroscopic dielectric boundary condition is satisfied for weak Coulomb fields. However, the parameter can be interpreted as the dielectric constant in the solvation shell, and hence letting it be a tunable parameter is an important question for future work. It is also interesting that the simple BIE modification eq 17 gives rise to a nonintuitive displacement boundary condition, suggesting that BIEs may have some advantages over PDE approaches when developing nonlinear models for hydration shell response. For instance, long-range interactions required introducing the second charge layer. Its physical interpretation is phenomenological at the present time; it is possible that this feature may improve predictions of pairwise interactions when molecules approach each other in encounter complexes. In these circumstances, the separate normalizing charges disappear as the Stern layers begin to overlap, changing the overall interaction. Numerous computational extensions are also possible, which would enhance speed and or efficiency for specific types of calculations. For example, to enhance speed one could derive a novel boundary-integral approach mixing the surface-charge model for the dielectric boundary, with a Green’s theorem approach at the Stern layer. It would also be valuable to

implement the purely second-kind BIE formulation of Juffer,109 and the Bordner−Huber110 model for multiprotein solutions. Solutes with overlapping Stern layers can be treated using SLIC versions of existing formulations.57 The hypersingular operator poses mathematical and computational challenges, which will require careful analysis and suitable basis and testing functions.67,111 We note that the use of modified boundary conditions has been explored recently to accelerate calculations of the linear Poisson−Boltzmann model, in which calculations involving a Stern layer (meaning both a dielectric interface and a Stern surface) can be replaced with a single-surface model with a modified boundary condition at the dielectric interface.58 Also, the modeling discussed here focuses on gaining insights via boundary-integral equations; further work will be needed to assess whether the modified boundary conditions affect relative performance of volumetric numerical methods (e.g., finiteelement methods) and boundary based methods.

■

APPENDIX A: BOUNDARY-INTEGRAL METHODS FOR CONTINUUM SOLVENT MODELS In this paper we consider potential fields satisfying either a Laplace-type PDE such as eq 2, or a linearized Poisson− Boltzmann (LPB) equation such as eq 1. To reformulate these coupled PDE problems using boundary-integral equations (BIEs), we make extensive use of the free-space Green’s functions for these equations. Each Green’s function represents the potential due to a discrete point charge (Dirac delta function). The Laplace Green’s function satisfies ∇2 G Laplace(r; r′) = −δ(r − r′)

(48)

where r′ is the location of the charge and r is the point in space, and the solution for zero boundary conditions is the familiar Coulomb potential: G Laplace(r; r′) =

1 4π∥r − r′∥

(49)

The LPB free-space Green’s function is known as the Yukawa or screened Coulomb potential: G LPB(r; r′) =

exp( −κ∥r − r′∥) 4π∥r − r′∥

(50)

The boundary-integral operators associated with these Green’s functions map either a single layer of surface charge, or a double layer, to the resulting potential at a surface, or normal gradient. The “source” and “destination” surfaces may be the same surface, or different surfaces, for example, a and b, or Stern layers of two proteins. For a free-space Green’s function G*(r;r′), the single-layer potential operator V*,xy maps a monopole charge distribution σ(r) on the boundary y to the resulting potential on the boundary x (note the subscript indices are ordered following matrix notation); that is, for r on x, we write ϕx(r) =

∫y G*(r, r′)σy(r′) d2r′

(51)

as

ϕx = V , xyσy (52) * Similarly, the double-layer potential operator K*,xy maps a dipole charge (a charge double-layer) on the boundary y, call it μy(r), to the potential on the boundary x, so that N


Article

Journal of Chemical Theory and Computation ψx(r) =

∫ ̵ y ∂G∂*n((rr;′)r′) μy (r′) d2r′

Green’s functions using a subscript denoting the appropriate region of space, we obtain

(53)

∂φIII ⎛1 ⎞ ⎜ I + K ⎟ = III,aa ⎠φIII − G III,aa ⎝2 ∂n

can be written more compactly as ψx = K

μ *,xy y

∂n(r)

∫̵ y

=

We next invoke Green’s theorem in region II, obtaining

= n(̂ r) ·∇

φII(r)=

∂n

+

2

(57)

φII(ra) = −

because the operator is adjoint to the double-layer potential operator K*,xy in an L2 sense. The normal gradient of the double-layer potential, denoted as W*,xy, poses special challenges analytically and numerically because of the increased order of the singularity.67,80,111 In this work, we describe a BIE formulation for SLIC that does not require using the hypersingular operator. We then invoke Green’s theorem in each region, take limits as the field point approaches the region’s bounding surface(s), and apply boundary conditions. For r in region III,

−

∫a

0=

∑ qiG III(r, ri)

∂φ (r′)

∫a GII(ra, r′) ∂nII(r′) ∫a

dr′

⎞ ∂G II(ra, r′) φII(r′) dr′⎟ ∂n(r′) ⎠

∂φ (r′)

+

∫b GII(ra, r′) ∂nII(r′)

−

∫b

dr′

∂G II(ra, r′) φII(r′) dr′ ∂n(r′)

(62)

∂φII(ra) ⎛1 ⎞ ⎜ I − K ⎟ + KII,abφII(rb) II,aa ⎠φII(ra) + G II,aa ⎝2 ∂n(ra) − G II,ab

In the limit as r → ra from inside (i.e., the side opposite the normal direction), then, we obtain

∫a GIII(ra, r′)

(61)

Rearranging terms,

(58)

φIII(ra) =

∂G II(r, r′) φII(r′) dr′ ∂n(r′)

⎛1 + ⎜ φII(ra) + ⎝2

∂φIII(r′)

dr′ ∂n(r′) ∂G III(r, r′) φIII(r′) dr′ + ∂n(r′)

∫b

dr′

This equation leads to two boundary-integral equations, one for the limit r → a and one for r → b. The signs on the first two terms of eq 61 reflect the fact that Green’s theorem is usually expressed for a potential inside a surface with the unit normal defined pointing outward. However, the normal on a points into region II (Figure 1). Taking the limit as r → a, we obtain

(56)

= K ′ ,xyσy *

∫a GIII(r, r′)

∂φ (r′)

∫b GII(r, r′) ∂nII(r′)

(55)

∂G (r; r′) * σy(r′) d2r′ ∂n(r)

φIII(r) =

∫a

−

∫y G*(r; r′)σy(r′) d r′

⎛ ∂φ (r′) ⎞ II ⎟ dr′ ∂n(r′) ⎠ ⎛ ∂G II(r, r′) ⎞ ⎜− ⎟φ (r′) dr′ ∂n(r′) ⎠ II ⎝

∫a GII(r, r′)⎜⎝− −

and we denote it as

∂ϕx

(60)

(54)

The symbol ∫ ̵ in eq 53 indicates that if x and y represent the same boundary, the boundary integral is to be understood in the Cauchy principal-value sense, to account for the singularity of the normal derivative of the Green’s functions as r′ → r. Some boundary-integral formulations also require the operator that maps a single-layer surface charge distribution to the normal gradient of the resulting potential. From eq 51 we see that the normal gradient of the single layer potential ϕx = V*,xyσy is ∂ϕx(r)

∑ qiG III

∂φII(rb) ∂n(rb)

(63)

We remind the reader that KII,ab denotes the double-layer operator that maps a double-layer distribution on b to the surface a, and the II subscript makes explicit that the associated Green’s function for the operator is the one associated with region II. Taking the limit of eq 61 as r → b,

⎛1 dr′ + ⎜ φIII(ra) − ∂n(r′) ⎝2

∂φIII(r′)

⎞

∫ ̵ a ∂G∂IIIn′((rra,′)r′) φIII(r′) dr′⎟⎠ + ∑ qiGIII(ra, ri)

φII(rb) = −

(59) 1 φ 2

In eq 59 the new term arises because the double-layer potential has a discontinuity of magnitude φ(ra) as r crosses the boundary at ra in the normal direction.73,108 The factor of 1 2 results from our assumption that the surface is sufficiently smooth (has a continuous normal) at ra, so that the surface can be considered as locally planar and a sphere of vanishing radius centered at ra has equal solid angles inside and outside the surface. Rearranging terms and specifying the operators’

∂φ (r′)

∫a GII(rb, r′) ∂nII(r′) d r′

+

∫a

+

∫b

∂G II(rb, r′) φII(r′)d r′ ∂n(r′) ∂φ (r′) G II(rb, r′) II d r′ ∂n(r′)

⎛1 + ⎜ φII(rb) − ⎝2

∫b

⎞ ∂G II(rb, r′) φII(r′)d r′⎟ ∂n(r′) ⎠

(64)

so O


Article

Journal of Chemical Theory and Computation 0 = − KII,baφII(ra) + G II,ba +

(7) Vlachy, V. Ionic effects beyond Poisson-Boltzmann theory. Annu. Rev. Phys. Chem. 1999, 50, 145−165. (8) Grochowski, P.; Trylska, J. Continuum molecular electrostatics, salt effects, and counterion binding−a review of the PoissonBoltzmann theory and its modifications. Biopolymers 2007, 89, 93− 113. (9) Borukhov, I.; Andelman, D.; Orland, H. Steric effects in electrolytes: a modified Poisson-Boltzmann equation. Phys. Rev. Lett. 1997, 79, 435−438. (10) Gong, H.; Freed, K. F. Langevin-Debye model for nonlinear electrostatic screening of solvated ions. Phys. Rev. Lett. 2009, 102, 057603 DOI: 10.1103/PhysRevLett.102.057603. (11) Hu, L.; Wei, G.-W. Nonlinear Poisson equation for heterogeneous media. Biophys. J. 2012, 103, 758−766. (12) Marten, B.; Kim, K.; Cortis, C.; Friesner, R. A.; Murphy, R. B.; Ringnalda, M. N.; Sitkoff, D.; Honig, B. New model for calculation of solvation free energies: correction of self-consistent reaction field continuum dielectric theory for short-range hydrogen-bonding effects. J. Phys. Chem. 1996, 100, 11775−11788. (13) Gallicchio, E.; Paris, K.; Levy, R. M. The AGBNP2 implicit solvation model. J. Chem. Theory Comput. 2009, 5, 2544−2564. (14) Latimer, W. M.; Pitzer, K. S.; Slansky, C. M. The free energy of hydration of gaseous ions, and the absolute potential of the normal calomel electrode. J. Chem. Phys. 1939, 7, 108−112. (15) Lynden-Bell, R. M.; Rasaiah, J. C.; Noworyta, J. P. Using simulation to study solvation in water. Pure Appl. Chem. 2001, 73, 1721−1731. (16) Rajamani, S.; Ghosh, T.; Garde, S. Size dependent ion hydration, its asymmetry, and convergence to macroscopic behavior. J. Chem. Phys. 2004, 120, 4457. (17) Grossfield, A. Dependence of ion hydration on the sign of the ion’s charge. J. Chem. Phys. 2005, 122, 024506. (18) Mukhopadhyay, A.; Fenley, A. T.; Tolokh, I. S.; Onufriev, A. V. Charge hydration asymmetry: the basic principle and how to use it to test and improve water models. J. Phys. Chem. B 2012, 116, 9776− 9783. (19) Vath, P.; Zimmt, M. B.; Matyushov, D. V.; Voth, G. A. A failure of continuum theory: temperature dependence of the solvent reorganization energy of electron transfer in highly polar solvents. J. Phys. Chem. B 1999, 103, 9130−9140. (20) Bockris, J. O.; Devanathan, M. A. V.; Müller, K. On the structure of charged interfaces. Proc. R. Soc. London, Ser. A 1963, 274, 55−79. (21) Kang, Y. K.; Némethy, G.; Scheraga, H. A. Free energies of hydration of solute molecules. 1. Improvement of the hydration shell model by exact computations of overlapping volumes. J. Phys. Chem. 1987, 91, 4105−4109. (22) Yu, Z.; Jacobson, M. P.; Josovitz, J.; Rapp, C. S.; Friesner, R. A. First-shell solvation of ion pairs: Correction of systematic errors in implicit solvent models. J. Phys. Chem. B 2004, 108, 6643−6654. (23) Corbeil, C. R.; Sulea, T.; Purisima, E. O. Rapid prediction of solvation free energy. 2. The first-shell hydration (FiSH) continuum model. J. Chem. Theory Comput. 2010, 6, 1622−1637. (24) Bonthuis, D. J.; Netz, R. R. Beyond the continuum: how molecular solvent structure affects electrostatics and hydrodynamics at solid-electrolyte interfaces. J. Phys. Chem. B 2013, 117, 11397−11413. (25) Das, S.; Chakraborty, S.; Mitra, S. K. Redefining electrical double layer thickness in narrow confinements: effect of solvent polarization. Phys. Rev. E 2012, 85, 051508. (26) Rashin, A. A.; Honig, B. Reevaluation of the Born model of ion hydration. J. Phys. Chem. 1985, 89, 5588−5593. (27) Purisima, E. O.; Sulea, T. Restoring charge asymmetry in continuum electrostatic calculations of hydration free energies. J. Phys. Chem. B 2009, 113, 8206−8209. (28) Mukhopadhyay, A.; Aguilar, B. H.; Tolokh, I. S.; Onufriev, A. V. Introducing charge hydration asymmetry into the Generalized Born model. J. Chem. Theory Comput. 2014, 10, 1788−1794. (29) Bardhan, J. P.; Jungwirth, P.; Makowski, L. Affine-response model of molecular solvation of ions: Accurate predictions of asymmetric charging free energies. J. Chem. Phys. 2012, 137, 124101.

∂φII(ra) ∂n(ra)

∂φII(rb) ⎛1 ⎞ ⎜ I + K ⎟ φ (r ) − G II,bb b II,bb II ⎝2 ⎠ ∂n(rb)

(65)

The steps for region I are essentially the same and lead to one boundary-integral equation: 0=

∂φ (rb) ⎛1 ⎞ ⎜ − KI,bb⎟φI(rb) + G I,bb I ⎝2 ⎠ ∂n(rb)

(66)

Applying the boundary conditions eqs 4, 5, 9, and 10, we reduce the number of unknowns from eight to four, so that eqs 60, 63, 65, and 66 lead to the given system of four BIEs with four unknowns57,86 in eq 11.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Christopher D. Cooper: 0000-0003-0282-8998 Jaydeep P. Bardhan: 0000-0002-6051-8499 Funding

M.G.K. was partially supported by the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research, under Contract DE-AC02-06CH11357, and also NSF Grant OCI-1450339. C.D.C. acknowledges support from Conicyt through the Basal Project FB 0821. J.P.B., A.M.T., and A.M.R. have been supported in part by the National Institute of General Medical Sciences (NIGMS) of the National Institutes of Health (NIH) under Award No. R21GM102642, and also NSF Grant CBET-1604369. Notes

The authors declare no competing financial interest. MATLAB source code and surface discretizations for running the nonlinear boundary-condition calculations and generating the figures, are freely and publicly available online at https:// bitbucket.org/bardhanlab/si-slic-lpb.

■

ACKNOWLEDGMENTS The content is solely the responsibility of the authors, and does not necessarily represent the official views of the National Institutes of Health. The authors gratefully acknowledge the use of computing facilities at the Mathematics and Computer Science Division of Argonne National Laboratory, and the Discovery Cluster at Northeastern University.

■

REFERENCES

(1) Roux, B.; Simonson, T. Implicit Solvent Models. Biophys. Chem. 1999, 78, 1−20. (2) Davis, M. E.; McCammon, J. A. Electrostatics in biomolecular structure and dynamics. Chem. Rev. 1990, 90, 509. (3) Zhou, H.-X. Macromolecular electrostatic energy within the nonlinear Poisson-Boltzmann equation. J. Chem. Phys. 1994, 100, 3152−3162. (4) Beglov, D.; Roux, B. Solvation of complex molecules in a polar liquid: an integral equation theory. J. Chem. Phys. 1996, 104, 8678− 8689. (5) Stultz, C. M. An Assessment of potential of mean force calculations with implicit solvent models. J. Phys. Chem. B 2004, 108, 16525−16532. (6) Kollman, P. Free energy calculations: Applications to Chemical and Biochemical Phenomena. Chem. Rev. 1993, 93, 2395−2417. P


Article


cavity in classical and quantum water models. J. Phys. Chem. Lett. 2014, 5, 2767−2774. (51) Rodgers, J. M.; Weeks, J. D. Local molecular field theory for the treatment of electrostatics. J. Phys.: Condens. Matter 2008, 20, 494206. (52) Nguyen, C. N.; Young, T. K.; Gilson, M. K. Grid inhomogeneous solvation theory: hydration structure and thermodynamics of the miniature receptor cucurbit[7]uril. J. Chem. Phys. 2012, 137, 044101. (53) Huggins, D. J. Application of inhomogeneous fluid solvation theory to model the distribution and thermodynamics of water molecules around biomolecules. Phys. Chem. Chem. Phys. 2012, 14, 15106−15117. (54) Pomogaeva, A.; Chipman, D. M. Field-extremum model for short-range contributions to hydration free energy. J. Chem. Theory Comput. 2011, 7, 3952−3960. (55) Rizzo, R. C.; Aynechi, T.; Case, D. A.; Kuntz, I. D. Estimation of Absolute Free energies of hydration using continuum methods: accuracy of partial charge models and optimization of nonpolar contributions. J. Chem. Theory Comput. 2006, 2, 128−139. (56) Bockris, J. O.; Reddy, A. K. N. Modern Electrochemistry: An Introduction to an interdisciplinary area; Plenum Press: New York, 1973. (57) Altman, M. D.; Bardhan, J. P.; White, J. K.; Tidor, B. Accurate Solution of Multi-region Continuum Electrostatic Problems Using the Linearized Poisson-Boltzmann Equation and Curved Boundary Elements. J. Comput. Chem. 2009, 30, 132−153. (58) Manzin, A.; Bottauscio, O.; Ansalone, D. P. Application of the thin-shell formulation to the numerical modeling of Stern layer in biomolecular electrostatics. J. Comput. Chem. 2011, 32, 3105−3113. (59) Kathmann, S. M.; Kuo, I.-F. W.; Mundy, C. J.; Schenter, G. K. Understanding the surface potential of water. J. Phys. Chem. B 2011, 115, 4369−4377. (60) Bardhan, J. P.; Tejani, D.; Wieckowski, N.; Ramaswamy, A.; Knepley, M. G. A Nonlinear Boundary Condition for Continuum Models of Biomolecular Electrostatics. Proceedings of the Progress in Electromagnetics Research Symposium (PIERS). Cambridge, MA, USA, 2015; pp 1215−1221, (Prague, Czech Republic: July 6−9, 2015). (61) Miertus, S.; Scrocco, E.; Tomasi, J. Electrostatic Interactions of a solute with a continuum - a direct utilization of ab initio molecular potentials for the prevision of solvent effects. Chem. Phys. 1981, 55, 117−129. (62) Tomasi, J.; Persico, M. Molecular Interactions in Solution: An Overview of Methods based on continuous descriptions of the solvent. Chem. Rev. 1994, 94, 2027−2094. (63) Shaw, P. B. Theory of the Poisson Green’s-function for discontinuous dielectric media with an application to protein biophysics. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 32, 2476−2487. (64) Zauhar, R. J.; Morgan, R. S. The Rigorous Computation of the Molecular Electric Potential. J. Comput. Chem. 1988, 9, 171−187. (65) Altman, M. D.; Bardhan, J. P.; White, J. K.; Tidor, B. An Efficient and Accurate Surface Formulation for Biomolecule Electrostatics in Non-ionic Solution. Engineering in Medicine and Biology Conference (EMBC). Piscataway, NJ, USA, 2005; pp 7591−7595, (Shanghai, China: September 1−4, 2005). (66) Bardhan, J. P. Numerical solution of boundary-integral equations for molecular electrostatics. J. Chem. Phys. 2009, 130, 094102. (67) Hsiao, G. C.; Kleinman, R. E. Error analysis in numerical solution of acoustic integral equations. Int. J. Num. Methods Eng. 1994, 37, 2921−2933. (68) Bardhan, J. P.; Knepley, M. G. Mathematical analysis of the boundary-integral based electrostatics estimation approximation for molecular solvation: Exact results for spherical inclusions. J. Chem. Phys. 2011, 135, 124107. (69) Rocchia, W.; Alexov, E.; Honig, B. Extending the Applicability of the Nonlinear Poisson-Boltzmann Equation: Multiple Dielectric Constants and Multivalent Ions. J. Phys. Chem. B 2001, 105, 6507− 6514.

(30) Ashbaugh, H. S. Convergence of Molecular and Macroscopic Continuum Descriptions of Ion Hydration. J. Phys. Chem. B 2000, 104, 7235−7238. (31) Cerutti, D. S.; Baker, N. A.; McCammon, J. A. Solvent reaction field potential inside an uncharged globular protein: a bridge between implicit and explicit solvent models? J. Chem. Phys. 2007, 127, 155101. (32) Bardhan, J. P.; Knepley, M. G. Modeling Charge-Sign Asymmetric Solvation Free Energies With Nonlinear Boundary Conditions. J. Chem. Phys. 2014, 141, 131103. (33) Decherchi, S.; Colmenares, J.; Catalano, C. E.; Spagnuolo, M.; Alexov, E.; Rocchia, W. Between algorithm and model: different molecular surface definitions for the Poisson-Boltzmann based electrostatic characterization of biomolecules in solution. Commun. Comput. Phys. 2013, 13, 61−89. (34) Blum, L.; Wei, D. Q. Analytical solution of the mean spherical approximation for an arbitrary mixture of ions in a dipolar solvent. J. Chem. Phys. 1987, 87, 555−565. (35) Fawcett, W. R.; Blum, L. The role of dipole-dipole interactions in the solvation of monoatomic monovalent ions in water on the basis of the mean spherical approximation. J. Electroanal. Chem. 1993, 355, 253−263. (36) Molavi Tabrizi, A.; Knepley, M.; Bardhan, J. P. Generalizing The Mean Spherical Approximation as a Multiscale, Nonlinear Boundary Condition at the Solute-Solvent Interface. Mol. Phys. 2016, 114, 2558− 2567. (37) Molavi Tabrizi, A.; Goossens, S.; Mehdizadeh Rahimi, A.; Knepley, M.; Bardhan, J. P. Predicting solvation free energies and thermodynamics in polar solvents and mixtures using a solvation-layer interface condition. J. Chem. Phys. 2017, 146, 094103. (38) Bardhan, J. P. Biomolecular electrostatics−I want your solvation (model). Comput. Sci. Discovery 2012, 5, 013001. (39) Kovalenko, A.; Hirata, F. Three-dimensional density profiles of water in contact with a solute of arbitrary shape: A RISM approach. Chem. Phys. Lett. 1998, 290, 237−244. (40) Kovalenko, A.; Ten-no, S.; Hirata, F. Solution of threedimensional reference interaction site model and hypernetted chain equations for simple point charge water by modified method of direct inversion in iterative subspace. J. Comput. Chem. 1999, 20, 928−936. (41) Luchko, T.; Gusarov, S.; Roe, D. R.; Simmerling, C.; Case, D. A.; Tuszynski, J.; Kovalenko, A. Three-dimensional molecular theory of solvation coupled with molecular dynamics in Amber. J. Chem. Theory Comput. 2010, 6, 607−624. (42) Abrashkin, A.; Andelman, D.; Orland, H. Dipolar PoissonBoltzmann equation: ions and dipoles close to charge interfaces. Phys. Rev. Lett. 2007, 99, 077801. (43) Azuara, C.; Orland, H.; Bon, M.; Koehl, P.; Delarue, M. Incorporating dipolar solvents with variable density in PoissonBoltzmann electrostatics. Biophys. J. 2008, 95, 5587−5605. (44) Guo, Z.; Li, B.; Dzubiella, J.; Cheng, L.-T.; McCammon, J. A.; Che, J. Evaluation of hydration free energy by level-set variational implicit-solvent model with Coulomb-field approximation. J. Chem. Theory Comput. 2013, 9, 1778−1787. (45) Zhao, Y.; Kwan, Y.; Che, J.; Li, B.; McCammon, J. A. Phase-field approach to implicit solvation of biomolecules with Coulomb field approximation. J. Chem. Phys. 2013, 139, 024111. (46) Fennell, C. J.; Dill, K. A. Physical modeling of aqueous solvation. J. Stat. Phys. 2011, 145, 209. (47) Kodama, R.; Roth, R.; Harano, Y.; Kinoshita, M. Morphometric approach to thermodynamic quantities of solvation of complex molecules: extension to multicomponent solvent. J. Chem. Phys. 2011, 135, 045103. (48) Harano, Y.; Roth, R.; Chiba, S. A morphometric approach for the accurate solvation thermodynamics of proteins and ligands. J. Comput. Chem. 2013, 34, 1969−1974. (49) Remsing, R. C.; Weeks, J. D. The Role of local response in ion solvation: Born theory and beyond. J. Phys. Chem. B 2016, 120, 6238− 6249. (50) Remsing, R. C.; Baer, M. D.; Schenter, G. K.; Mundy, C. J.; Weeks, J. D. The role of broken symmetry in solvation of a spherical Q


Article

Journal of Chemical Theory and Computation (70) Baker, N. A.; Sept, D.; Holst, M. J.; McCammon, J. A.; Joseph, S. Electrostatics of Nanoysystems: Application to microtubules and the ribosome. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10037−10041. (71) Yu, S. N.; Zhou, Y. C.; Wei, G. W. Matched interface and boundary (MIB) method for elliptic problems with sharp-edged interfaces. J. Comput. Phys. 2007, 224, 729−756. (72) Cai, Q.; Hsieh, M.-J.; Wang, J.; Luo, R. Performance of nonlinear finite-diference Poisson-Boltzmann solvers. J. Chem. Theory Comput. 2010, 6, 203−211. (73) Jackson, J. D. Classical Electrodynamics, 3rd ed.; Wiley: New York, 1998. (74) Mobley, D. L.; Dill, K. A.; Chodera, J. D. Treating entropy and conformational changes in implicit solvent simulations of small molecules. J. Phys. Chem. B 2008, 112, 938−946. (75) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (76) Lu, B.; Cheng, X.; Hou, T.; McCammon, J. A. Calculation of the Maxwell stress tensor and the Poisson-Boltzmann force on a solvated molecular surface using hypersingular boundary integrals. J. Chem. Phys. 2005, 123, 084904. (77) Ganesh, M.; Steinbach, O. The numerical solution of a nonlinear hypersingular boundary integral equation. J. Comput. Appl. Math. 2001, 131, 267−280. (78) Hsiao, G. C.; Wendland, W. L. Encyclopedia of Computational Mechanics; John Wiley & Sons, Ltd, 2004. (79) Feistauer, M.; Hsiao, G. C.; Kleinman, R. E. Asymptotic and a Posteriori Error Estimates for Boundary Element Solutions of Hypersingular Integral Equations. SIAM J. Numer. Anal. 1996, 33, 666−685. (80) McLean, W. Strongly Elliptic Systems and Boundary Integral Equations; Cambridge University Press: New York, 2000. (81) Bardhan, J. P.; Eisenberg, R. S.; Gillespie, D. Discretization of the Induced-Charge Boundary Integral Equation. Phys. Rev. E 2009, 80, 011906 DOI: 10.1103/PhysRevE.80.011906. (82) Berti, C.; Gillespie, D.; Bardhan, J. P.; Eisenberg, R. S.; Fiegna, C. Comparison of three-dimensional Poisson solution methods for particle-based simulation and inhomogeneous dielectrics. Phys. Rev. E 2012, 86, 011912. (83) Saad, Y.; Schultz, M. GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM J. Sci. Stat. Comput. 1986, 7, 856−869. (84) Altman, M. D.; Bardhan, J. P.; Tidor, B.; White, J. K. FFTSVD: A Fast Multiscale Boundary-Element Method Solver Suitable for BioMEMS and Biomolecule Simulation. IEEE T. Comput.-Aid. D 2006, 25, 274−284. (85) Bardhan, J. P.; Hildebrandt, A. A Fast Solver for Nonlocal Electrostatic Theory in Biomolecular Science and Engineering. Proceedings of the 48th Design Automation Conference (DAC). New York, 2011; pp 801−805, San Diego, CA, 5−9 June, 2011. (86) Cooper, C. D.; Bardhan, J. P.; Barba, L. A. A biomolecular electrostatics solver using Python, GPUs and boundary elements that can handle solvent-filled cavities and Stern layers. Comput. Phys. Commun. 2013, 185, 720−729. (87) Barnes, J.; Hut, P. A hierarchical O(N log N) force-calculation algorithm. Nature 1986, 324, 446−449. (88) Connolly, M. L. Analytical Molecular Surface Calculation. J. Appl. Crystallogr. 1983, 16, 548−558. (89) Sanner, M.; Olson, A. J.; Spehner, J. C. Reduced Surface: An Efficient Way to Compute Molecular Surfaces. Biopolymers 1996, 38, 305−320. (90) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 1983, 4, 187−217. (91) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.,

III; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (92) Lamoureux, G.; Roux, B. Absolution Hydration Free Energy scale for alkali and halide ions established from simulations with a polarizable force field. J. Phys. Chem. B 2006, 110, 3308−3322. (93) Bardhan, J. P.; Altman, M. D.; White, J. K.; Tidor, B.; Willis, D. J.; Lippow, S. M. Numerical Integration Techniques for Curved-Panel Discretizations of Molecule-Solvent Interfaces. J. Chem. Phys. 2007, 127, 014701. (94) Nina, M.; Beglov, D.; Roux, B. Atomic radii for continuum electrostatics calculations based on molecular dynamics free energy simulations. J. Phys. Chem. B 1997, 101, 5239−5248. (95) Mobley, D. L.; Bayly, C. I.; Cooper, M. D.; Shirts, M. R.; Dill, K. A. Small molecule hydration free energies in explicit solvent: an extensive test of fixed-charge atomistic simulations. J. Chem. Theory Comput. 2009, 5, 350−358. (96) Mobley, D. L.; Bayly, C. I.; Cooper, M. D.; Dill, K. A. Predictions of hydration free energies from all-atom molecular dynamics simulations. J. Phys. Chem. B 2009, 113, 4533. (97) Mobley, D. L.; Guthrie, J. P. FreeSolv:a database of experimenetal and calculated hydration free energies, with input files. J. Comput.-Aided Mol. Des. 2014, 28, 711−720. (98) Helland, R.; Otlewski, J.; Sundheim, O.; Dadlez, M.; Smalas, A. O. The crystal structures of the complexes between bovine betatrypsin and ten P1 variants of BPTI. J. Mol. Biol. 1999, 287, 923−942. (99) Kreienkamp, A.; Liu, L. Y.; Minkara, M. S.; Knepley, M. G.; Bardhan, J. P.; Radhakrishnan, M. L. Analysis of fast boundary-integral approximations for modeling electrostatic contributions of molecular binding. Mol. Based Math. Biol. 2013, 1, 124−150. (100) Bardhan, J. P. Nonlocal continuum electrostatic theory predicts surprisingly small energetic penalties for charge burial in proteins. J. Chem. Phys. 2011, 135, 104113. (101) Marenich, A. V.; Kelly, C. P.; Thompson, J. D.; Hawkins, G. D.; Chambers, C. C.; Giesen, D. J.; Winget, P.; Cramer, C. J.; Truhlar, D. G. Minnesota Solvation Database, version 2012; University of Minnesota, 2012 (102) Sitkoff, D.; Sharp, K. A.; Honig, B. Accurate Calculation of Hydration Free Energies Using Macroscopic Solvent Models. J. Phys. Chem. 1994, 98, 1978−1988. (103) Park, I.; Jang, Y. H.; Hwang, S.; Chung, D. S. PoissonBoltzmann continuum solvation models for nonaqueous solvents. I. 1octanol. Chem. Lett. 2003, 32, 376−377. (104) Zhao, Y. H.; Abraham, M. H. Octanol/water partition of ionic species, including 544 cations. J. Org. Chem. 2005, 70, 2633−2640. (105) Baer, M. D.; Stern, A. C.; Levin, Y.; Tobias, D. J.; Mundy, C. J. Electrochemical surface potential due to classical point charge models drives anion adsorption to the air-water interface. J. Phys. Chem. Lett. 2012, 3, 1565−1570. (106) Kharkats, Y. I.; Kornyshev, A. A.; Vorotyntsev, M. A. Electrostatic models in the theory of solutions. J. Chem. Soc., Faraday Trans. 2 1976, 72, 361−371. (107) Luo, R.; Moult, J.; Gilson, M. K. Dielectric screening treatment of electrostatic solvation. J. Phys. Chem. B 1997, 101, 11226−11236. (108) Atkinson, K. E. The Numerical Solution of Integral Equations of the Second Kind; Cambridge University Press: New York, 1997. (109) Juffer, A. H.; Botta, E. F. F.; van Keulen, B. A. M.; van der Ploeg, A.; Berendsen, H. J. C. The Electric Potential of a Macromolecule in a Solvent: A Fundamental Approach. J. Comput. Phys. 1991, 97, 144−171. (110) Bordner, A. J.; Huber, G. A. Boundary element solution of the linear Poisson-Boltzmann equation and a multipole method for the rapid calculation of forces on macromolecules in solution. J. Comput. Chem. 2003, 24, 353−367. (111) Andjelic, Z.; Of, G.; Steinbach, O.; Urthaler, P. Boundary element methods for magnetostatic field problems: a critical view. Comput. Visualization Sci. 2011, 14, 117−130.

R


Extending the Solvation-Layer Interface Condition Continum

Recommend Documents