Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2885744

Formats

Article sections

- Abstract
- 1 Introduction
- 2 Three-layer solvation model in non-ionic solvents and its series solution
- 3 Semi-analytical multiple image approximation for non-ionic solvents
- 4 Least squares multiple image approximation
- 5 Multi-image approximation of reaction fields in ionic solvents
- 6 Numerical examples
- 7 Conclusions
- References

Authors

Related links

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.p955PMCID: PMC2885744

NIHMSID: NIHMS201932

Peihua Qin: nc.ca.cc.cesl@niqhp; Zhenli Xu: ude.ccnu@7uxz; Wei Cai: ude.ccnu@iacw; Donald Jacobs: ude.ccnu@1sbocajd

See other articles in PMC that cite the published article.

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.

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 (*h*^{2}/*a*^{2}) 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.

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 (*r* ≤ *a*), 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 (

$$\nabla \xb7({\epsilon}_{i}\nabla \mathrm{\Phi}(\mathbf{r}))=-{\rho}_{1}(\mathbf{r}),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}r\le a,$$

(2.1)

$$\nabla \xb7(\epsilon (r)\nabla \mathrm{\Phi}(\mathbf{r}))=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}a<r\le b,$$

(2.2)

$$\mathrm{\Delta}\mathrm{\Phi}(\mathbf{r})=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}r>b,$$

(2.3)

where *ρ*_{1}(**r**)= Σ* _{j} q_{j}δ*(|

Without loss of generality, using the linear superposition principle, we consider a single source charge *q* located on the *x*-axis at a distance *r _{s}* from the center of the sphere (see Fig. 1); i.e., in Eqs. (2.1)–(2.3), the charge distribution

$$\mathrm{\Phi}(r,\theta )=\sum _{n=0}^{\infty}\frac{{A}_{n}}{{r}^{n+1}}{P}_{n}(cos\theta ),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{for}\phantom{\rule{0.16667em}{0ex}}r>b,$$

(2.4)

where *P _{n}*(

On the other hand, for a field point **r** in the interior layer, i.e., 0≤*r*≤*a*, 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 (|

$$\mathrm{\Phi}(r,\theta )=\frac{q}{4\pi {\epsilon}_{i}R}+\sum _{n=0}^{\infty}{B}_{n}{r}^{n}{P}_{n}(cos\theta ),$$

(2.5)

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

$$R=\phantom{\rule{0.16667em}{0ex}}\mid \mathbf{r}-{\mathbf{r}}_{s}\mid \phantom{\rule{0.16667em}{0ex}}=\sqrt{{r}_{s}^{2}+{r}^{2}-{2rr}_{s}cos\theta}.$$

By using the expansion of the reciprocal distance

$$\frac{q}{4\pi {\epsilon}_{i}R}=\{\begin{array}{ll}\frac{q}{4\pi {\epsilon}_{i}r}\sum _{n=0}^{\infty}{\left(\frac{{r}_{s}}{r}\right)}^{n}{P}_{n}(cos\theta ),\hfill & {r}_{s}\le r\le a,\hfill \\ \frac{q}{4\pi {\epsilon}_{i}{r}_{s}}\sum _{n=0}^{\infty}{\left(\frac{r}{{r}_{s}}\right)}^{n}{P}_{n}(cos\theta ),\hfill & 0\le r\le {r}_{s},\hfill \end{array}$$

(2.6)

the total potential inside the interior sphere can be rewritten as

$$\mathrm{\Phi}(r,\theta )=\{\begin{array}{ll}\sum _{n=0}^{\infty}\left[\frac{q}{4\pi {\epsilon}_{i}r}{\left(\frac{{r}_{s}}{r}\right)}^{n}+{B}_{n}{r}^{n}\right]{P}_{n}(cos\theta ),\hfill & {r}_{s}\le r\le a,\hfill \\ \sum _{n=0}^{\infty}\left[\frac{q}{4\pi {\epsilon}_{i}{r}_{s}}{\left(\frac{r}{{r}_{s}}\right)}^{n}+{B}_{n}{r}^{n}\right]{P}_{n}(cos\theta ),\hfill & 0\le r\le {r}_{s}.\hfill \end{array}$$

(2.7)

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

$$\mathrm{\Phi}(r,\theta )=\sum _{n=0}^{\infty}\left({C}_{n}(r){r}^{n}+{D}_{n}(r){r}^{-(n+1)}\right){P}_{n}(cos\theta ),$$

(2.8)

where coefficients *C _{n}*(

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.

$${\mathrm{\Phi}}^{l}(r,\theta )=\sum _{n=0}^{\infty}\left({C}_{n}^{l}{r}^{n}+{D}_{n}^{l}{r}^{-(n+1)}\right){P}_{n}(cos\theta ),$$

(2.9a)

$${\mathrm{\Phi}}^{l+1}(r,\theta )=\sum _{n=0}^{\infty}\left({C}_{n}^{l+1}{r}^{n}+{D}_{n}^{l+1}{r}^{-(n+1)}\right){P}_{n}(cos\theta ).$$

(2.9b)

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

$${\mathrm{\Phi}}^{l}({r}^{+},\theta )={\mathrm{\Phi}}^{l+1}({r}^{-},\theta ),$$

(2.10a)

$${{\epsilon}^{l}\frac{\partial {\mathrm{\Phi}}_{l}(r,\theta )}{\partial r}|}_{r={r}^{+}}={{\epsilon}^{l+1}\frac{\partial {\mathrm{\Phi}}^{l+1}(r,\theta )}{\partial r}|}_{r={r}^{-}},$$

(2.10b)

and the orthogonality of the Legendre polynomials, we have

$${C}_{n}^{l}{r}^{n}+{D}_{n}^{l}{r}^{-(n+1)}={C}_{n}^{l+1}+{D}_{n}^{l+1}{r}^{-(n+1)},$$

(2.11a)

$${\epsilon}^{l}\left(n{C}_{n}^{l}{r}^{n-1}-{D}_{n}^{l}(n+1){r}^{-(n+1)}\right)={\epsilon}^{l+1}\left(n{C}_{n}^{l+1}{r}^{n-1}-(n+1){D}_{n}^{l+1}{r}^{-(n+2)}\right).$$

(2.11b)

For an infinitely thin shell,
${C}_{n}^{l+1}$ can be replaced by *C _{n}*(

$$f(r+dr)=f(r)+{f}^{\prime}(r)dr,$$

two coupled differential equations can then be obtained:

$$(2n+1)\epsilon (r){C}_{n}^{\prime}(r)=-n{\epsilon}^{\prime}(r){C}_{n}(r)+(n+1){r}^{-(2n+1)}{\epsilon}^{\prime}(r){D}_{n}(r),$$

(2.12a)

$$(2n+1)\epsilon (r){D}_{n}^{\prime}(r)=n{r}^{2n+1}{\epsilon}^{\prime}(r){C}_{n}(r)-(n+1){\epsilon}^{\prime}(r){D}_{n}(r).$$

(2.12b)

By multiplying Eq. (2.12a) by *r*^{2}^{n}^{+1}, together with Eq. (2.12b), we obtain

$${D}_{n}^{\prime}(r)=-{r}^{2n+1}{C}_{n}^{\prime}(r).$$

(2.13)

By plugging Eq. (2.13) into Eq. (2.12), two second-order differential equations of *C _{n}*(

$$r\epsilon (r){\epsilon}^{\prime}(r){C}_{n}^{\u2033}(r)+\left((2n+1)\epsilon (r){\epsilon}^{\prime}(r)+2r{{\epsilon}^{\prime}}^{2}(r)-r\epsilon (r){\epsilon}^{\u2033}(r)\right){C}_{n}^{\prime}(r)+n{{\epsilon}^{\prime}}^{2}(r){C}_{n}(r)=0,$$

(2.14a)

$$r\epsilon (r){\epsilon}^{\prime}(r){D}_{n}^{\u2033}(r)+\left(2r{{\epsilon}^{\prime}}^{2}(r)-(2n+1)\epsilon (r){\epsilon}^{\prime}(r)-r\epsilon (r){\epsilon}^{\u2033}(r)\right){D}_{n}^{\prime}(r)-(n+1){{\epsilon}^{\prime}}^{2}(r){D}_{n}(r)=0.$$

(2.14b)

Based on these equations, a linear dielectric profile *ε*(*r*) = *d*_{1} +*d*_{2}*r* is natural for making a continuous transition of the permittivity, where constants *d*_{1} and *d*_{2} can be determined by interpolating the values at interface *r* = *a* and *r* = *b*. However, the obtained functions *C _{n}*(

Therefore, we aim to construct a dielectric function in the layer such that all of *C _{n}*(

$$\mathrm{\Delta}\left(\sqrt{\epsilon (r)}\mathrm{\Phi}(\mathbf{r})\right)=0.$$

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

$$\mathrm{\Phi}(r,\theta )=\frac{1}{\sqrt{\epsilon (r)}}\sum _{n=0}^{\infty}\left({C}_{n}{r}^{n}+{D}_{n}{r}^{-(n+1)}\right){P}_{n}(cos\theta ).$$

(2.15)

To match this result, the interpolation with two harmonic functions 1 and ${\scriptstyle \frac{1}{r}}$ are used to obtain

$$\sqrt{\epsilon (r)}=\alpha +\frac{\beta}{r},$$

(2.16)

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

Linear (dot line) and
${(\alpha +{\scriptstyle \frac{\beta}{r}})}^{2}$ (solid line) dielectric profiles in the intermediate layer.

$$\begin{array}{l}\alpha =\frac{-a\sqrt{{\epsilon}_{i}}+b\sqrt{{\epsilon}_{o}}}{h},\\ \beta =\frac{ab(\sqrt{{\epsilon}_{i}}-\sqrt{{\epsilon}_{o}})}{h}.\end{array}$$

Together with Eqs. (2.5), (2.7) and (2.15), the expansion coefficients *A _{n}*,

$$\mathrm{\Phi}({a}^{+},\theta )=\mathrm{\Phi}({a}^{-},\theta ),$$

(2.17a)

$${\epsilon ({a}^{+})\frac{\partial \mathrm{\Phi}(r,\theta )}{\partial r}|}_{r={a}^{+}}={{\epsilon}_{i}\frac{\partial \mathrm{\Phi}(r,\theta )}{\partial r}|}_{r={a}^{-}},$$

(2.17b)

$$\mathrm{\Phi}({b}^{+},\theta )=\mathrm{\Phi}({b}^{-},\theta ),$$

(2.17c)

$${{\epsilon}_{o}\frac{\partial \mathrm{\Phi}(r,\theta )}{\partial r}|}_{r={b}^{+}}={\epsilon ({b}^{-})\frac{\partial \mathrm{\Phi}(r,\theta )}{\partial r}|}_{r={b}^{-}}.$$

(2.17d)

Using the orthogonality of Legendre polynomials, we can obtain *A _{n}*,

$${B}_{n}=\frac{{T}_{1}}{{T}_{2}},$$

(2.18)

and, *T*_{1} and *T*_{2} are given by

$${T}_{1}=-\left[{a}^{-1-n}(\sqrt{{\epsilon}_{i}}-\sqrt{{\epsilon}_{o}})\left({a}^{3+2n}(\sqrt{{\epsilon}_{i}}-\sqrt{{\epsilon}_{o}})-{a}^{3}(\sqrt{{\epsilon}_{i}}-\sqrt{{\epsilon}_{o}}){(a+h)}^{2n}-\sqrt{{\epsilon}_{o}}{h}^{3}{(a+h)}^{2n}(1+2n)+{a}^{2+2n}h(-\sqrt{{\epsilon}_{o}}+2\sqrt{{\epsilon}_{i}}(1+n))-{a}^{2}h{(a+h)}^{2n}(2\sqrt{{\epsilon}_{i}}+\sqrt{{\epsilon}_{o}}(-1+2n))-a{h}^{2}{(a+h)}^{2n}(\sqrt{{\epsilon}_{i}}+\sqrt{{\epsilon}_{o}}(1+4n))\right)\right]q{\left(\frac{{r}_{s}}{a}\right)}^{n},{T}_{2}=4\pi {\epsilon}_{i}(a+h)\left({a}^{2+2n}{(\sqrt{{\epsilon}_{i}}-\sqrt{{\epsilon}_{o}})}^{2}-{a}^{2}{(\sqrt{{\epsilon}_{i}}-\sqrt{{\epsilon}_{o}})}^{2}{(a+h)}^{2n}+2a{(\sqrt{{\epsilon}_{i}}-\sqrt{{\epsilon}_{o}})}^{2}h{(a+h)}^{2n}+\sqrt{{\epsilon}_{o}}{h}^{2}{(a+h)}^{2n}(1+2n)(\sqrt{{\epsilon}_{o}}+2\sqrt{{\epsilon}_{i}}n)\right).$$

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

$${B}_{n}=\frac{{a}^{-1-n}({\epsilon}_{i}-{\epsilon}_{o})(1+n)q{({\scriptstyle \frac{{r}_{s}}{a}})}^{n}}{4\pi {\epsilon}_{i}({\epsilon}_{o}+{\epsilon}_{i}n+{\epsilon}_{o}n)}.$$

(2.19)

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.

To obtain a line image, we approximate the reaction field inside the sphere by truncating the Taylor series expansion at (*h*^{2}) in terms of the thickness *h* of the buffer layer as

$${\mathrm{\Phi}}_{\text{RF}}(r,\theta )=\sum _{n=0}^{\infty}{B}_{n}{r}^{n}{P}_{n}(cos\theta )=\frac{q}{4\pi {\epsilon}_{i}a}\sum _{n=0}^{\infty}{\left(\frac{r}{{r}_{i}}\right)}^{n}\left\{{C}_{1}+{C}_{2}n+\frac{{C}_{3}}{n+{\scriptstyle \frac{1-\gamma}{2}}}+\frac{{C}_{4}}{{(n+{\scriptstyle \frac{1-\gamma}{2}})}^{2}}\right\}{P}_{n}(cos\theta )=\mathcal{O}\left(\frac{{h}^{2}}{{a}^{2}}\right),$$

(3.1)

where

$$\gamma =\frac{{\epsilon}_{i}-{\epsilon}_{o}}{{\epsilon}_{i}+{\epsilon}_{o}},\phantom{\rule{0.38889em}{0ex}}{r}_{i}=\frac{{a}^{2}}{{r}^{s}}$$

and *C _{k}* for

$$\begin{array}{l}{C}_{1}=\gamma -\frac{h}{3a}\left(1+\gamma -\sqrt{1-{\gamma}^{2}}\right)\phantom{\rule{0.16667em}{0ex}}\left(\frac{5}{2}-{\gamma}^{2}+(\frac{1}{2}+\gamma )\sqrt{1-{\gamma}^{2}}\right),\\ {C}_{2}=-\frac{h}{3a}\left(1+\gamma -\sqrt{1-{\gamma}^{2}}\right)\left(2-\gamma +\sqrt{1-{\gamma}^{2}}\right),\\ {C}_{3}=\frac{\gamma (1-\gamma )}{2}-\frac{h}{12a}\left(1+\gamma -\sqrt{1-{\gamma}^{2}}\right)\phantom{\rule{0.16667em}{0ex}}\left(1+4\gamma -3{\gamma}^{3}+\sqrt{1-{\gamma}^{2}}(3{\gamma}^{2}-1)\right),\\ {C}_{4}=-\frac{h}{24a}(\gamma -{\gamma}^{3})\phantom{\rule{0.16667em}{0ex}}{\left(1+\gamma -\sqrt{1-{\gamma}^{2}}\right)}^{2}.\end{array}$$

By using formula (2.6), the first term in Eq. (3.1) is the contribution of an image charge at Kelvin image point **r*** _{i}* = (

$$\frac{q}{4\pi {\epsilon}_{i}a}\sum _{n=0}^{\infty}\frac{{r}^{n}}{{r}_{i}^{n}}{P}_{n}(cos\theta )=\frac{{qr}^{i}}{4\pi {\epsilon}_{i}a}\frac{1}{\mid \mathbf{r}-{\mathbf{r}}_{i}\mid}.$$

(3.2)

For the second term in (3.1), we have

$$\sum _{n=0}^{\infty}{\left(\frac{r}{{r}_{i}}\right)}^{n}n{P}_{n}(cos\theta )=\sum _{n=0}^{\infty}{\left(\frac{r}{{r}_{i}}\right)}^{n}n{P}_{n}(cos\theta )=r\frac{\partial}{\partial r}\left(\sum _{n=1}^{\infty}{\left(\frac{r}{{r}_{i}}\right)}^{n}{P}_{n}(cos\theta )\right)=r\frac{\partial}{\partial r}\left(\sum _{n=0}^{\infty}{\left(\frac{r}{{r}_{i}}\right)}^{n}{P}_{n}(cos\theta )\right)=r{r}_{i}\frac{\partial}{\partial r}\left(\frac{1}{\mid \mathbf{r}-{\mathbf{r}}_{i}\mid}\right).$$

(3.3)

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

$$f(r,x)=\sum _{n=0}^{\infty}{\left(\frac{r}{x}\right)}^{n}\frac{1}{n+{\scriptstyle \frac{1-\gamma}{2}}}{P}_{n}(cos\theta )\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{for}\phantom{\rule{0.16667em}{0ex}}x\ge r.$$

(3.4)

Multiplying both sides by
${x}^{-{\scriptstyle \frac{1-\gamma}{2}}}$ and differentiating with respect to *x*, we have

$$\frac{\partial}{\partial x}\left(f(r,x){x}^{-{\scriptstyle \frac{1-\gamma}{2}}}\right)=-\sum _{n=0}^{\infty}\frac{1}{{x}^{{\scriptstyle \frac{1-\gamma}{2}}+1}}{\left(\frac{r}{x}\right)}^{n}{P}_{n}(cos\theta )=-{x}^{-{\scriptstyle \frac{1-\gamma}{2}}}\frac{1}{\mid \mathbf{r}-\mathbf{x}\mid}.$$

(3.5)

Taking the integration from *r _{i}* to ∞, due to

$$\sum _{n=0}^{\infty}{\left(\frac{r}{{r}_{i}}\right)}^{n}\frac{1}{n+{\scriptstyle \frac{1-\gamma}{2}}}{P}_{n}(cos\theta )=f(r,{r}_{i})={\int}_{{r}_{i}}^{\infty}{\left(\frac{x}{{r}_{i}}\right)}^{-{\scriptstyle \frac{1-\gamma}{2}}}\frac{1}{\mid \mathbf{r}-\mathbf{x}\mid}dx.$$

(3.6)

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

$$\varphi (r,x)=\sum _{n=0}^{\infty}\frac{1}{{\left(n+{\scriptstyle \frac{1-\gamma}{2}}\right)}^{2}}{\left(\frac{r}{x}\right)}^{n}{P}_{n}(cos\theta ),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{for}\phantom{\rule{0.16667em}{0ex}}x\ge r,$$

(3.7)

$$g(r,x)=\frac{\partial}{\partial x}\left(\varphi (r,x){x}^{-{\scriptstyle \frac{1-\gamma}{2}}}\right).$$

(3.8)

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

$$g(r,x)=-\sum _{n=0}^{\infty}\frac{{x}^{-(n+{\scriptstyle \frac{1-\gamma}{2}}+1)}}{n+{\scriptstyle \frac{1-\gamma}{2}}}{r}^{n}{P}_{n}(cos\theta ).$$

(3.9)

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

$$\frac{\partial}{\partial x}(xg(r,x))=\sum _{n=0}^{\infty}{x}^{-(n+1+{\scriptstyle \frac{1-\gamma}{2}})}{r}^{n}{P}_{n}(cos\theta )={x}^{-{\scriptstyle \frac{1-\gamma}{2}}}\frac{1}{\mid \mathbf{r}-\mathbf{x}\mid}.$$

(3.10)

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

$$\varphi (r,x),{\varphi}_{x}^{\prime}(r,x)\to 0\phantom{\rule{0.38889em}{0ex}}\text{as}\phantom{\rule{0.16667em}{0ex}}x\to \infty ,$$

we get

$$g(r,x)=-\frac{1}{x}{\int}_{x}^{\infty}{t}^{-{\scriptstyle \frac{1-\gamma}{2}}}\frac{1}{\mid \mathbf{r}-\mathbf{t}\mid}dt,$$

(3.11)

with vector ** t** = (

$$\begin{array}{l}\varphi (r,{r}_{i}){r}_{i}^{-{\scriptstyle \frac{1-\gamma}{2}}}={\int}_{{r}_{i}}^{\infty}\frac{1}{x}{\int}_{x}^{\infty}{t}^{-{\scriptstyle \frac{1-\gamma}{2}}}\frac{1}{\mid \mathbf{r}-\mathbf{t}\mid}\mathit{dtdx}\\ ={\int}_{{r}_{i}}^{\infty}\frac{{t}^{-{\scriptstyle \frac{1-\gamma}{2}}}}{\mid \mathbf{r}-\mathbf{t}\mid}{\int}_{{r}_{i}}^{t}\frac{1}{x}\mathit{dxdt}={\int}_{{r}_{i}}^{\infty}\frac{{t}^{-{\scriptstyle \frac{1-\gamma}{2}}}ln{\scriptstyle \frac{t}{{r}_{i}}}}{\mid \mathbf{r}-\mathbf{t}\mid}dt.\end{array}$$

(3.12)

We thus obtain the line image from the Kelvin point *r _{i}* to infinity:

$$\sum _{n=0}^{\infty}\frac{1}{{\left(n+{\scriptstyle \frac{1-\gamma}{2}}\right)}^{2}}{\left(\frac{r}{{r}_{i}}\right)}^{n}{P}_{n}(cos\theta )={\int}_{{r}_{i}}^{\infty}\frac{{\left({\scriptstyle \frac{x}{{r}_{i}}}\right)}^{-{\scriptstyle \frac{1-\gamma}{2}}}ln{\scriptstyle \frac{x}{{r}_{i}}}}{\mid \mathbf{r}-\mathbf{x}\mid}dx.$$

(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,

$${\mathrm{\Phi}}_{RF}(r,\theta )={W}_{1}\frac{{q}_{i1}}{4\pi {\epsilon}_{i}\mid \mathbf{r}-{\mathbf{r}}_{i}\mid}-{W}_{2}\frac{{q}_{i1}}{4\pi {\epsilon}_{i}}r\frac{\partial}{\partial r}\left(\frac{1}{\mid \mathbf{r}-{\mathbf{r}}_{i}\mid}\right)+{W}_{3}{\int}_{{r}_{i}}^{\infty}\frac{{q}_{i2}(x)}{4\pi {\epsilon}_{i}\mid \mathbf{r}-\mathbf{x}\mid}dx-{W}_{4}{\int}_{{r}_{i}}^{\infty}\frac{{q}_{i2}(x)ln{\scriptstyle \frac{x}{{r}_{i}}}}{4\pi {\epsilon}_{i}\mid \mathbf{r}-\mathbf{x}\mid}dx+\mathcal{O}\left(\frac{{h}^{2}}{{a}^{2}}\right),$$

(3.14)

where

$${q}_{i1}=\gamma \frac{a}{{r}_{s}}q,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{q}_{i2}=\frac{q}{a}\frac{\gamma (1+\gamma )}{2}{\left(\frac{x}{{r}_{i}}\right)}^{-{\scriptstyle \frac{1-\gamma}{2}}},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\gamma =\frac{{\epsilon}_{i}-{\epsilon}_{o}}{{\epsilon}_{i}+{\epsilon}_{o}},$$

and

$$\begin{array}{l}{W}_{1}=1-\frac{1}{3\gamma}\left(1+\gamma -\sqrt{1-{\gamma}^{2}}\right)\phantom{\rule{0.16667em}{0ex}}\left(\frac{5}{2}-{\gamma}^{2}+\left(\gamma +\frac{1}{2}\right)\sqrt{1-{\gamma}^{2}}\right)\frac{h}{a},\\ {W}_{2}=\frac{1}{3\gamma}\left(1+\gamma -\sqrt{1-{\gamma}^{2}}\right)\phantom{\rule{0.16667em}{0ex}}\left(2-\gamma +\sqrt{1-{\gamma}^{2}}\right)\frac{h}{a},\\ {W}_{3}=1-\frac{1+\gamma -\sqrt{1-{\gamma}^{2}}}{3\gamma (1+\gamma )}\xb7\frac{1+4\gamma -3{\gamma}^{3}+(3{\gamma}^{2}-1)\sqrt{1-{\gamma}^{2}}}{2}\frac{h}{a},\\ {W}_{4}=\frac{1}{12}\left(1+\gamma -\sqrt{1-{\gamma}^{2}}\right)\phantom{\rule{0.16667em}{0ex}}\left(1-{\gamma}^{2}+(\gamma -1)\sqrt{1-{\gamma}^{2}}\right)\frac{h}{a}.\end{array}$$

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

$$I={\int}_{{r}_{i}}^{\infty}\frac{1}{\mid \mathbf{r}-\mathbf{x}\mid}{\left(\frac{x}{{r}_{i}}\right)}^{-{\scriptstyle \frac{1-\gamma}{2}}}dx.$$

(3.15)

By changing variable
${\scriptstyle \frac{{r}_{i}}{x}}={\left({\scriptstyle \frac{1-s}{2}}\right)}^{\tau}$ with *τ* >0, we have,

$$I=\tau {2}^{{\scriptstyle \frac{\gamma -1}{2}}\tau}{\int}_{-1}^{1}{(1-s)}^{\alpha}\xb7h(r,s;\tau )ds,$$

(3.16)

where

$$\alpha =\frac{(1-\gamma )\tau}{2}-1,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}h(r,s;\tau )=\frac{{2}^{\tau}{r}_{i}}{\mid {(1-s)}^{\tau}\mathbf{r}-{2}^{\tau}{\mathbf{r}}_{i}\mid}.$$

(3.17)

We use the Gauss-Radau quadrature [9,15] which is based on Jacobi polynomials ${P}_{n}^{(\alpha ,\beta )}(s)$ orthogonal on the interval (−1,1) with the weight function

$${w}^{(\alpha ,\beta )}(s)={(1-s)}^{\alpha}{(1+s)}^{\beta},$$

i.e.,

$${\int}_{-1}^{1}{(1-s)}^{\alpha}{(1+s)}^{\beta}{P}_{j}^{(\alpha ,\beta )}(s){P}_{k}^{(\alpha ,\beta )}(s)ds={\delta}_{jk},$$

(3.18)

where *α* >−1, *β* >−1. Take *s*_{0} = −1, *α* = (1 − *γ*)*τ*/2 − 1 and *β* = 0, and denote the Jacobi-Gauss-Radau points and weights by *s _{m}*,

$$I\approx \tau {2}^{{\scriptstyle \frac{\gamma -1}{2}}\tau}\sum _{m=0}^{M}{\omega}_{m}h(r,{s}_{m};\tau ).$$

(3.19)

Thus, for the line image charges, we have,

$${\int}_{{r}_{i}}^{\infty}\frac{{q}_{i2}(x)}{4\pi {\epsilon}_{i}\mid \mathbf{r}-\mathbf{x}\mid}dx\approx \sum _{m=0}^{M}\frac{{q}_{m}}{4\pi {\epsilon}_{i}\mid \mathbf{r}-{\mathbf{x}}_{m}\mid},$$

(3.20)

$${\int}_{{r}_{i}}^{\infty}\frac{{q}_{i2}(x)ln{\scriptstyle \frac{x}{{r}_{i}}}}{4\pi {\epsilon}_{i}\mid \mathbf{r}-\mathbf{x}\mid}dx\approx \sum _{m=0}^{M}\frac{{q}_{m}ln{\scriptstyle \frac{{x}_{m}}{{r}_{i}}}}{4\pi {\epsilon}_{i}\mid \mathbf{r}-{\mathbf{x}}_{m}\mid},$$

(3.21)

where

$${q}_{m}={2}^{{\scriptstyle \frac{\gamma -1}{2}}\tau -1}\gamma (1+\gamma )\tau {\omega}_{m}q\frac{{x}_{m}}{a},\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.16667em}{0ex}}{x}_{m}={r}_{i}{\left(\frac{2}{1-{s}_{m}}\right)}^{\tau}.$$

