J . Phys. Chem. 1993,97,4855-4864
Interaction of a Biomolecule with Mobile Ions in Aqueous Solution. Comparison of Three Fast Approximate Methods with the Direct Solution of the Nonlinear Poisson-Boltzmann Equation Yury N. Vorobjevt and Harold A. Scheraga' Baker Laboratory of Chemistry, Cornel1 University, Ithaca, New York 14853-1301 Received: December 11, 1992;In Final Form: February 16, I993
To treat the electrostatic effect of the solvent and mobile ions on a biomolecule solute, it is necessary to have a collection of approximate methods and to choose the particular method according to the available computer resources and the complexity of the system being modeled. In the present paper, calculation of the solvation energy, the solvation forces, and the interactions with mobile ions are considered, and three fast approximate methods are introduced. The first one is the thin ion layer (TIL) approximation, which approximates the real profile of the mobile ion charge density by a hypothetical set of surface-charge ions condensed on a particular surface around the biomolecule solute. The TIL approximation has one empirical parameter, the distance ATIL between the molecular surface and the surface carrying the ionic charge. The second method involves an approximation of the Debye diffuse ion layer (DDIL), which does not involve any empirical parameters. The DDIL method treats the mobile ion density profile as a diffuse ion shell of the mobile ions as if it were condensed on the molecular surface. Both the TIL and DDIL methods are based on the calculation of a surface charge density of the condensed ions by the fast iterative boundary element method (IBE). The third method (the N S A or normalized superposition approximation) is a modification of the procedure for calculating the mobile ion density, based on a simplification of a formal statistical mechanical expression for the conditional probability of the mobile ion density by applying the Kirkwood superposition approximation (SA), and an analytical approximation of the pair correlation functions between the ions and the atoms of a biomolecule solute. The approximate methods are tested by comparing the electrostatic interaction energy between the solute molecule and the mobile ions, and the self-electrostatic energy of the ion-charge-density profile, with the results of the solution of the nonlinear Poisson-Boltzmann (NLPB) equation for a spherical Tanford-Kirkwood model and for oligopeptides of lysine with ionized side chains. It is found that the TIL method gives good agreement with the results of the NLPB equation, and the dependence of the empirical parameter ATILcan be approximated by a simple expression, depending on the Debye length and the partial charges of the mobile ions and the biomolecule solute. The fast approximate TIL and DDIL methods are able to estimate the electrostatic energy difference between conformations of peptides with errors within 25%, and the accuracy of the N S A method is within 15%, relative to the solution of the NLPB equation. We demonstrate that the solvent contribution from a polar ionic solution stabilizes the a-helical conformation over the extended one, for decalysine with ionized side chains, and that a less polar solvent with mobile ions stabilizes the a-helix more effectively than a more polar solvent. The quantitative aspect of this behavior is related to the change of charge distribution and accounts for the observed stabilization of a charged poly@-lysine) a-helix in alcohol-water mixtures, whereas the helix is disrupted in water. The CPU times required by the TIL, DDIL, and NSA methods and the NLPB equation method to calculate the energy of interaction of a charged decalysine with the mobile ions are, on the average, in the ratio 1:1.5:2:15, respectively.
Introduction
Computation of ionic effects on charged biomolecules has been a subject of continuous interest because of the importance of the solvent and ionic environment for structural stability and conformational transitions as well as for intermolecular interactions and substrate-binding equilibria. For example, at high pH in water, a poly(L-lysine) chain appears to be in an a-helical conformation, and at low pH adopts a statistical-coil conformati0n.l But, in methanol and other alcohols, as shown by Epand and Scheraga,* poly( L-lysine) seems to be completely helical, even at neutral pH where this polymer is highly charged, as shown by Joubert et al.3a and by Liem et al.3b Stabilization of the a-helical conformation for charged poly(L-lysine) in alcohol solution is probably due in the main to interactionswith the solvent and counterions. Ionic effects also play a central role in the E Z and E A conformational transitions of DNA.4 The approaches proposed for computing ionic effects on biomolecules vary widely according to the level of structural Author to whom correspondence and reprint requests should beaddressed.
+
On leave from the Novosibirsk Institute of Bioorganic Chemistry,
Novosibirsk 630090,Russia, 1992-1993.
modeling of the biomolecule-solvent-ion system. In early studies, the water solvent was treated as a dielectric continuum with a uniform dielectric constant, and the ions were considered either as simple point particles or hard spheres. The solvent- and ionaccessible surface of globular proteins has been modeled as the surface of a sphere with a low value of the dielectric constant in~ide53~ or as an infinitely-long cylinder in the case of a-helical proteins7 or DNA.8 The ion-density distribution around the biomolecule solute has been calculated by analytical5,6,9 or numerica17J0J1 solutionsof the Poisson-Boltzmann (PB) equation, as well as by Monte Carlo (MC) techniques.I1-I5 Monte Carlo calculations of the ion density profile and the potential of mean force between an ion and a protein provide evidence that, for low ionic strength (51 M), the nonlinear PB equation is an excellent approximation and gives results that are in good agreement with direct computer simulations for 1:l and 2: 1 electrolytes.llJ5 An explicit microscopic treatment of both the biomolecule and the solvent molecules and ions is possible in large-scale MC or molecular dynamics (MD) simulations, as reported by Clementi and CorongiuI6and by Berendsen et alei7 But such calculations are extremely costly, and can be carried
0022-3654/93/2097-4855%04.00/0 0 1993 American Chemical Society
4856 The Journal of Physical Chemistry, Vol. 97, No. 18, 1993
Vorobjev and Scheraga
only the solvent reaction potential and the solvation energy.20q22,24 out practically only for simplified systems (biomolecule-waterTherefore, it is reasonable to develop a collection of methods to counterions, i.e. without added salts), because, to simulate a 0.1 provide an approximate treatment of the interaction of a M NaCl solution around a protein solute, one has to consider biomolecule with mobile ions in aqueous solution; thus, one may several thousand water molecules and ions. In principle, a largechoose a particular method depending on the available computer scale computer simulation can provide such an exact treatment resources and the complexity of the system being modeled. of electrostatic interactions in solution but, at the present time, The goal of the present work is to test several computationallythey are prohibitively computationally expensive. fast approximate methods for calculating the energy of interaction In order to surmount the high cost of explicit models of a between a biomolecule solute and the mobile ions in solutions of solvent and mobile ions, much renewed interest has developed in polar solvents by comparing the results with a solution of the the last five years in the use of simple continuum models18-22 NLPB equation. We introduce three new approaches which are based on the PB mean-field theory for the ion distribution, with computationally considerably faster than solving the NLPB a continuum homogeneous solvent surrounding the solute molecule equation. The first one is an approximation of the real profile that is characterized by an explicitly-treated atomic structure of the mobile ion charge density by a hypothetical set of surfaceand a solute-solvent interface. The inside molecular volume is charge ions condensed on a particular surface (not the molecular considered to have a homogeneous low dielectric constant, while surface) around the biomolecule solute [a thin ion layer (TIL) the outside volume, filled by the solvent, usually has a high approximation]. The second one is an approximation of the Debye dielectric constant. The border between the two dielectric media is treated as a solvent-accessible molecular s ~ r f a c e , whose ~ ~ - ~ ~ diffuse ion layer (DDIL); it treats the mobile ion density profile as a diffuse ion shell of the mobile ions as if it were condensed area is calculated by the procedure of ref 25a. The interaction on the molecular surface. The third method to calculate the ion of the surface solvent-induced charges and the solvent ion charge density is the SA procedure of Klement et a1.37-38We show here density with the biomolecule solute charges determines the energy that the original versi0n~~J8 of the SA method must be corrected of interaction of the biomolecule solute with its environment and because it does not satisfy the electroneutrality condition for the the forces acting on the atoms of the biomolecule solute. To solve sum of the total ionic charge and the total charge of the the linear or nonlinear PB equation for a molecule embedded in biomolecule solute; our corrected version will be referred to as an electrolyte solution, two different sets of numerical methods a normalized superposition approximation (NSA). We found have been developed: the finite difference (FD) method of Honig that the TIL and DDIL approximations and the NSA approach et al.19 and the approaches based on the boundary element (BE) give results (for the total energies and forces acting on the atoms method elaborated by Rashin20a.bfor the nonlinear PB equation of a biomolecule) that agree within 20% with the results of a and by Zauhar and Morgan (for the Poisson equation, Le., only solution of the NLPB equation for a model spherical system and for the case of zero ion concentration),21 and by Yoon and for real biomolecules (charged oligolysines in both the a-helical Lenhoff,22aand by Juffer et a1.22bfor the linear PB equation. An and extended conformations). In the next sections, we describe effectiveversion of the iterative boundary element (IBE) method the theory for the different approaches used in this work and the for the solution of the nonlinear PB equation has been reported results of the calculations. by Vorobjev et al.24 Theory A set of very elegant calculations has recently shown that a continuum model reproduces the solutesolvent free energies Continuum Model. This is a model in which a biomolecule obtained by using a microscopic treatment of the ~ o l v e n t . *It~ . ~ ~ solute is represented by a cavity of low dielectric constant Di is likely that the success of such continuum models is in part due embedded in a continuum solvent medium of high dielectric to the unique nature of the behavior, at the molecular level, of constant Do. The biomolecule soluteis separated from the solvent the water in the region of the dielectric boundary.26 At the by a boundary which is defined to be the smooth surface S traced boundary, the dielectric constant varies very rapidly. According by the inward-facing part of a probe (with the radius of a water to direct calculations with the aid of liquid theory, the effective molecule) as it rolls over the atoms of the biomolecule solute.23-25 thickness of a shell around a solute molecule, containing the This boundary confines the biomolecule solute charge distribution, solvent-induced charges, is about the width of the first peak of which is typically represented as a set of fixed charges (qk]located the pair correlation function between the water molecules and at atomic centers (rk]. The polar solvent contains n ionic species, the atoms of the solute molecule.28 Because the first peak of this i = 1, ..., n, with bulk number densities p?. Now, we consider pair correlation function is quite sharp (the average width is about the different approaches to calculate the energy of interaction of 0.5 A), the continuum model with an infinitely-thin dielectric the mobile ions with the biomolecule solute. We begin by boundary has some statistical mechanical justification. Calcusummarizing our earlier treatment24 of the nonlinear Poissonlations of the potential of mean force between ions in aqueous Boltzmann equation. ~ o l u t i ohave n ~ ~also ~ shown ~ ~ that water completely screensvacuum Nonlinear Poissoehltzmann Equation. This equation treats Coulombic interactions within one hydration shell. Such behavior the mobile ions as a charge density distributed, within the is well represented by a simple continuum model with an infinitelyframework of the Boltzmann law, according to the mean-field thin dielectric boundary. These considerations were the basis for potential of the ions and the electrostatic potential of the using a shell model to treat hydration in conformational energy biomolecule solute. This total electrostatic potential can be calculations on polypeptides and p r ~ t e i n s . ~ I - ~ ~ writter12~as Finally, the ion density can be calculated by an approach based on thereferenceinteractingsitemodel (RISM)15134.350f thetheory of liquids. A numerical solution of the RISM integral equation where @,,I is the potential due to the solute charge distribution, is still quite expensive; that is why approximate methods to model 9,is the reaction potential, arising from the linear response of ion density on the basis of the Kirkwood superposition approxthe solvent dielectric medium to the biomolecule solute charge imation (SA)36and the mean spherical model (MSM) for the distribution, and @ion is the mean-field potential of the ion pair correlation functions between ions and the atoms of the distribution. The sum of the potentials, I ,@, + ar,satisfy the biomolecule solute have been proposed by Klement et a1.37J8 Poisson equation, namely For a biomolecule of arbitrary shape, the calculation of the ion density by a solution of the general three-dimensional NLPB equation is still a difficult task and requires more computer resources than the solution of the Poisson equation to calculate where rk is the position of a set of fixed charges qk within the
The Journal of Physical Chemistry, Vol. 97, No. 18, 1993 4857
Interaction of a Biomolecule with Mobile Ions solute cavity, and D(r) is the dielectric constant at point r. We solve eq 2 by the effective iterative boundary element method24 (IBE), which is a computationally more efficient version of the boundary element method of Zauhar and Morgan.z',39-40The is given straightforwardly by molecular potential (3)
is obtained from the surface
and the reaction potential integral
ssm a(s) ds
=
(4)
where the surface charge density a(s) at s is determined by the integral equation24
=fJs
(5)
It - 4 3
where the constant f is f=-- 1 Di--DO
277 Di
+ Do
and n(t) is the vector nornlal to the surface S at the point t. The last term in eq 5 is the normal component of the electrostatic field due to the solute charge distribution. In the IBE method, the integral appearing in the first term of eq 5 is replaced by a discrete sum over a finite number of boundary elements which tesselate the molecular surface. The mean-field potential due to the ion distribution is given by the NLPB equation; thus
V.D(r)V@ion(r)+ 477piOn(r) = 0 where pion(r) is the total ionic charge density given by
(7)
The index 1 in eq 8 runs over the different ion species in solution, and Zl and pp are the ion charge and bulk density, respectively, of the Ith ionic species. Unbi(r)is the hard-sphere potential energy of the nonbonded interaction of ion type 1 with the biomolecule solute. This nonbonded term serves in part to account for the finite size of the ions and is usually taken to be independent of ion type (see, for example, refs 19 and 20). The potential due to the ion distribution @io&) appearing in eq 8 can be decomposed into a sum of a direct term, plus a reaction term arising from contributions due to solvent response, so that
where the direct term is given by
equation similar to eq 5, viz.
where s and p are integration variables. The last term in eq 12 is the normal component of the field due to the mobile ions at a point t on the molecular surface. In using the IBE method, expressions for the surface charge densities in both eqs 5 and 12 can be written as a set of linear equations in matrix form;41thus
Ka = b (13) The matrix K is a function only of the geometric properties of the boundary elements that cover the molecular surface S and needs to be computed only once. A complete definition of the matrix elements of K,as it is used in this work, is given in ref 41. Theright-handsideofeq 13 isa functionofthenormalcomponent of the field, due toeither the biomoleculesolutecharge distribution when solving eq 5 or the mobile ions when solving eq 12. Since @ion in eq 7 is coupled to pion(r)through eq 8, eqs 8-12 must be solved iteratively, for example by setting aion(r)initially to zero. A similar approach has been used by Rashin and Malinsky20bfor the calculation of the mobile ion distribution around a rodlike polyelectrolyte. Our implementation of the iteration method for solving the NLPB equation is general for an arbitrarily-shaped biomolecule solute; it is effective because we use a particularlyconstructed boundary-element-based grid (BEBG) to represent the volume outside the biomolecule solute as a collection of nonequivalent prismatic and pyramidal cells filling the shells between successively expanded surfaces around the biomolecule solute. Use of such a BEBG minimizes numerical errors in describing the mobile ions density and the total number of cells. Further details of the numerical implementation of the IBE method are given in ref 24. Energy and Forces The electrostatic free energy of solvation or the interaction of the solute molecule with a solvent (e.g. pure water), Wsolvrwhich is equal to the total work associated with moving the molecular charge distribution from vacuum to a solvent, is given (with the aid of eq 4) by
where the integration is carried out over the molecular surface, and the index k runs over all the point charges of the biomolecule. The forces acting on the atoms of a biomolecule solute due to the polarization of the solvent can be considered as solvation forces. These electrostatic solvation forces reflect the change of the electrostatic free energy of solvation in response to infinitesimaly small displacementsof the atoms of a biomoleculesolute. Thus, the electrostatic solvation forces acting on an Ith atom can be calculated as the derivative of the solvation free energy Ws0l, over the atomic coordinates (rl) of the biomolecule.
FF = -6 Wsolv/6r, where the limit of integration V,,, is used to indicate that numerically we carry out the integration over a finite volume. The reaction potential @ir(r)can be expressed in terms of a surface charge density ai,(s) induced on the molecular surface, such that
The surface charge density ai,(s) can be obtained from an integral
(15)
According to classical electrostatic theory,4*q4)the electrostatic forces acting on a body embedded in a dielectric medium may be considered by two equivalent and independent approaches, either as a derivative of the electrostatic energy, eq 15, or as the integral of the stress tensor over the surface around the body. These two approaches should not be mixed as is done in ref 40. If an Ith atom of the biomolecule moves over a distance brl, then the coordinates (sj of the surface elements of the surface S of the biomolecule, and the surface charge density u(s), will also
4858
The Journal of Physical Chemistry, Vol. 97, No. 18, 1993
be changed; hence, the total force Ftotcan be represented as
F:,, = F ; ~+~ F;
+ F,"
(16) where Ftiris the direct electrostatic interaction of the Ith atom with the induced polarization charges
where q/ is the charge on the Ith atom of the biomolecule solute. There are two additional indirect solvation-force components. The first is a force F; due to a change in the surface of the biomolecule and, hence, in the positions s of the induced charges
4s)
where the index k runs over all atoms of the biomolecule solute and q k is the charge on the kth atom. The second component F," is due to a change of the polarization charges a(s) on the surface of the biomolecule
1% - 91
(19)
Becausethe induced charge density a(s) satisfiesthe conservation law, application of Gauss' la^^^.^^,^^ results in
where Qmolis the total charge of the biomolecule. The derivatives of the charge density a(s) satisfy the condition
Js-
which itself changes position. It is easy to show that the derivatives (6s/6rl)ak do not exceed I/z for saddle elements of atoms (expressed in terms of their van der Waals radii) and a probe radius of the water molecule of 1.4A. Hence, eq 22 overestimates the derivatives for the saddle and the concave surface elements adjacent to the kth atom. We determine the partial surface A s k as a collection of points s on the surface S of a biomolecule, for which a distance 1s - rkl scaled by the atomic radius rk" is minimal, when k runs over all the atoms of the biomolecule. Such a determination of d s k obviously underestimates the area of the surface of the biomoleculewhich moves simultaneouslywith atom k and, hence, leads to partial compensation of errors, which are produced by the overestimation of the derivatives for saddle and concave surface elements, after integration over the surface S in eq 18. Introducing the approximation of eq 22 into eq 18, and combining the Ffr and F,S components of the total force Ftot, we obtain the final expression, which is suitable for a calculation within the framework of the IBE method.
a(s)(r/- s) ds
F:,, =
1 -4/xJ.., 2
k#/
-
ir, - si3
1
@J(s)/6rJ ds
1
F," = - ; F q k s S
Vorobjev and Scheraga
ds = 0
Thus, because of eq 2 1, it is reasonable to assume that the indirect solvation force acting on any atom, arising from the induced polarization charges of different parts of the surface of the biomolecule, will be partially canceled in eq 19; hence, it follows that the component F," of the indirect polarization force will be of second order compared to the direct one, at least for the most important case of charged surface atoms of a charged biomolecule, and we neglect it. The derivatives of the coordinates of the molecular surface elements (6s/6rJ in eq 18 can be approximated by the simple expression
where I is a diagonal 3 X 3 unit matrix, &/is the Kronecker delta, and A s k is that part of the surface contributed by the kth atom of the biomolecule, and it is assumed that the portion A S k moves over the same distance brk as the kth atom, while the other parts of the surface S are fixed. As follows from the anatomy of a molecular surface,23eq 22 is exact for convex elements of the surface but it is not valid for saddle and concave elements of a surface, whose coordinates depend on the movement of the kth atom in a complex manner. From a general point of view, the most flexible elements of a molecular surface are the convex elements, because their coordinates depend on the coordinates of one atom. The positions of saddle and concave elements of the surface are determined by the inward-facing part of the surface of a probe sphere as it rolls over the biomolecule solute tangent to two or three atoms of the biomolecule solute;23consequently, these positions are more immobile, because they depend on the coordinates of two and three atoms, respectively, only one of
-xqkJsI
2kfl
a(s)(rk - S) ds
Irk
-4 3
(23)
In the first term, the direct interaction of the charge qr on the Ith atom with the surface charge u(s) induced on the surface AS/of the Ith atom is canceled exactly after combining the FPrand F; components of the total force Ftot. In the last term, all =0 for k # I, and hence the integration is carried out only over S = As/.The IBE method has some advantage over the treatment within the framework of the FD in the main, because the IBE methodconsidersthedielectricboundary and the positions of the induced polarization charges more precisely. It can be seen that, even if the fixed charge q/ on the Ith atom of the biomolecule solute is equal to zero, the second term of eq 23 remains nonzero if q k on some of the atoms of the biomolecule solute is nonzero; hence eq 23 describes an indirect solvation force acting on an uncharged surface atom 1of a biomolecule due to the interaction between the polarization charges induced on the surface of the atom 1 and the fixed charges of other atoms of the biomolecule. Now consider the presence of mobile ions. An expression for the total Helmholtz electrostatic free energy of the mobile ions in the electrostatic field of the biomolecule is given in refs 45 and 46. This expression consists of three terms
where is the energy of electrostatic interaction of the ions with the biomolecule solute (including the effect of the reaction field from the solvent), ITislris a self-electrostatic energy of the ion cloud, andAidalis the free energy of noninteractingions treated as an ideal gas distributed with the same density profile pi&). The direct electrostatic interaction E,,,, of the ion cloud with the biomolecule is
Thedirect electrostaticforceof the ion cloud actingon thecharges q k of the biomolecule is
The expression for the self-electrostatic energy E,,,, of the ion
Interaction of a Biomolecule with Mobile Ions
The Journal of Physical Chemistry, Vol. 97, No. 18, 1993 4859
cloud is
The ideal part, Aiidell,of the total free energy is the free energy of the noninteracting mobile ions in the ion
where XI is the thermal wavelength and p/(r) is the density of ions of type I. If the bulk density of the mobile ions is constant and independent of. the conformation of the biomolecule solute,20b then the expression for the excess free energy Mildeal over that for a uniform bulk ion distribution, taken as a reference state, is
Obviously, this value of Miidea, is due to a change of the configurational entropy of the distribution of the mobile ions. If the ion densities pion(r)are the solution of the NLPB equation, then, taking eqs 8,25, and 21 into account, and recognizing that p/(r)U,,b/(r) = 0 for a hard-sphere potential, eq 29 becomes
Finally, the total free energy of interaction of the mobile ions with the biomolecule solute, in the NLPB approximation, is
Normalized Superposition Approximation Klement et a1.37,38suggested an approximate method for calculating the total ion density around a biomolecule solute by using the formal expression of the conditional probability density of the mobile ions in the field of a set of charged particles fixed in a particular configuration (atoms of a biomolecule) and simplifying it by applyingKirkwood'ssuperpositionapproximation to express the many-particle ion-atom correlation function for the biomolecule solute as a product of ion-atom pair correlation functions for the biomolecule solute. This led to a simple expression for the ion charge density pionSA(r)
where the index I runs over all ion species in the solvent,the index k runs over the atoms of the biomolecule solute, ZIand pP are the ion charge and bulk density, respectively, and g/k is the spherically symmetrical correlation function between mobile ion 1 and atom k. The approximate equation (32) becomes exact only far from the biomolecule solute, where all distances (r - rk) become large. In principle, the pair correlation functions, g/k, may be calculated by the RISM integral equation method,28s34,35 which, however, is quite expensive for a large biomolecule solute. Another possibility is to adopt the approximate analytical treatment of the ion-atom pair correlation functions, g/kr calculated by a mean spherical approximation (MSA) for a restricted primitive model of electrolyte^.^^ According to the MSA, the pair correlation function, g/k, between two charged particles has the form gIk(r)= g,,h"(r) e x p [ - z T
h,;'(r)]
for r Id
(33)
g,k(r) = 0 for r > d where d is the sum of the hard sphere radii of the two particles, g/khs(r)is the pair correlation function of a system of hard spheres, and the function h/kel(r)describes the electrostatic correlations.
The correlation function g/khs(r)can be calculated on the basis of an analytical solution of the Percus-Ye~ick~~ equation with a correction suggested by Verlet and Weis.49 An analytical expression for the function hjkel(r) has been evaluated within the framework of the MSA theory of Waisman and Lebowitz by Blum et al.50251
where
x = [ ( l + k,d)"2-
11
(35) and kD-l is the Debye length. Because the ion charge density, pionSA(r) calculated according to eqs 32-34, is only an approximation, it does not satisfy the general electroneutrality condition. The ion charge density pio&) must satisfy two basic physical conditions. The first is the condition that the ion charge density pion(r)must be equal to zero at large distances from the solute molecule. The superposition approximation, eq 32, satisfies this condition if the pair correlation function, glk(r), has the form of eq 33. The second is the general electroneutrality condition
where Qmolis the total charge of the biomolecule solute. The superposition approximation, eq 32, with thecorrelation function, g/k(r),does not satisfy the electroneutrality condition; hence, to achieve physical consistency, it is very important to normalize the ion density, pionSA(r),to satisfy the electroneutrality condition, eq 36, if one is to obtain a reliable calculation of the energy of interaction between the mobile ions and the biomolecule solute. Klement et a1.37,38did not mention normalization of their SA approximation. We refer to such normalized ion charge density as a normalized superposition approximation (NSA).
(37) The normalization factor knrNSAcan be found by directly integrating the mobile ion charge density, pionSA(r)over the outside volume of the biomolecule solute using our fine BEBGZ4 to minimize the numerical errors of the volume integration and dividing this result into -Qm0l. Another important modification of the method of Klement et a1.37,38 for calculating the ion density, tested in this work, consists of replacing the real distance f l k , between an ion 1 at the point ri and an atom k at the point rk, by the distance rls, where rls is the minimal distance between the ion I and the molecular surface S. Hence, instead of eq 34 for the function h/kel(r),we use
Such a modification of the pair correlation function takes into account the effect of exclusion of the mobile ions from the inside volume of the biomolecule solute.
Model of a Thin Ion Layer The effective thicknessof theion cloud in theoutward direction from the surface of the biomolecule solute is about kD-' in the linear PB approximation and less for both the NLPB equation and a MC simulation of the ion density profile.ll Multivalent ions and low temperature increase the tendency of the ion profile to be sharp. For the purpose. of calculating the energy of interaction of the biomolecule solute with the mobile ions, we can approximate the exact ion charge profile pion@)by an infinitely thinion layer (TIL) condensed on a particular surfaceSTll, which is outside the molecular surface. Thus, consider the ion charge
Vorobjev and Scheraga
4860 The Journal of Physical Chemistry, Vol. 97, No. 18, 1993
density model
where s is a point on the surface ST~L and is the surface charge density of the condensed ions. Formally, an ionic solution can be considered as a conducting medium or, from the point of view of electrostatics, as a medium with an infinite static dielectric constant DO.Hence, the surface charge density uic(s) of the condensed ions can be found as a solution of the Poisson equation for the continuum model in which thesurfaceSTlLserves as a border between the conducting medium beyond the surface S T ~and L the polar solvent between STIL and the molecular surface. Thus, there is a shell of thickness ATIL between the ion-condensed surface ST~L and the surface S of the biomolecule solute. Hence, the TIL surface charge density at a point t of the condensed ions uic(t) is the solution of the integral equation
where E, is the electrostatic field from the atomic charges of the biomolecule and the solvent charge density u(s) induced on the molecular surface S , and a(s) is a solution of eq 5 . The constant Jon is Jon
= -1/27r (41) satisfies the electroneutrality
The TIL ion charge density u&) condition
model can be considered as a diffuse ion density representation of the TIL model. Such a distributed ion density profile is an extension of the TIL model toward the more realistic treatment of the ion charge density profile
where r is an arbitrary point in the volume outside of the solute molecule, sr is the nearest point to the point r on the surface Si,,,, knrDDlLis a normalization coefficient calculated from the electroneutrality condition eq 42, uiaW(s)is a surface ion charge density on the ion-accessible surface Sia,/1( TK NLPB NSA (knrNsA; eq 37) TIL (ATIL, A) DDIL
charges at the two positions" -0.3
-1.0
-2.0
-3.0
F,d
-19.47 -19.55
-216.33 -217.19
-865.32 -867.51
-1946.94 -1952.07
-0.153 -0.155 -0.165 (2.65) -0.158 (5.50) -0.207 -1.697 -1.872 -1.986 (2.51) -1.921 (4.75) -2.274 -9.51 1 -10.55 -12.49 (1.31) -10.28 (2.80) -9.218 -15.280 -30.750 -49.63 (0.32) -28.31 (1.80) -26.74
F,d
F:d
0.706' 0.831f 0.738g 0.750h
salt type concentration, M ko-' 1:l 0.1 M 9.741 A
-1.822 -1.789 -1.784 -1.637
0.0 0.063 0.0 0.12
-20.257 -19.881 -19.881 -18.187
0.0
1:l 0.2 M 6.888 8,
-80.97 -79.22 -79.31 -72.76
0.0
2.78 0.0 5.64
31.32' 37.011 32.789 33.35h
1:l 0.5 M 4.355 A
-182.32 -187.92 -178.39 -163.53
0.0 6.216 0.0 12.69
70.56' 83.28, 73.809 75.03h
7.839' 0.696 9.253' 0.0 8.197g 1.41 8.337h
Two charges q of equal sign and magnitude inside a sphere of radius of 5 A; charge positions are (-3.5, 0.0, 0.0 A) and (-2.475, 0.0, 2.475 A). The thickness of the mobile ion shell around the surface of the biomolecule solute was 22 8,. Solvent dielectric constant DO = 80.0, dielectric constant inside the sphere D, = 1.0, salt 1:1, concentration of 0.2 M, Debye constant k~ = 0.1452 A-' or kD-' = 6.89 8, at 25 O C . * E,,I, is the energy of solvation of the biomolecule in kcal/mol, computed by the analytical solution for the TK model and by the numerical IBE method, is the energy of interaction of the biomolecule eq 14, respectively. E,,,tl solute with the mobile ions in kcal/mol, computed from eq 25 for the ion calculated by the analytical solution of the TK model and density plon(r) by the NLPB, NSA, TIL, and DDIL methods, respectively. Fx,F,, and F,are the components of the solvation forces, in (kcal/mol)/A, calculated from eq 23, where footnotes e and g are the forces for the analytical solution of the TK model, for atoms 1 and 2, respectively, and footnotes fand h are the forces of the IBE method for atoms 1 and 2, respectively.
equation. It is seen that the value of the energy E,,,,, calculated by the IBE method, coincides well with the analytical results of the TK model. When the values of the charges of the solute molecule are increased, the condition for linearizing the PB equation is less satisfied, and the values of the energy E,,,,, calculated by the solution of the NLPB equation, are always higher (more negative) than the results of the TK model. The TIL model leads to a reasonable value of the empirical parameter AT]^, i.e., less than the Debye length kD-l, for each case considered. The results of some calculations are shown in Tables I and 11. When the absolute value of the charges of the biomolecule solute increases and, consequently, the electrostatic potential around the solute molecule rises, the value of the parameter ATILdecreases (Table I). The same behavior is observed in Table I1 when the charge of the mobile ions increases. Such behavior is in accordance with the physical meaning of ATIL as a distance between the molecular surface and the surface of the most probable location of the mobile ions around the solute molecule. We fitted the values of the parameters ATILin Tables I and I1 to a simple empirical expression ATIL
TABLE Ik Resulh of Calculations for the Dependence on Ionic Strength for a Pair of Identical Charges in a Spherical Tanford-Kirkwood Model.
= kD-' exp(-0.545qa,,Zi)
(46)
where qaveis the average of the absolute values of the partial charges of the solute molecule and Zi is the valence of the counterions. This TIL model allows us to estimate the energy Ei,, for the direct electrostatic interaction of the solute molecule
1:l 1.0 M 3.076 8, 2:2 0.05 M 6.888 8,
Et,": TK NLPB NSA (knrsA) TIL (ATIL,A) DDIL -1.28 1 -1.544 -1.827 (3.01) -1.572 (6.50) -2.045 -1.697 -1.872 -1.986 (2.51) -1.921 (4.75) -2.274 -2.378 -2.367 -2.362 (2.09) -2.481 (3.00) -2.748 -2.968 -2.772 -2.735 (1.49) -2.976 (2.05) -3.157 -1.697 -2.638 -3.123 (1.31) -2.712 (2.50) -2.304
EisLl,C
W d , d
NLPB NSA TIL DDIL
NLPB NSA
0.43 1 0.6 17 0.702 0.8 15
0.649 0.878
DDIL
0.987
0.579 0.68 1 0.890 0.899
0.730 0.730
0.807 0.854 1.154 1.116
0.784 0.621
1.008 1.035 1.386 0.9 14
0.783 0.576
0.803 1.103 1.247 0.9 14
0.875
0.831
0.840 1.053 1.342 1.470
Two equal charges q = -1 e inside the sphere. Positions of charges are the same as in Table I. The thickness of the mobile ion shell around the surface of the biomolecule solute was 22 A. Same as in Table I (kcal/mol). E,Ul, is the self-electrostatic energy of the mobile ion cloud, eq 27, kcal/mol. A,4,,dLa,is the free energy change due to a change of the configurational entropy of the mobile ion distribution, eq 29, kcal/ mol.
with the mobile ions with an error within 10% relative to the results of the solution of the NLPB equation in a wide range of ionic strength for the model spherical system. The DDIL model describes the ion density profile as a diffuse layer of the ions primarily condensed on the ion-accessiblesurface S,,,.Such a model can be regarded as a suitable trial function for the density functional which in general leads to the NLPB equation. Results of calculations of the energies E,,,, and E,,,, by the DDIL model are presented in Tables I and 11. From these results, we can conclude that theDDIL model gives correct enough results for these energies; the deviations from the NLPB are within 20%. In the main, these deviations indicate that the effective width of the diffuse ion layer depends on the charge density of the biomolecule solute. The deviations increase when the charges of the biomolecule increase (see Table I). The use of the value of the parameter ATILT'instead of k D in eq 45 takes into account the dependence of the width of the diffuse ion layer on the electrostatic properties of the biomolecule solute; AT~L-' was used in the DDIL calculation on decalysine. Another possibility to choose an optimal value of the parameter for the effective width of the diffuse ion layer in the DDIL model is to consider the parameter l c in ~ eq 45 as a variable parameter kDv. The optimal value of kDv can be calculated by minimization of the total free energy, E,,,,(kDV)+ E,,lr(kDV)+ M,,dul, with respect to kDV.To calculate the configurationalfree energy M,,,, according to eq 29, we need the partial ion densities pl(r), which are not provided directly by the DDIL method. We found that the total ion charge density plon(r)of the DDIL method can be decomposed into the partial ion densities p&) by satisfying all three of the following conditions
4862
The Journal of Physical Chemistry, Vol. 97, No. 18. 1993
Vorobjev and Scheraga
TABLE IIk Results of Calculations for Ionized JXlysh&
&(r) /
= min, where min 1 cp:
(47)
salt type concentration. M
Em~v
I
&/Di conformatione
We found that for symmetrical electrolytes this decomposition procedure works fine and the calculated values of the configurational free energy, AAi,d,lDDIL, of the DDIL method agree fairly well with those of the NLPB method (see Table 11). The NSA model based on eqs 32, 37,and 38 leads to fairly good results for the energies EimOl, Eiwl,., and Mi,,,, for a spherical biomolecule solute (see Tables I and 11). The values of the normalization coefficient k,,NSA in eq 37 vary considerably as the bulk density of the mobile ions or the charges of the biomolecule solute and the mobile ions vary; hence, normalization of the SA is essential. The only exception to the good agreement between the results of the NSA and the NLPB models is found in the case of a very highly charged biomolecule solute (see the data for a charge of -3.0 in Table I). In Table I11 are shown the results of the calculations of the solvation energy, E,I,, and energies EimOl, Ei,,, and Mi,,,, for ionized dilysine in the a-helical and extended conformations. We used eq 46 to determine the parameter A T ~ Lfor charged oligipeptides and found that the TIL model with eq 46 gives reasonable values for Eim0, for these biomolecules. It is seen that, for different concentrations I 1 M of the mobile ions of a 1:l electrolyte, all three approximate models (the TIL, the DDIL, and the NSA) give the sum EimOl + Ei,,, of the energies in a good agreement (within 20%) with the results of the NLPB equation. In Table IV are presented the results of the calculations of the solvationenergies and the energies of interactions with the mobile ions for a totally ionized decalysine. As for ionized dilysine, the results of all three approximate methods agree quite well (within 20%) with the results for the solution of the NLPB equation. From Tables I11 and IV, we can see that the variations of the configurational entropy or free energy AAAi,k, of the mobile ion cloud, for the different conformations of the biomolecule solute, are small compared to the variations in Mi,,,, and Mi,,,. Thus, reliable enough estimations of the change of the total free energy of interaction of a biomolecule solute with the mobile ions, as far as the dependence on the biomolecule conformation is concerned, can be made by considering only two terms, Mi,,, and Mi,,,.It should also be noted that the NSA, the DDIL, and the TIL models simulate the differences in the energies (Mi,,, MjEJ between the a-helical and extended conformations of decalysine with a precision within 1 kcal/mol and 25% for the polar and the lowpolar solvents, respectively, compared to the results of the NLPB equation. Thus, the approximate NSA, DDIL, and TIL models can be used to model the dependence of the total electrostatic energy of a biomolecule solute and its environment on the conformation of the biomolecule solute. We now consider the comparison of the dependence of the solvation energies and the energies of interaction of the mobile ions with the solute molecule on the solvent polarity, Le., on t h e value of the dielectric constant DO.The interaction with ions in a very polar aqueous solvent is less important than the polarization interaction in a highly polar aqueous solvent. For water solution, the difference M,,I,~+ between the solvation energies of the a-helical and extended conformations of decalysine is about -50 kcal/mol, while the difference of the energies of interaction with ions, Mimol + Mi,ll, is only about -2.0 kcal/mol. When we consider a less polar solvent, the energy of interaction of a biomolecule solute with ions, EimOl + increases (becomes more negative) quickly, almost inversely proportional to the dielectricconstantof the solvent,DJDo, while the solvationenergy, E,I,, decreases (becomes less negative) slowly, because it is proportional to (1 - (Di/D0)). In the bottom half of Table IV
+
1:l 0.05 M
-182.19
801 1 a-helix 1:l 0.1 M
-182.19
8011 a-helix 1:l 0.20M
-182.19
8011 a-helix 1:l
-1 82.19
0.50M 8011 a-helix 1:l 1.00 M
-182.19
8011 a-helix 1:l 0.20 M
-44.96
8014 a-helix 1:l 0.20M
-43.74
40/4 a-helix 1:l 0.20M
-43.63
NLPB NSA TIL DDIL
Elh,,, NLPB NSA TIL DDIL
-0.900 -1.118 -1.025 -1.121
0.277 0.428 0.504 0.439
0.326 0.574
-1.068 -1.199 -1.286 -1.199
0.368 0.460 0.630 0.471
0.320 0.431
-1.271 -1.381 -1.580 -1.371
0.461 0.554 0.775 0.560
0.357 0.408
-1.520 -1.546 -2.009 -1.528
0.594 0.635 0.985 0.643
0.369 0.311
-1.704 -1.718 -2.355 -1.676
0.703 0.732 1.153 0.731
0.329 0.274
-1.267 -1.382 -1.580 -1.372
0.460 0.554 0.775 0.560
0.339 0.408
-2.987 -3.075 -4.280 -2.970
1.169 1.256 2.093 1.237
0.645 0.685
-1.264
0.459 0.549 0.766 0.568
0.356 0.401
El,",
-1.373 -1.562 -1.396
8014 extended
AAi,dc.,l
NLPB NSA DDIL
0.979
0.677
0.618
0.412
0.314
0.618
0.856
0.599
a Solvation energy and energies E,",",, E,
