PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Commun Comput Phys. Author manuscript; available in PMC 2010 June 15.
Published in final edited form as:
Commun Comput Phys. 2009 November; 6(5): 955–977.
doi:  10.4208/cicp.2009.v6.p955
PMCID: PMC2885744
NIHMSID: NIHMS201932

Image Charge Methods for a Three-Dielectric-Layer Hybrid Solvation Model of Biomolecules

Abstract

This paper introduces a three dielectric layer hybrid solvation model for treating electrostatic interactions of biomolecules in solvents using the Poisson-Boltzmann equation. In this model, an interior spherical cavity will contain the solute and some explicit solvent molecules, and an intermediate buffer layer and an exterior layer contain the bulk solvent. A special dielectric permittivity profile is used to achieve a continuous dielectric transition from the interior cavity to the exterior layer. The selection of this special profile using a harmonic interpolation allows an analytical solution of the model by generalizing the classical Kirkwood series expansion. Discrete image charges are used to speed up calculations for the electrostatic potential within the interior and buffer layer regions. Semi-analytical and least squares methods are used to construct an accurate discrete image approximation for the reaction field due to solvent with or without salt effects. In particular, the image charges obtained by the least squares method provide accurate approximations to the reaction field independent of the ionic concentration of the solvent. Numerical results are presented to validate the accuracy and effectiveness of the image charge methods.

Keywords: Image charge methods, distance-dependent dielectric permittivity, hybrid implicit/explicit model, reaction field, Poisson-Boltzmann equation, protein

1 Introduction

The study of electrostatic interactions using macroscopic continuum models [5, 11, 17, 26, 40] has been widely conducted for investigating structure, function, and properties of protein molecules in an aqueous environment. These implicit models characterize the solvent in terms of macroscopic physical quantities such as dielectric constants and Debye lengths, and thus greatly reduce the degrees of freedom in comparison with explicit atomistic solvent models. In macroscopic models, the solute molecule is usually considered as a uniform low-dielectric medium (with a dielectric constant between one and four) with a fixed charge distribution, and the solvent is treated as a homogeneous medium with a high dielectric constant, such as 80 for water. A Linearized Poisson-Boltzmann (LPB) equation is then solved to obtain the electrostatic potential of the system.

Historically, Born [8] first studied the solvation effects for an ion placed at the center of a spherical region with a low dielectric constant embedded in a high-dielectric medium, and derived the electrostatic free energy. Onsager [32] extended this study to a dipole. Both Born and Onsager models are special cases of the results of Kirkwood [25] and Tan-ford and Kirkwood [38], which represented the solute molecule as a collection of fixed charges within a spherical cavity. The LPB equation was solved analytically for spherical geometries in these early studies. More recent work has considered ellipsoids [4,13]. For an irregular boundary, numerical methods [6,24,28] such as finite difference and finite element methods in three-dimensional grids must be used. However, numerically solving the LPB equation is computationally expensive. To analytically treat biomolecules with general geometries, the hybrid implicit/explicit model using a spherical cavity has been developed [7, 31]. In addition to the biomolecule solute, such as a protein, the spherical cavity contains several layers of explicit water molecules to model the interactions between the solute and solvent molecules.

In previous work, we developed a multiple image charge approximation [9] to the Kirkwood series solution using numerical quadratures of the line image charge representation [18, 27, 29, 30] of the reaction field from pure water. The locations of resulting image charges are related to Gauss quadrature points. Although less accurate, a single image charge [19] and its improved version [1] have been widely applied in molecular dynamic and Monte Carlo simulations [21,33,41]. Employing more images improves the approximation of the Kirkwood solution, especially closer to the boundary. The method of multiple discrete images has also been extended to treat the reaction field for ionic solvents [14,15,43] in the case of low ionic strength.

Unfortunately, continuum solvation models with piecewise constant dielectric functions produce unphysical reaction fields within the spherical cavity near the boundary. This artifact strongly affects the charged solvent molecules (i.e. water) near the boundary. Inaccurate estimation of the pairwise electrostatic interaction [35] results because of improper electrostatic screening by the high-dielectric solvent. To overcome this drawback, pairwise electrostatic interactions near the interface have been modeled using a distance-dependent effective dielectric constant differing from the two homogeneous phases [16]. The simplest of these models assigns an effective dielectric constant between all pairs of charged particles as a linear function [42] of their separation distance, but this creates an inconsistency with the bulk dielectric constant for long-range interactions. Improvements along these lines are obtained by imposing sigmoidal forms [22,23,34,36] of dielectric functions, which approximate the low dielectric constant for short-range interactions and the bulk dielectric constant of solvent for the long-range interactions. These models attempt to compensate for the source of inaccuracies, which derive from the sharp interface between two dielectric constants used in the spherical hybrid solvation model [7,31].

The sharp jump in the dielectric permittivity induces a singular reaction field when an explicit charge within the cavity approaches the interface, which causes the self contribution of the electrostatic solvation energy to diverge. Simply increasing the size of the spherical cavity and the amount of explicit water within does not eliminate unphysical long range effects that propagate throughout the cavity due to the influence of the boundary. Details of how to minimize these unwanted effects within molecular dynamic simulations will be published elsewhere. However, the unphysical boundary effects can be reduced using a three-layer model [12, 35], which was originally developed by using an intermediate dielectric constant in a buffer layer. This idea can be further extended to eliminate any artificial discontinuity. In the present work, a smooth transition layer from low to high dielectric mediums is proposed to form a new three-layer dielectric model. The first inner layer defines a spherical cavity with a low dielectric constant, containing the biomolecule under study and explicit solvent molecules. The second intermediate layer defines a spherical shell as a buffer zone of solvent described by a dielectric permittivity with a radial distance dependent function. The third layer consists of the bulk solvent characterized by a high dielectric constant.