(3.22)

Therefore, multiple image charges approximation to the reaction field due to the point charge at **r*** _{s}* is obtained:

$${\mathrm{\Phi}}_{RF}(\mathbf{r})\approx {W}_{1}\frac{{q}_{i1}}{4\pi {\epsilon}_{i}\mid \mathbf{r}-{\mathbf{r}}_{i}\mid}-{W}_{2}\frac{{q}_{i1}}{4\pi {\epsilon}_{i}}r\frac{\partial}{\partial r}\left(\frac{1}{\mid \mathbf{r}-{\mathbf{r}}_{i}\mid}\right)+\sum _{m=0}^{M}\frac{({W}_{3}-{W}_{4}ln{\scriptstyle \frac{{x}_{m}}{{r}_{i}}}){q}_{m}}{4\pi {\epsilon}_{i}\mid \mathbf{r}-{\mathbf{x}}_{m}\mid},$$

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

$$r\frac{\partial}{\partial r}\left(\frac{1}{\mid \mathbf{r}-{\mathbf{r}}_{i}\mid}\right)=\frac{{r}_{i}^{2}-{r}^{2}}{2\mid \mathbf{r}-{\mathbf{r}}_{i}{\mid}^{3}}-\frac{1}{2\mid \mathbf{r}-{\mathbf{r}}_{i}\mid}.$$

