|Home | About | Journals | Submit | Contact Us | Français|
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  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  extended this study to a dipole. Both Born and Onsager models are special cases of the results of Kirkwood  and Tan-ford and Kirkwood , 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  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  and its improved version  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  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 . The simplest of these models assigns an effective dielectric constant between all pairs of charged particles as a linear function  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  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 (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.
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 (r≥b), 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 = b−a, 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:
where ρ1(r)= Σj qjδ(|r−rj|) is the charge distribution inside the interior cavity. Here rj is the location of charge qj, and Δ are gradient and Laplace operators, respectively.
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)=qδ(|r−rs|), where rs =(rs,0,0) under a spherical coordinate system r= (r,θ,). Because of the azimuthal symmetry, the potential Φ(r,θ) in the exterior layer can be expressed as
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≤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 (|r|>a):
where R is the distance between the field and source points,
By using the expansion of the reciprocal distance
the total potential inside the interior sphere can be rewritten as
In the intermediate layer with continuous dielectric constants ε(r), the potential at r is represented in the following form:
where coefficients Cn(r) and Dn(r) are continuous functions of r, which can be determined by two differential equations . 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,
By using the continuities of the potential and the flux normal to the interface
and the orthogonality of the Legendre polynomials, we have
For an infinitely thin shell, can be replaced by Cn(r+dr) and by Cn(r). The same steps can be taken for and εl, εl+1. As the function f (r+dr) can be expressed as
two coupled differential equations can then be obtained:
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 in Eq. (2.2), then the potential holds
Then, by a simple transform, the potential in the layer can be expressed by
To match this result, the interpolation with two harmonic functions 1 and are used to obtain
where α,β can be determined by the continuity of the dielectric permittivity at the interface (Fig. 3), say,
Using the orthogonality of Legendre polynomials, we can obtain An,Bn,Cn and Dn. We are only interested in the interior layer, where
and, T1 and T2 are given by
In particular, when h=0, it reproduces the Kirkwood’s result :
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 (h2) in terms of the thickness h of the buffer layer as
and Ck for k = 1,···,4 are constants depending on a, h and γ:
For the second term in (3.1), we have
To obtain an image formula of the third term in (3.1), we denote function
Multiplying both sides by and differentiating with respect to x, we have
Taking the integration from ri to ∞, due to f (r,x) → 0 as x → ∞, we transform the series into a line image:
Similarly, the fourth term in (3.1) can also be represented by a line image as follows. Let
Calculating the differentiation in g(r,x), we have
Multiplying g by x and taking the derivative with respect to x yields:
Again, taking the integration from x to ∞, because of
with vector t = (t,0,0). With a simple calculation, we also have
We thus obtain the line image from the Kelvin point ri to infinity:
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,
In order to discretize the line image charges into multiple point charges, we use a numerical quadrature rule to approximate the following integral
By changing variable with τ >0, we have,
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 . The integral I is then discretized as
Thus, for the line image charges, we have,
Therefore, multiple image charges approximation to the reaction field due to the point charge at rs is obtained:
where the derivative of the Kelvin image with respective to r can be expressed as:
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:
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,
where rn for n = 1,···, N are N field points inside the interior sphere. Our least squares procedure amounts to finding , such that
Taking the derivative with respective to ql,
we obtain a linear algebraic system with M equations:
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.
Extension of the three-layer model to solvents with salt effects is based on the following model,
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 rs = (rs,0,0) for rs < a is considered, i.e., ρ1(r) = qδ(r − rs). Due to the azimuthal symmetry, the potential at an observation point r= (r,θ,) can be represented by
where in(r) and kn(r) are the first and the second kind of modified spherical Bessel functions  defined recursively by
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 ε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 , and the cutoff error for calculating the exact reaction field is chosen to be 10−8, i.e.,
where the (N+1)-term and (N+2)-term are both less than this number.
First, we compute the L2-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
for ΦRF(ri, θj) and 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 . 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 (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.
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.
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.
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 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 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.
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).