The electrostatic potential can be solved everywhere within the first two layers of the three-layer model by generalizing the Kirkwood series expansion [25] in terms of Legendre polynomials. From a computational point of view, direct summation of this series expansion is very expensive, more than the Kirkwood series [3,39]. However, this computational cost can be by-passed by using a multiple image charge method for the three-layer model. In this paper, we will present a fast and accurate method to calculate the reaction field for a single spherical cavity having a boundary layer to diminish the un-physical artifacts due to the abrupt transition between dielectric constants inside and outside the cavity. We will show that the complexity of solving this three-layer hybrid model using the method of multiple image charges is no greater than for the hybrid model with discontinuous change in dielectric at the boundary of the sphere. Specifically, compared with using exact expansions, the multiple image charge method is computationally much less expensive when used together with the order N fast multipole method [9,10].

To proceed in this direction, the first step is to develop the formulas for the image charges. Herein, two multiple image approaches are developed for the proposed three-layer model. First, a semi-analytical approach extends the approach of using numerical quadratures to line images, giving an approximation of the reaction field in an order An external file that holds a picture, illustration, etc.
Object name is nihms201932ig1.jpg(h2/a2) where h is the thickness of the buffer layer and a is the inner spherical radius. Second, a least squares approach minimizing the error of the reaction field is used to optimize image charges at pre-defined image locations, which is shown to produce more accurate results than the semi-analytical approach. Specifically, the image charges obtained by the least squares method provide accurate approximations to the reaction field independent of the ionic concentration in the solvent. In contrast, the accuracy of our previously obtained image charges for the reaction field depends on the solvent’s ionic strength [14,15,43].

The rest of the paper is organized as follows: In Section 2, we introduce the three-layer hybrid solvation model for electrostatics of biomolecules in an aqueous environment, and discuss the series solution in terms of Legendre polynomials. In Sections 3 and 4, semi-analytical and least squares methods for finding multiple image charge approximations are studied. In Section 5, extension to solvents with salt effects is considered. In Section 6, numerical results are given to show the performance of the proposed image methods. In Section 7, conclusions are made.

2 Three-layer solvation model in non-ionic solvents and its series solution

A three-layer hybrid solvation model is considered for the electrostatic interactions of biomolecules in a solvent, which partitions the solute/solvent system into three non-overlapping parts, as shown in Fig. 1. Inside the interior sphere (ra), the solute and solvent molecules are explicitly described by an array of fixed charges, and the dielectric constant εi is assumed to be the free-space permittivity or slightly greater (from 1 to 4). Outside the larger spherical cavity (rb), the bulk solvent is represented by a continuum medium with a high dielectric constant εo. The intermediate layer of implicit solvent between the two spherical surfaces is considered as a buffer zone [3,12] with a thickness h = ba, and its dielectric permittivity ε(r) is assumed to be dependent on the radial distance from the origin. With this setting, the electrostatic potential Φ is given by the following Poisson equation:

·(εiΦ(r))=ρ1(r),ra,
(2.1)

·(ε(r)Φ(r))=0,a<rb,
(2.2)

ΔΦ(r)=0,r>b,
(2.3)

where ρ1(r)= Σj qjδ(|rrj|) is the charge distribution inside the interior cavity. Here rj is the location of charge qj, [nabla] and Δ are gradient and Laplace operators, respectively.

Figure 1
Schematic of a point charge q inside the interior sphere of the three-layer dielectric model.

Without loss of generality, using the linear superposition principle, we consider a single source charge q located on the x-axis at a distance rs from the center of the sphere (see Fig. 1); i.e., in Eqs. (2.1)(2.3), the charge distribution ρ1(r)=(|rrs|), where rs =(rs,0,0) under a spherical coordinate system r= (r,θ,[var phi]). Because of the azimuthal symmetry, the potential Φ(r,θ) in the exterior layer can be expressed as

Φ(r,θ)=n=0Anrn+1Pn(cosθ),forr>b,
(2.4)

where Pn(x) are Legendre polynomials of order n, and the constant coefficients An are to be determined.

On the other hand, for a field point r in the interior layer, i.e., 0≤ra, the potential due to the charge is the potential in free space q/(4πεi R), plus a reaction potential induced by the polarization of the dielectric medium (|r|>a):

Φ(r,θ)=q4πεiR+n=0BnrnPn(cosθ),
(2.5)

where R is the distance between the field and source points,

R=rrs=rs2+r22rrscosθ.

By using the expansion of the reciprocal distance