(3.23)

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:

$$\mathrm{\Phi}(\mathbf{r})=\frac{q}{4\pi {\epsilon}_{i}\mid \mathbf{r}-{\mathbf{r}}_{s}\mid}+\sum _{m=0}^{M}\frac{{q}_{m}}{4\pi {\epsilon}_{i}\mid \mathbf{r}-{\mathbf{x}}_{m}\mid}.$$

(4.1)

The least squares method is used to find the charge strengths by minimizing the *L*_{2}-error of the reaction potential induced by the discrete images,

$$\text{Err}(\mathbf{q})=\sum _{n=1}^{N}{\left(\sum _{m=0}^{M}\frac{{q}_{m}}{4\pi {\epsilon}_{i}\mid {\mathbf{r}}_{n}-{\mathbf{x}}_{m}\mid}-{\mathrm{\Phi}}_{RF}({\mathbf{r}}_{n})\right)}^{2},$$

(4.2)

where **r*** _{n}* for

$$\text{Err}({\mathbf{q}}^{0})=min(\text{Err}(\mathbf{q})).$$

Taking the derivative with respective to *q _{l}*,

$$\frac{\partial \text{Err}}{\partial {q}_{l}}=0,\phantom{\rule{0.38889em}{0ex}}\text{for}\phantom{\rule{0.16667em}{0ex}}l=0,\cdots ,M,$$

we obtain a linear algebraic system with *M* equations:

$$\sum _{m=0}^{M}\left(\sum _{n=1}^{N}\frac{1}{\mid {\mathbf{r}}_{n}-{\mathbf{x}}_{m}\mid \phantom{\rule{0.16667em}{0ex}}\mid {r}_{n}-{x}_{l}\mid}\right)\frac{{q}_{m}}{4\pi {\epsilon}_{i}}-\sum _{n=1}^{N}\frac{{\mathrm{\Phi}}_{RF}({\mathbf{r}}_{n})}{\mid {\mathbf{r}}_{n}-{\mathbf{x}}_{l}\mid}=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}l=0,\cdots ,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 *r _{s}* can be interpolated from the data in the table.

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

$$\nabla \xb7({\epsilon}_{i}\nabla \mathrm{\Phi}(\mathbf{r}))=-{\rho}_{1}(\mathbf{r}),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}r\le a,$$

(5.1)

$$\nabla \xb7(\epsilon (r)\nabla \mathrm{\Phi}(\mathbf{r}))-{\lambda}^{2}{\epsilon}_{o}\mathrm{\Phi}(\mathbf{r})=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}a<r\le b,$$

(5.2)

$$\mathrm{\Delta}\mathrm{\Phi}(\mathbf{r})-{\lambda}^{2}\mathrm{\Phi}(\mathbf{r})=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}r\ge b,$$