q4πεiR={q4πεirn=0(rsr)nPn(cosθ),rsra,q4πεirsn=0(rrs)nPn(cosθ),0rrs,
(2.6)

the total potential inside the interior sphere can be rewritten as

Φ(r,θ)={n=0[q4πεir(rsr)n+Bnrn]Pn(cosθ),rsra,n=0[q4πεirs(rrs)n+Bnrn]Pn(cosθ),0rrs.
(2.7)

In the intermediate layer with continuous dielectric constants ε(r), the potential at r is represented in the following form:

Φ(r,θ)=n=0(Cn(r)rn+Dn(r)r(n+1))Pn(cosθ),
(2.8)

where coefficients Cn(r) and Dn(r) are continuous functions of r, which can be determined by two differential equations [37]. To obtain these equations, the buffer layer is decomposed into thin shells (Fig. 2) and at each shell the coefficients Cn, Dn, and dielectric permittivity ε(r) can be approximated by piecewise constants. Consider the l-th and (l+1)-th shells,

Figure 2
The intermediate layer (a<r <b) is divided into several thin shells. The dielectric permittivity ε(r) in each shell can be considered as a constant.
Φl(r,θ)=n=0(Cnlrn+Dnlr(n+1))Pn(cosθ),
(2.9a)
Φl+1(r,θ)=n=0(Cnl+1rn+Dnl+1r(n+1))Pn(cosθ).
(2.9b)

By using the continuities of the potential and the flux normal to the interface

Φl(r+,θ)=Φl+1(r,θ),
(2.10a)

εlΦl(r,θ)r|r=r+=εl+1Φl+1(r,θ)r|r=r,
(2.10b)

and the orthogonality of the Legendre polynomials, we have

Cnlrn+Dnlr(n+1)=Cnl+1+Dnl+1r(n+1),
(2.11a)
εl(nCnlrn1Dnl(n+1)r(n+1))=εl+1(nCnl+1rn1(n+1)Dnl+1r(n+2)).
(2.11b)

For an infinitely thin shell, Cnl+1 can be replaced by Cn(r+dr) and Cnl by Cn(r). The same steps can be taken for Dnl,Dnl+1 and εl, εl+1. As the function f (r+dr) can be expressed as

f(r+dr)=f(r)+f(r)dr,

two coupled differential equations can then be obtained:

(2n+1)ε(r)Cn(r)=nε(r)Cn(r)+(n+1)r(2n+1)ε(r)Dn(r),
(2.12a)

(2n+1)ε(r)Dn(r)=nr2n+1ε(r)Cn(r)(n+1)ε(r)Dn(r).
(2.12b)

By multiplying Eq. (2.12a) by r2n+1, together with Eq. (2.12b), we obtain

Dn(r)=r2n+1Cn(r).
(2.13)

By plugging Eq. (2.13) into Eq. (2.12), two second-order differential equations of Cn(r) and Dn(r) are obtained [37]

rε(r)ε(r)Cn(r)+((2n+1)ε(r)ε(r)+2rε2(r)rε(r)ε(r))Cn(r)+nε2(r)Cn(r)=0,
(2.14a)
rε(r)ε(r)Dn(r)+(2rε2(r)(2n+1)ε(r)ε(r)rε(r)ε(r))Dn(r)(n+1)ε2(r)Dn(r)=0.
(2.14b)

Based on these equations, a linear dielectric profile ε(r) = d1 +d2r is natural for making a continuous transition of the permittivity, where constants d1 and d2 can be determined by interpolating the values at interface r = a and r = b. However, the obtained functions Cn(r) and Dn(r), which can be represented by a Taylor series expansion, are inefficient for computations.

Therefore, we aim to construct a dielectric function in the layer such that all of Cn(r) and Dn(r) are linearly dependent. Note that if ε(r) satisfies Δε(r)=0 in Eq. (2.2), then the potential holds

Δ(ε(r)Φ(r))=0.

Then, by a simple transform, the potential in the layer can be expressed by

Φ(r,θ)=1ε(r)n=0(Cnrn+Dnr(n+1))Pn(cosθ).
(2.15)

To match this result, the interpolation with two harmonic functions 1 and 1r are used to obtain

ε(r)=α+βr,
(2.16)

where α,β can be determined by the continuity of the dielectric permittivity at the interface (Fig. 3), say,

Figure 3
Linear (dot line) and (α+βr)2 (solid line) dielectric profiles in the intermediate layer.
α=aεi+bεoh,β=ab(εiεo)h.

Together with Eqs. (2.5), (2.7) and (2.15), the expansion coefficients An,Bn,Cn and Dn can be determined by the continuities of the potentials and the fluxes normal to the boundaries, i.e.,

Φ(a+,θ)=Φ(a,θ),
(2.17a)
ε(a+)Φ(r,θ)r|r=a+=εiΦ(r,θ)r|r=a,
(2.17b)
Φ(b+,θ)=Φ(b,θ),
(2.17c)
εoΦ(r,θ)r|r=b+=ε(b)Φ(r,θ)r|r=b.
(2.17d)

Using the orthogonality of Legendre polynomials, we can obtain An,Bn,Cn and Dn. We are only interested in the interior layer, where

Bn=T1T2,
(2.18)

and, T1 and T2 are given by

T1=[a1n(εiεo)(a3+2n(εiεo)a3(εiεo)(a+h)2nεoh3(a+h)2n(1+2n)+a2+2nh(εo+2εi(1+n))a2h(a+h)2n(2εi+εo(1+2n))ah2(a+h)2n(εi+εo(1+4n)))]q(rsa)n,T2=4πεi(a+h)(a2+2n(εiεo)2a2(εiεo)2(a+h)2n+2a(εiεo)2h(a+h)2n+εoh2(a+h)2n(1+2n)(εo+2εin)).

In particular, when h=0, it reproduces the Kirkwood’s result [25]:

Bn=a1n(εiεo)(1+n)q(rsa)n4πεi(εo+εin+εon).
(2.19)

3 Semi-analytical multiple image approximation for non-ionic solvents

We now turn to the problem of finding image charges outside the interior sphere to accurately approximate the reaction field expressed by the series expansion of the Legendre polynomials. It will be shown that with only a few image charges, rapid and accurate calculation of the reaction field can be achieved.

3.1 Line image approximation

To obtain a line image, we approximate the reaction field inside the sphere by truncating the Taylor series expansion at An external file that holds a picture, illustration, etc.
Object name is nihms201932ig1.jpg(h2) in terms of the thickness h of the buffer layer as

ΦRF(r,θ)=n=0BnrnPn(cosθ)=q4πεian=0(rri)n{C1+C2n+C3n+1γ2+C4(n+1γ2)2}Pn(cosθ)=O(h2a2),
(3.1)

where

γ=εiεoεi+εo,ri=a2rs

and Ck for k = 1,···,4 are constants depending on a, h and γ:

C1=γh3a(1+γ1γ2)(52γ2+(12+γ)1γ2),C2=h3a(1+γ1γ2)(2γ+1γ2),C3=γ(1γ)2h12a(1+γ1γ2)(1+4γ3γ3+1γ2(3γ21)),C4=h24a(γγ3)(1+γ1γ2)2.

By using formula (2.6), the first term in Eq. (3.1) is the contribution of an image charge at Kelvin image point ri = (ri,0,0), i.e.,

q4πεian=0rnrinPn(cosθ)=qri4πεia1rri.
(3.2)

For the second term in (3.1), we have

n=0(rri)nnPn(cosθ)=n=0(rri)nnPn(cosθ)=rr(n=1(rri)nPn(cosθ))=rr(n=0(rri)nPn(cosθ))=rrir(1rri).
(3.3)

To obtain an image formula of the third term in (3.1), we denote function

f(r,x)=n=0(rx)n1n+1γ2Pn(cosθ)forxr.
(3.4)

Multiplying both sides by x1γ2 and differentiating with respect to x, we have

x(f(r,x)x1γ2)=n=01x1γ2+1(rx)nPn(cosθ)=x1γ21rx.
(3.5)

Taking the integration from ri to ∞, due to f (r,x) → 0 as x → ∞, we transform the series into a line image:

n=0(rri)n1n+1γ2Pn(cosθ)=f(r,ri)=ri(xri)1γ21rxdx.
(3.6)

Similarly, the fourth term in (3.1) can also be represented by a line image as follows. Let

ϕ(r,x)=n=01(n+1γ2)2(rx)nPn(cosθ),forxr,
(3.7)
g(r,x)=x(ϕ(r,x)x1γ2).
(3.8)

Calculating the differentiation in g(r,x), we have

g(r,x)=n=0x(n+1γ2+1)n+1γ2rnPn(cosθ).
(3.9)

Multiplying g by x and taking the derivative with respect to x yields:

x(xg(r,x))=n=0x(n+1+1γ2)rnPn(cosθ)=x1γ21rx.
(3.10)

Again, taking the integration from x to ∞, because of

ϕ(r,x),ϕx(r,x)0asx,

we get

g(r,x)=1xxt1γ21rtdt,
(3.11)

with vector t = (t,0,0). With a simple calculation, we also have

ϕ(r,ri)ri1γ2=ri1xxt1γ21rtdtdx=rit1γ2rtrit1xdxdt=rit1γ2lntrirtdt.
(3.12)

We thus obtain the line image from the Kelvin point ri to infinity:

n=01(n+1γ2)2(rri)nPn(cosθ)=ri(xri)1γ2lnxrirxdx.
(3.13)

Finally, we obtain the image charge approximation to the reaction field, which consists of two Kelvin images plus two line images from the Kelvin point to the infinity, as,

ΦRF(r,θ)=W1qi14πεirriW2qi14πεirr(1rri)+W3riqi2(x)4πεirxdxW4riqi2(x)lnxri4πεirxdx+O(h2a2),
(3.14)

where

qi1=γarsq,qi2=qaγ(1+γ)2(xri)1γ2,γ=εiεoεi+εo,

and

W1=113γ(1+γ1γ2)(52γ2+(γ+12)1γ2)ha,W2=13γ(1+γ1γ2)(2γ+1γ2)ha,W3=11+γ1γ23γ(1+γ)·1+4γ3γ3+(3γ21)1γ22ha,W4=112(1+γ1γ2)(1γ2+(γ1)1γ2)ha.

3.2 Discretization of the line image charges

In order to discretize the line image charges into multiple point charges, we use a numerical quadrature rule to approximate the following integral

I=ri1rx(xri)1γ2dx.
(3.15)

By changing variable rix=(1s2)τ with τ >0, we have,

I=τ2γ12τ11(1s)α·h(r,s;τ)ds,
(3.16)

where

α=(1γ)τ21,h(r,s;τ)=2τri(1s)τr2τri.
(3.17)

We use the Gauss-Radau quadrature [9,15] which is based on Jacobi polynomials Pn(α,β)(s) orthogonal on the interval (−1,1) with the weight function

w(α,β)(s)=(1s)α(1+s)β,

i.e.,

11(1s)α(1+s)βPj(α,β)(s)Pk(α,β)(s)ds=δjk,
(3.18)

where α >−1, β >−1. Take s0 = −1, α = (1 − γ)τ/2 − 1 and β = 0, and denote the Jacobi-Gauss-Radau points and weights by sm,ωm, for m = 0,1,2,···, M, which can be generated by the package ORTHPOL [20]. The integral I is then discretized as

Iτ2γ12τm=0Mωmh(r,sm;τ).
(3.19)

Thus, for the line image charges, we have,

riqi2(x)4πεirxdxm=0Mqm4πεirxm,
(3.20)

riqi2(x)lnxri4πεirxdxm=0Mqmlnxmri4πεirxm,
(3.21)

where

qm=2γ12τ1γ(1+γ)τωmqxma,andxm=ri(21sm)τ.
(3.22)

Therefore, multiple image charges approximation to the reaction field due to the point charge at rs is obtained:

ΦRF(r)W1qi14πεirriW2qi14πεirr(1rri)+m=0M(W3W4lnxmri)qm4πεirxm,

where the derivative of the Kelvin image with respective to r can be expressed as:

rr(1rri)=ri2r22rri312rri.
(3.23)

4 Least squares multiple image approximation

An alternative to the semi-analytical approach is the least squares method to obtain the multiple image approximation. In this procedure, we fix the locations of all image charges as given in (3.22), and their strengths are to be determined by minimizing the least squares error with respect to the exact series solution. The advantage of this method is that the details of the functional form of an integrand does not affect the choice of the Gauss-quadrature points. Only the overall scaling properties of the function are important, but these are captured well knowing the analytical solution that was presented above (i.e. Jacobi-Gauss-Radau quadrature). Moreover, the least squares error method attempts to spread the error uniform throughout the cavity. The question we will answer, and demonstrated by numerical examples (see below) is: How will the least squares solution for the image charges compare to the semi-analytical method?

We first assume that the potential is given as follows:

Φ(r)=q4πεirrs+m=0Mqm4πεirxm.
(4.1)

The least squares method is used to find the charge strengths by minimizing the L2-error of the reaction potential induced by the discrete images,

Err(q)=n=1N(m=0Mqm4πεirnxmΦRF(rn))2,
(4.2)

where rn for n = 1,···, N are N field points inside the interior sphere. Our least squares procedure amounts to finding q0=(q00,,qM0), such that

Err(q0)=min(Err(q)).

Taking the derivative with respective to ql,

Errql=0,forl=0,,M,

we obtain a linear algebraic system with M equations:

m=0M(n=1N1rnxmrnxl)qm4πεin=1NΦRF(rn)rnxl=0,l=0,,M.

The strengths of the image charges can then be determined by solving the linear system for each source charge inside the interior sphere. A table of image charges and their locations for discrete source points can be made for practical simulations, where the image charges at any location rs can be interpolated from the data in the table.

5 Multi-image approximation of reaction fields in ionic solvents

Extension of the three-layer model to solvents with salt effects is based on the following model,

·(εiΦ(r))=ρ1(r),ra,
(5.1)

·(ε(r)Φ(r))λ2εoΦ(r)=0,a<rb,
(5.2)

ΔΦ(r)λ2Φ(r)=0,rb,
(5.3)

where the Poisson equations for bulk solvents have been replaced by the Poisson-Boltzmann equations, and λ is the inverse Debye screening length.

The dielectric function ε(r) = (α+β/r)2 in the buffer layer (arb) is still used as before. Again, only one source charge at location rs = (rs,0,0) for rs < a is considered, i.e., ρ1(r) = (rrs). Due to the azimuthal symmetry, the potential at an observation point r= (r,θ,[var phi]) can be represented by

Φ(r,θ)={n=0[q4πεi1rs(rrs)n+Bnrn]Pn(cosθ),0rrs,n=0[q4πεi1r(rsr)n+Bnrn]Pn(cosθ),rsra,1ε(r)n=0[Cnkn(λεor)+Dnin(λεor)]Pn(cosθ),arb,n=0Ankn(λr)Pn(cosθ),rb,
(5.4)

where in(r) and kn(r) are the first and the second kind of modified spherical Bessel functions [2] defined recursively by

{i0(r)=sinh(r)r,i1(r)=sinh(r)r2+cosh(r)r,in(r)=in2(r)2n1rin1(r),forn>1,
(5.5)

and

{k0(r)=π2rer,k1(r)=π2rer(1+1r),kn(r)=kn2(r)+2n1rkn1(r),forn>1.
(5.6)

The expansion coefficients Bn are determined by the continuity of the potentials and the fluxes normal to the boundaries (2.17a)–(2.17d), which are given by

Bn=qa1(rsa2)n4πεiS1S2,
(5.7)

where

S1=in(t2)(αkn(u)b(α+bβ)λkn(u))s11+((2+n)α+a(1+n)β)in(t1)s12t2(α+bβ)in(t2)kn(u)s11+t1(α+aβ)in(t1)s13,
(5.8)
S2=in(t2)(αkn(u)b(α+bβ)λkn(u))s21+(a+nα+anβ)in(t1)s12+t2(α+bβ)in(t2)kn(u)s21t1(α+aβ)in(t1)s13.
(5.9)

Here,

s11=((2+n)α+a(1+n)β)kn(t1)+t1(α+aβ)kn(t1),
(5.10a)

s12=b(α+bβ)λkn(t2)kn(u)+kn(u)(αkn(t2)+t2(α+bβ)kn(t2)),
(5.10b)

s13=kn(t2)(αkn(u)b(α+bβ)λkn(u))+t2(α+bβ)kn(u)kn(t2),
(5.10c)

s21=(αnαanβ)kn(t1)+t1(α+aβ)kn(t1),
(5.10d)

and

u=λb=λ(a+h),t1=aεoλ,t2=bεoλ.

With the analytical formula, the algorithm of using the least squares error approach to find the multiple image charges is the same as that of the case of pure water solvent. For the ionic solvent, we will show in the numerical examples that the accuracy of the image charge approximation based on the least squares approach does not depend on ionic strength, which was a limitation in our previous image method for ionic solvent using an analytical approach [15,43].

6 Numerical examples

In this section, we will numerically investigate the performance of the multiple image approximations of both semi-analytical and least squares approaches to the reaction field. Unless otherwise specified, the dielectric constants for the interior and exterior layers are set to be εi = 2 and εo = 80. Also, the single unit point charge (q = 1) is located on the x-axis inside the interior sphere at a distance rs < a from the common center of the two spheres. The radius of the interior sphere is assumed to be dimensionless, a = 1. For the least squares approach, we choose field points (r,θ,0) with uniform distribution for r =0, 0.1,···,0.8 and θ=0,π5,,9π5, and the cutoff error for calculating the exact reaction field is chosen to be 10−8, i.e.,

ΦRFexactn=0NBnrnPn(cosθ),

where the (N+1)-term and (N+2)-term are both less than this number.

6.1 Numerical results for pure water solvents

First, we compute the L2-norm relative error of both approaches with the thickness of the buffer zone set as h= ba=0.1 and 0.01, where the error is defined by

||E||2=(i,jΦRF(ri,θj)Φ^RF(ri,θj)2)12(i,jΦRF(ri,θj)2)12
(6.1)

for ΦRF(ri, θj) and [Phi w/ hat]RF(ri, θj) being the reaction fields computed from the series solution and image approximation, respectively, and the points (ri, θj) chosen from [0,0.8] with a step 0.05 and [0,2π) with a step 2π20. The results are illustrated in Figs. 4 and and5.5. It can be seen that for relatively large h even using M+1 = 10 image charges, the semi-analytical approximation has error between 0.1% and 1%, and its performance is improved with the decrease of the thickness h according to the An external file that holds a picture, illustration, etc.
Object name is nihms201932ig1.jpg(h2/a2) truncation error. In contrast, the error using the least squares approach is sharply reduced with the increase of the image number; when M = 9, it reaches a quite small magnitude around 10−5% for both h = 0.1 and 0.01 cases. These results show that the semi-analytical approach is usually acceptable in simulations, but in order to simulate a system requiring high accuracy, the least squares method should be used.

Figure 4
Semi-analytical vs. least squares image charge methods. Accuracy of the reaction field for the source charge located rs = (rs, 0, 0) with discrete image charges number M+1 = 2 and 10 including the Kelvin image. Left: h=0.1; right: h=0.01.
Figure 5
Semi-analytical vs. least squares image charge methods. Left: accuracy with the thickness of the buffer layer h = ba with the source location rs = 0.4, using 2 and 10 discrete image charges; right: accuracy with the number of discrete image ...

As shown in Fig. 5, we see that a small number of discrete image charges is enough to provide a high-accurate approximation to the reaction field, which shows the advantage of using the image method in comparison with directly truncating the series solution. For M =2, the error magnitude of the least squares approach has already reached 0.01%, and for M = 3, it is 0.001%. In Table 1, we list data of the positions and strengths of the image charges for different exterior dielectric constants and source locations. In Fig. 6, we plot the spatial distribution of the relative errors of the least squares image charge approximation along x- and y-axes, where two image charges are used and we can see the maximum error appears near the spherical boundary. Although the absolute deviation between the exact solution from the series expansion and the multiple image charge method based on least squares is always small, the sign of the deviations are sometimes positive or negative. Since the approximate solution is not always over-estimating or always under-estimating the exact solution, there will inevitably be places within the cavity that the two solutions will be exactly equal. When this happens the relative error is zero. The peaks (extra high accuracy) seen in Fig. 6 reflect regions where this accidental zero error situation occurs.

Figure 6
Spatial distribution of the errors by using the least squares approximation with 2 discrete image charges. Left: 21 observation points on the x-axis; right: 21 observation points on the y-axis.
Table 1
Locations and strengths of image charges by the least squares approach with a=1, h=0.1, εi = 2. The first image charge q1 is located at the Kelvin point x1, and the second image charge q2 with location x2.

The proposed harmonic form for ε(r) shown in Fig. 3 allowed finding an analytical solution relatively easy, and its form is close to a linear interpolation. To observe how much difference there is from a model that linearly interpolates, the linear form is approximated using a certain number of uniform steps. The relative L2 errors for both the harmonic and discretized-stepped linear forms with respect to the analytical potential with no buffer zone are shown in Fig. 7. It is found that the relative errors within the inner sphere (r <a) are relatively small with respect to the original no-buffer layer model. In this sense, the exact details of the form of ε(r) seem not to matter much. On the other hand, on the scale of overall relative errors between the various models, it is clear that removing finite jumps in the dielectric constant is influential. More importantly, it is the unphysical characteristics near the boundaries that we are interested in removing. These results suggest that a continuous harmonic function has the advantage to maintain a continuous form for ε(r) within the buffer zone, which will eliminate unwanted artifacts near the boundary caused by discrete jumps, while providing similar results in the interior of the cavity as the no-buffer zone case.

Figure 7
Relative L2 error of various models with different dielectric profiles in the buffer zone compared to the no-buffer zone model. For rs =0.4, the harmonic profile; and linear profile using n piecewise dielectric constants.

6.2 Numerical results for ionic solvents

For the case of solvents with salt effects, we choose field points (r,θ,0) with uniform distribution for r=0,0.05,···,0.8 and θ=0,π10,,19π10 during the least squares procedure. And the cutoff error of the exact reaction fields is chosen as 10−7. The results are illustrated in Figs. 810 where two image charges are used, the inverse Deybe length is fixed to equal 1 (i.e., u = λb = 1), and the buffer thickness is 0.1. Fig. 8 plots the error distribution in space along x- and y-axes and due to the same reason as in pure water, the approximate numerical potential matches the analytical potential at the spikes. Figs. 9 and and1010 show the error versus the source location, the image charge number, the buffer thickness, and the ionic concentration. These results agree with those of the pure water case, and show the attractive accuracy of the least squares approach in approximating the series solution of the reaction field model. In particular, it can be seen in the right picture of Fig. 10 that the error of the multiple image charge approximation remains small for arbitrary ionic concentration, which dramatically improves our previous results on image charge approximations with the semi-analytical approaches [14, 15, 43], for two-layer dielectric model, where the accuracy of the semi-analytical image charge approximation is limited by not only the number of image charges but also by the ionic concentration. As a final note, the L2 errors we show is for a single point source. This error is the most relevant in comparing models because it is independent of system details. Nevertheless, we have not found any accumulated errors when the cavity consists of many source charges. In real applications, if anything, random error distribution works in the favor of the image method.

Figure 8
Spatial distribution of the errors by using the least squares approximation with 2 discrete image charges for the three-layer model. Left: 21 observation points on the x-axis; right: 21 observation points on the y-axis.
Figure 9
Relative error of the least squares approximation compared with the exact solution for the three-layer model. Left: accuracy vs. source charge location rs for different dielectric constants εo of bulk solvents; right: accuracy vs. the number of ...
Figure 10
Error of the least squares approximation. Left: accuracy vs. the thickness of the of the buffer layer using 2 discrete image charges, right: accuracy vs. the ionic strength u.

7 Conclusions

In this paper, we propose a three-layer dielectric solvation model for the electrostatic computation of biomolecules in solvents (with or without salt effects), where the dielectric profile within the buffer zone depends on the radial distance. The analytical solution in series form is obtained by generalizing the classical Kirkwood expansion. Two approaches for finding multiple image charges to approximate the reaction field are developed, and numerical examples are performed to validate the accuracy and effectiveness of the image charge methods. These mathematical investigations are the first step of the applicability of the three-layer model, for which a range of test and calibration, including selecting the optimal thickness of the buffer layer and a comparative study of accuracy and speed performances, are underway. Regardless of what is determined to be the optimal dielectric profile for minimizing artifacts near the boundary of the cavity, the method of multiple image charges will provide a computationally accurate and efficient representation of the reaction field. Moreover, the least squares method allows numerical solutions to be obtained with high accuracy when exact analytical forms are not possible to obtain nor easy to implement numerically.

Acknowledgments

The authors would like to thank the financial support provided by the National Institutes of Health (Grant No. 1R01GM083600-02). Z. Xu is also partially supported by the Charlotte Research Institute through a Duke Postdoctoral Fellowship. W. Cai and Z. Xu are also partially supported by the Department of Energy (Grant No. DEFG0205ER25678).

References

1. Abagyan R, Totrov M. Biased probability Monte Carlo conformational searches and electrostatic calculations for peptides and proteins. J Mol Biol. 1994;235:983–1002. [PubMed]
2. Abramowitz M, Stegun IA. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover; New York: 1964.
3. Alper H, Levy RM. Dielectric and thermodynamic response of a generalized reaction field model for liquid state simulations. J Chem Phys. 1993;99:9847–9852.
4. Ambia-Garrido J, Pettitt BM. Free energy calculations for DNA near surfaces using an ellipsoidal geometry. Commun Comput Phys. 2008;3:1117–1131. [PMC free article] [PubMed]
5. Baker NA. Improving implicit solvent simulations: a Poisson-centric view. Curr Opin Struct Biol. 2005;15:137–143. [PubMed]
6. Baker NA, Holst MJ, Wang F. Adaptive multilevel finite element solution of the Poisson-Boltzmann equation II. Refinement at solvent-accessible surfaces in biomolecular systems. J Comput Chem. 2000;21(15):1343–1352.
7. Belgov D, Roux B. Finite representation of an infinite bulk system - Solvent boundary potential for computer simulations. J Chem Phys. 1994;100:9050–9063.
8. Born M. Volumes and heats of hydration of ions. Z Phys. 1920;1:45–48.
9. Cai W, Deng S, Jacobs D. Extending the fast multipole method to charges inside or outside a dielectric sphere. J Comput Phys. 2007;223:846–864.
10. Cheng H, Greengard L, Roklin V. A fast adaptive multipole algorithm in three dimensions. J Comput Phys. 1999;155:468–498.
11. Cramer CJ, Truhlar DG. Implicit solvation models: Equilibria, structure, spectra and dynamics. Chem Rev. 1999;99:2161–2200. [PubMed]
12. Dai J, Tsukerman I, Rubinstein A, Sherman S. New computational models for electrostatics of macromolecules in solvents. IEEE Trans Magn. 2007;43:1217–1220.
13. Deng S. Electrostatic potential of point charges inside dielectric prolate spheroids. J Electrostat. 2008;66:549–560. [PMC free article] [PubMed]
14. Deng S, Cai W. Discrete image approximations of ionic solvent induced reaction field to charges. Commun Comput Phys. 2007;2:1007–1026.
15. Deng S, Cai W. Extending the fast multipole method for charges inside a dielectric sphere in an ionic solvent: High-order image approximations for reaction fields. J Comput Phys. 2007;227:1246–1266. [PMC free article] [PubMed]
16. Edsall JT, McKenzie HA. Water and protein. II. The location and dynamics of water in proteins and its relation to their stability and properties. Adv Biophys. 1983;16:53–183. [PubMed]
17. Feig M, Brooks CL. III Recent advances in the development and application of implicit solvent models in biomolecule simulations. Curr Opin Struct Biol. 2004;14:217–224. [PubMed]
18. Finkelstein AV. Electrostatic interactions of charged groups in an aqueous medium and their effect on the formation of polypeptide chain secondary structure. Mol Biol. 1977:627–634. [PubMed]
19. Friedman HL. Image approximation to the reaction field. Mol Phys. 1975;29:1533–1543.
20. Gautschi W. Algorithm 726; ORTHPOL- a package of routines for generating orthogonal polynomials and Gauss-type quadrature rules. ACM Trans Math Softw. 1994;20:21–62.
21. Havranek JJ, Harbury PB. Tanford-Kirkwood electrostatics for protein modeling. Proc Natl Acad Sci USA. 1999;96:11145–11150. [PubMed]
22. Hingerty BE, Ritchie RH, Ferrell TL, Turner JE. Dielectric effects in bio-polymers – the theory of ionic saturation revisited. Biopolymers. 1985;24:427–439.
23. Jha AK, Freed KF. Solvation effect on conformations of 1,2:dimethoxyethane: charge-dependent nonlinear response in implicit solvent models. J Chem Phys. 2008;128:034501. [PMC free article] [PubMed]
24. Juffer AH, Botta EFF, Vankeulen BAM, Vanderploeg A, Berendsen HJC. The electric-potential of a macromolecule in a solvent - a fundamental approach. J Comput Phys. 1991;97(1):144–171.
25. Kirkwood JG. Theory of solutions of molecules containing widely separated charges with special applications to awitterions. J Chem Phys. 1934;2:351–361.
26. Koehl P. Electrostatics calculations: latest methodological advances. Curr Opin Struct Biol. 2006;16:142–151. [PubMed]
27. Lindell IV. Electrostatic image theory for the dielectric sphere. Radio Sci. 1992;27:1–8.
28. Lu BZ, Zhou YC, Holst MJ, McCammon JA. Recent progress in numerical methods for the Poisson-Boltzmann equation in biophysical applications. Commun Comput Phys. 2008;3:973–1009.
29. Neumann C. Hydrodynamische untersuchungen: Nebst einem anhange uber die probleme der elektrostatik und der magnetischen induktion. Teubner; Leipzig: 1883. pp. 279–282.
30. Norris WT. Charge images in a dielectric sphere. IEE Proc-Sci Meas Technol. 1995;142:142–150.
31. Okur A, Simmerling C. Hybrid explicit/implicit solvation methods. In: Spellmeyer D, editor. Annu Rep Comput Chem. Vol. 2. 2006.
32. Onsager L. Electric moments of molecules in liquids. J Am Chem Soc. 1936;58:1486.
33. Petraglio G. Nonperiodic boundary conditions for solvated systems. J Chem Phys. 2005;123:044103. [PubMed]
34. Ramstein J, Lavery R. Energetic coupling between DNA bending and base pair opening. P Natl Acad Sci USA. 1988;85:7231–7235. [PubMed]
35. Rubinstein A, Sherman S. Influence of the solvent structure on the electrostatic interactions in proteins. Biophys J. 2004;87:1544–1557. [PubMed]
36. Shen MY, Freed KF. All-atom fast protein folding simulations: The villin headpiece. Proteins. 2002;49:439–445. [PubMed]
37. Sihvola A, Lindell IV. Polarizability and effective permittivity of layered and continuously inhomogeneous dielectric spheres. J Electrom Waves Appl. 1989;3:37–60.
38. Tanford C, Kirkwood JG. Theory of protein titration curves. I. General equations for impenetrable spheres. J Am Chem Soc. 1957;79:5333–5339.
39. Tironi IG, Sperb R, Smith PE, van Gunsteren WF. Generalized reaction field method for molecular dynamics simulations. J Chem Phys. 1995;102:5451–5459.
40. Wang J, Tan C, Tan YH, Lu Q, Luo R. Poisson-Boltzmann solvents in molecular dynamics simulations. Commun Comput Phys. 2008;3:1010–1031.
41. Wang L, Hermans J. Reaction field molecular dynamics simulation with Friedman’s image method. J Phys Chem. 1995;99:12001–12007.
42. Warshel A, Levitt M. Theoretical studies of enzymic reactions: dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J Mol Biol. 1976;103:227–249. [PubMed]
43. Xu Z, Deng S, Cai W. Image charge approximations of reaction fields in solvents with arbitrary ionic strength. J Comput Phys. 2009;228:2092–2099. [PMC free article] [PubMed]