(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 (*a* ≤ *r* ≤ *b*) is still used as before. Again, only one source charge at location **r*** _{s}* = (

$$\mathrm{\Phi}(r,\theta )=\{\begin{array}{ll}\sum _{n=0}^{\infty}\left[\frac{q}{4\pi {\epsilon}_{i}}\frac{1}{{r}_{s}}{\left(\frac{r}{{r}_{s}}\right)}^{n}+{B}_{n}{r}^{n}\right]{P}_{n}(cos\theta ),\hfill & 0\le r\le {r}_{s},\hfill \\ \sum _{n=0}^{\infty}\left[\frac{q}{4\pi {\epsilon}_{i}}\frac{1}{r}{\left(\frac{{r}_{s}}{r}\right)}^{n}+{B}_{n}{r}^{n}\right]{P}_{n}(cos\theta ),\hfill & {r}_{s}\le r\le a,\hfill \\ \frac{1}{\sqrt{\epsilon (r)}}\sum _{n=0}^{\infty}[{C}_{n}{k}_{n}(\lambda \sqrt{{\epsilon}_{o}}r)+{D}_{n}{i}_{n}(\lambda \sqrt{{\epsilon}_{o}}r)]{P}_{n}(cos\theta ),\hfill & a\le r\le b,\hfill \\ \sum _{n=0}^{\infty}{A}_{n}{k}_{n}(\lambda r){P}_{n}(cos\theta ),\hfill & r\ge b,\hfill \end{array}$$

(5.4)

where *i _{n}*(

$$\{\begin{array}{l}{i}_{0}(r)=\frac{sinh(r)}{r},\hfill \\ {i}_{1}(r)=-\frac{sinh(r)}{{r}^{2}}+\frac{cosh(r)}{r},\hfill \\ {i}_{n}(r)={i}_{n-2}(r)-\frac{2n-1}{r}{i}_{n-1}(r),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{for}\phantom{\rule{0.16667em}{0ex}}n>1,\hfill \end{array}$$

(5.5)

and

$$\{\begin{array}{l}{k}_{0}(r)=\frac{\pi}{2r}{e}^{-r},\hfill \\ {k}_{1}(r)=\frac{\pi}{2r}{e}^{-r}\left(1+\frac{1}{r}\right),\hfill \\ {k}_{n}(r)={k}_{n-2}(r)+\frac{2n-1}{r}{k}_{n-1}(r),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{for}\phantom{\rule{0.16667em}{0ex}}n>1.\hfill \end{array}$$

(5.6)

The expansion coefficients *B _{n}* are determined by the continuity of the potentials and the fluxes normal to the boundaries (2.17a)–(2.17d), which are given by

$${B}_{n}=\frac{{qa}^{-1}{\left({\scriptstyle \frac{rs}{{a}^{2}}}\right)}^{n}}{4\pi {\epsilon}_{i}}\frac{{S}_{1}}{{S}_{2}},$$

(5.7)

where

$${S}_{1}=-{i}_{n}({t}_{2})\phantom{\rule{0.16667em}{0ex}}\left(\alpha {k}_{n}(u)-b(\alpha +b\beta )\lambda {k}_{n}^{\prime}(u)\right){s}_{11}+\left((2+n)\alpha +a(1+n)\beta \right){i}_{n}({t}_{1}){s}_{12}-{t}_{2}(\alpha +b\beta ){i}_{n}^{\prime}({t}_{2}){k}_{n}(u){s}_{11}+{t}_{1}(\alpha +a\beta ){i}_{n}^{\prime}({t}_{1}){s}_{13},$$

(5.8)

$${S}_{2}={i}_{n}({t}_{2})\phantom{\rule{0.16667em}{0ex}}\left(\alpha {k}_{n}(u)-b(\alpha +b\beta )\lambda {k}_{n}^{\prime}(u)\right){s}_{21}+(-a+n\alpha +an\beta ){i}_{n}({t}_{1}){s}_{12}+{t}_{2}(\alpha +b\beta ){i}_{n}^{\prime}({t}_{2}){k}_{n}(u){s}_{21}-{t}_{1}(\alpha +a\beta ){i}_{n}^{\prime}({t}_{1}){s}_{13}.$$

(5.9)

Here,

$${s}_{11}=\left((2+n)\alpha +a(1+n)\beta \right){k}_{n}({t}_{1})+{t}_{1}(\alpha +a\beta ){k}_{n}^{\prime}({t}_{1}),$$

(5.10a)

$${s}_{12}=-b(\alpha +b\beta )\lambda {k}_{n}({t}_{2}){k}_{n}^{\prime}(u)+{k}_{n}(u)\phantom{\rule{0.16667em}{0ex}}\left(\alpha {k}_{n}({t}_{2})+{t}_{2}(\alpha +b\beta ){k}_{n}^{\prime}({t}_{2})\right),$$

(5.10b)

$${s}_{13}={k}_{n}({t}_{2})\phantom{\rule{0.16667em}{0ex}}\left(\alpha {k}_{n}(u)-b(\alpha +b\beta )\lambda {k}_{n}^{\prime}(u)\right)+{t}_{2}(\alpha +b\beta ){k}_{n}(u){k}_{n}^{\prime}({t}_{2}),$$

(5.10c)

$${s}_{21}=(\alpha -n\alpha -an\beta ){k}_{n}({t}_{1})+{t}_{1}(\alpha +a\beta ){k}_{n}^{\prime}({t}_{1}),$$

(5.10d)

and

$$u=\lambda b=\lambda (a+h),\phantom{\rule{0.38889em}{0ex}}{t}_{1}=a\sqrt{{\epsilon}_{o}}\lambda ,\phantom{\rule{0.38889em}{0ex}}{t}_{2}=b\sqrt{{\epsilon}_{o}}\lambda .$$

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].

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

$${\mathrm{\Phi}}_{RF}^{\mathit{exact}}\approx \sum _{n=0}^{N}{B}_{n}{r}^{n}{P}_{n}(cos\theta ),$$

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

First, we compute the *L*^{2}-norm relative error of both approaches with the thickness of the buffer zone set as *h*= *b*−*a*=0.1 and 0.01, where the error is defined by

$${\left|\right|E\left|\right|}_{2}=\frac{{\left({\sum}_{i,j}\mid {\mathrm{\Phi}}_{RF}({r}_{i},{\theta}_{j})-{\widehat{\mathrm{\Phi}}}_{RF}({r}_{i},{\theta}_{j}){\mid}^{2}\right)}^{{\scriptstyle \frac{1}{2}}}}{{\left({\sum}_{i,j}\mid {\mathrm{\Phi}}_{RF}({r}_{i},{\theta}_{j}){\mid}^{2}\right)}^{{\scriptstyle \frac{1}{2}}}}$$

(6.1)

for Φ* _{RF}*(

Semi-analytical vs. least squares image charge methods. Accuracy of the reaction field for the source charge located **r**_{s} = (*r*_{s}, 0, 0) with discrete image charges number *M*+1 = 2 and 10 including the Kelvin image. Left: *h*=0.1; right: *h*=0.01.

Semi-analytical vs. least squares image charge methods. Left: accuracy with the thickness of the buffer layer *h* = *b* − *a* with the source location *r*_{s} = 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.

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.

Locations and strengths of image charges by the least squares approach with *a*=1, *h*=0.1, *ε*_{i} = 2. The first image charge *q*_{1} is located at the Kelvin point *x*_{1}, and the second image charge *q*_{2} with location *x*_{2}.

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 *L*^{2} 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.

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
$\theta =0,{\scriptstyle \frac{\pi}{10}},\cdots ,{\scriptstyle \frac{19\pi}{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. 8–10 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 *L*^{2} 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.

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.

Relative error of the least squares approximation compared with the exact solution for the three-layer model. Left: accuracy vs. source charge location *r*_{s} for different dielectric constants *ε*_{o} of bulk solvents; right: accuracy vs. the number of **...**

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.

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).

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]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |