J. Phys. Chem. B 2001, 105, 6507-6514
6507
Extending the Applicability of the Nonlinear Poisson-Boltzmann Equation: Multiple Dielectric Constants and Multivalent Ions† W. Rocchia,§ E. Alexov,‡ and B. Honig*,‡ Research Center “E. Piaggio”, Via DiotisalVi, 2 - Pisa, Italy and HHMI and Department of Biochemistry and Molecular Biophysics, Columbia UniVersity, 630 West 168 Street, New York, New York 10032 ReceiVed: February 7, 2001; In Final Form: April 23, 2001
A new version of the DelPhi program, which provides numerical solutions to the nonlinear Poisson-Boltzmann (PB) equation, is reported. The program can divide space into multiple regions containing different dielectric constants and can treat systems containing mixed salt solutions where the valence and concentration of each ion is different. The electrostatic free energy is calculated by decomposing the various energy terms into Coulombic interactions so that that the calculated free energies are independent of the lattice used to solve the PB equation. This, together with algorithms that optimally position polarization charges on the molecular surface, leads to a significant decrease in the dependence of the electrostatic free energy on the resolution of the lattice used to solve the PB equation and, hence, to a remarkable improvement in the precision of the calculated values. The Gauss-Seidel algorithm used in the current version of DelPhi is retained so that the new program retains many of the optimization features of the old one. The program uses dynamic memory allocation and can easily handle systems requiring large grid dimensionssfor example a 3003 system can be conveniently treated on a single SGI R12000 processor. An algorithm that estimates the best relaxation parameter to solve the nonlinear equation for a given system is described, and is implemented in the program at run time. A number of applications of the program are presented.
1. Introduction Numerical solutions to the Poisson-Boltzmann equation have found multiple applications in biology and chemistry.1-7 Many properties associated with molecules in solution are electrostatic in origin and, perhaps surprisingly, the continuum approximations implicit in the PB equation have proved remarkably accurate in treating a variety of such phenomena. Different numerical techniques, including finite difference,8-12 boundary element13,14 and finite-element15-17 methods, have been used to solve the equation. Programs such as DelPhi,11 UHBD,10 MEAD,18 ITPACT,19 and Mainfold Code (MC)16,17 are widely used, attesting the robustness of the continuum approximations that provide the physical basis of these methods. Although most applications of the PB equation have been limited to the linearized form, the nonlinear PB equation is more accurate for highly charged systems such as DNA, or the surfaces of many biological membranes. Solutions to the nonlinear PB equation are also available as is an expression for the electrostatic free energy consistent with this form of the equation.20 Expressions appropriate to mixtures of salts of different valence have also been derived.21 An important feature of the work reported here involves the derivation of an accurate expression for the electrostatic free energy that is also consistent with a number of techniques that have been incorporated into DelPhi to increase the precision of the finite difference solutions. A second feature introduced in this work is the ability to assign different dielectric constants to different regions of space. In a †
Part of the special issue “Bruce Berne Festschrift”. * To whom correspondence should be addressed. Fax: (212) 305-6926. E-mail:
[email protected]. § Research Center “E. Piaggio”, Via Diotisalvi, 2 - Pisa, Italy. ‡ HHMI and Department of Biochemistry and Molecular Biophysics, Columbia University.
biological context, this would allow assigning different dielectric constants to different regions of a membrane22 or to different regions of a protein with varying degrees of flexibility. The various expressions derived in this work have been implemented in the form of computer algorithms that have been incorporated into a new version of the DelPhi program.23 The new version of DelPhi thus provides solutions to the nonlinear PB equation for the general case of charged and polar molecules with variable internal dielectric constants embedded in media which themselves may have different dielectric constants and which may contain variable concentrations of ions of different valence. Such features expand the range of applications that are possible with DelPhi, two of which will be described below. 2. Methods Multiple Dielectric Regions. Finite difference solutions to the PB equation involve solving the equation in a lattice representation.8 The charges and dielectric properties of the system are mapped onto this lattice, the latter being accomplished by discriminating the dielectric constant interior to a molecule from that of the exterior. An important task in this regard is the correct assignment of boundary points. A given boundary point might be located on the interface between an object (molecule) and a solvent, but it might be located at the interface between two or more objects as well. We define an “external” boundary point as one that is in contact with the solvent, whereas “internal” points separate different dielectric regions which are not exposed to the solvent. This distinction is important because only external points can be in contact with mobile ions in solution. In the present implementation of the Finite Difference Poisson-Boltzmann (FDPB) method, space is subdivided in grid points, where the potential is calculated,
10.1021/jp010454y CCC: $20.00 © 2001 American Chemical Society Published on Web 06/05/2001
6508 J. Phys. Chem. B, Vol. 105, No. 28, 2001
Rocchia et al.
and into midpoints (a midpoint is the central point between two adjacent grid points). To make it possible to handle different dielectric regions, each midpoint has to be associated with the dielectric constant of the region in which it is located. This is accomplished, while at the same time limiting memory requirements, by labeling each medium of a different dielectric constant with an integer. A look-up table that is read immediately prior to the iteration procedure correlates each label to the actual dielectric constant value of the midpoint. This avoids the need for storing nonessential information. Free Energies. Given the complexity of the system, it is necessary to partition the total energy into an easily and accurately computable form. The total electrostatic free energy of a system consisting of fixed charges and mobile ions is24,25
∆Gelec )
∫R
(
3
Ffixed φ -
1
)
∑i ∆πi - 2BEDB dV
[ (
) ]
3
1 free F φdV ) 2
∫R
3
1 B ED B dV 2
∑i
1 ∆πi - Fsolvφ 2
(2)
(3)
(4)
In the above equation, the second term is an osmotic pressure and the third is an electrostatic stress.20 It can be shown that in the case of the linearized PB equation the last two terms cancel so that the free energy is given by the first term. This can be used a as criterion of nonlinearity of a system, i.e., in a system with weak nonlinearity these two terms should be small and have similar magnitudes. In the finite-difference method, the system is discretized and the free energy term, 1/2 Ffixedφ, can be rewritten and expanded in the form
1 2
Ffixedφ )
1 2
∑j qjφ(br j)
The Coulombic potential, generated by other fixed charges at the j-th charge position, is calculated in a hypothetical homogeneous media with the dielectric constant i of the region where the charge is located
φcoul(b r j) )
qi
∑ r -b r| i*j 4π |b j
(7)
j
The “corrected reaction field” term, arising from the polarization of the boundary between different media is given by
r j) ) φreact(b
δp
∑p 4π |br - br | 0
and the definition that free (non polarization) charge is the sum of the fixed charge on the molecule and the ionic charge in the solvent, Ffree ) Ffixed + Fsolv, one can rewrite the energy density expression as
1 ∆g ) Ffixedφ 2
(6)
(1)
is the local difference in concentration of the i-th ion compared to the bulk. Equation 2 accounts for the osmotic work of introducing the excess (or deficit) ions into the solution.25,24 Using the well-known expression,26
∫R
φ(b r j) ) φcoul(b r j) + φreact(b r j) + φion(b r j)
0 i
where Ffixed is the fixed charge distribution of the polyelectrolyte, φ is the electrostatic potential, E is electrostatic field, D is the displacement vector in solution, and
eφ -1 ∆πi ) ci exp -zi KBT
duced here avoids this problem. The potential at a given point arises from the direct effect of real charges, from surface polarization charges (the reaction field term) and from mobile ions in solution. Thus, the potential at the position of charge j can be written as
(5)
where the potential is the one generated by all the charges, of any kind, except for the one located at b rj. The standard way used to extract free energies from eq 5 is to multiply the charge that is located at a grid point by the potential at that point. The grid charge is obtained with some method that distributes point charges located at atomic nuclei at neighboring grid points. Because the grid charge is sensitive to this details of this algorithm, and to the grid spacing, this procedure results in a loss of precision in the free energy calculation. An alternative grid-independent procedure intro-
j
(8)
j
where p goes over the locations of polarization surface charges and δp is the polarization charge at position rp on the surface. These polarization charges are calculated using Gauss’s law and the potential map obtained from the finite difference algorithm. Originally, the induced charges are located at grid points. These charges are then repositioned to lie on the true molecular surface at a location defined by the intersection of a normal vector originating at the grid point and the molecular surface.27 The last term accounts for the potential which is generated by the mobile ions in solution. Under the Boltzmann approximation, these charges are a function of the potential, obtained by solving the PB equation, and are screened by the polarizable solvent, whose dielectric constant is ) solv
r j) ) φsolv(b
h3dFion k
∑k 4π
rj 0 solv|b
-b r j|
(9)
where h is the grid spacing, δFion k is the net (excess) ion charge density at the k-th grid point in solution, and k goes over the grid points in solution. The procedure outlined above makes it possible to use real charges rather than their grid representation in calculating the free energy of the system. Once the surface and excess ion charges are calculated, then the potential at the position of all fixed charges can be obtained from eqs 7-9 and the free energy is calculated from eqs 4-5. This leads to major improvements in precision because the calculation of free energies involves only using Coulomb’s law for real, polarization and mobile charges. This also eliminates the so-called grid energy, which is an artificial energy term that comes from the partitioning of point charges at neighboring lattice points. The only drawback of this procedure is that the evaluation of the contribution of mobile ions can be time-consuming. An option has been incorporated in DelPhi which turns this term off, for example for cases where salt effects are expected to be small or nonexistent. In addition, there is also an option to evaluate the free energy using the grid-dependent method. For systems with several dielectric regions, the reference medium is vacuum ( ) 1). To be treat such a system consistently, one has to account for the transfer free energy of a charge from vacuum to the medium with the dielectric
Nonlinear Poisson-Boltzmann Equation
J. Phys. Chem. B, Vol. 105, No. 28, 2001 6509
constant, r, in which it is located. This free energy corresponds to the interaction between the charge and the volumetric polarization charge of the medium in question. We have termed this contribution a “self-reaction energy”. Its value is always negative because transferring a charge from vacuum to some other medium is always a favorable process. The polarization charge can now be written in the form26
( ) ( ) ( )
Fpol ) -∇ B ‚P B)∇ B‚
1 1 D B -1 D B) -1 ∇ B ‚D B - 2∇ B r ) r r - 1-
r
0E 1 free B F ∇ B (10) r r r
Multivalent Salts. The charge density of mobile ions in solution can be described as21,12,28 n
r )e Fsolv(b)
n
r i ) e∑Cbi zi exp(-ziΦ(b)) r ∑1 Ci(b)z 1
(11)
where Cbi ) Ci(∞) is the bulk concentration of the i-th ion type and Zi is its valence. Substituting eq 11 into the PB equation results in
∇ B ‚[(b)∇ r B φ(b)] r ) -
[
1 fixed F (b) r +e 0
(
n
ezi
)]
r ∑1 Cbi zi exp - kTφ(b)
(12)
For symmetric salts, which have the same concentration of negative and positive ions of the same valence, one can define effective Debye-Huckel parameters for each of the ions. In this case, eq 12 could be rewritten as a sum of hyperbolic sin functions, whereas in general this is not possible. DelPhi treats different combinations of ions; i.e., (1:1), (2:2), (1:2), (2:1) salts, where the numbers in the brackets indicate the ion’s valence. The contribution of each ion to the total free energy can be calculated. In this way, it is possible to consider for example, the interplay between mono and divalent ions in stabilizing nucleic acid conformations. Convergence Issues for the Nonlinear Poisson-Boltzmann Equation. The solution to the linearized PB implemented in DelPhi converges very quickly due in part to its over-relaxation procedure.11 The efficiency of the procedure is due to the fact that the Gauss-Seidel method provides an analytical prescription for the best value of the relaxation parameter given the variables that describe the system. For the nonlinear PB equation, this is no longer possible. First, as opposed to equations involving linear algebraic systems, convergence is no longer independent of the initial value. This can be seen from the following onedimensional example. Consider the nonlinear equation
φ ) F(φ) ) tφ + q - k sinh(φ)
φ
) F(φ ) (n)
Figure 1), then the system diverges. The threshold value in this one-dimensional case is φ* ) FF(φ*). On the basis of the assumption that the same behavior holds true for the multidimensional case, the nonlinear solution we have implemented starts with linear iterations and then gradually adds nonlinearity. This is accomplished by adding the nonlinear term multiplied by 0.05 to the initially linear equation. Thus, after 20 blocks of 10 iterations each, the algorithm begins to solve the full nonlinear equation. In this way, the “initial” guess for the nonlinear iteration process is “pushed” toward the “convergent basin”. A second important difference between the linear and nonlinear case is the value of relaxation parameter used in the Gauss-Seidel algorithm. It was shown11 that over-relaxation is the best strategy in the case of the linear PB equation. In contrast, solution of the nonlinear equation requires underrelaxation,11 as will be shown below. Here, we present a qualitative approach to the problem starting from the linear case. Given a linear system of algebraic equations in the matrix form Aφ + Q ) 0, its Jacobi form is
φ ) Tφ + Q
(14)
As can be seen in Figure 1, the choice of different initial values can result in divergence in one case and convergence in the other. Thus, if the initial results in the next step being closer to the correct solution (solid line in Figure 1), the process converges. In contrast, if the second step turns out to be farther from the correct solution than the first one (dashed curve in
(15)
where T ) A + I, and I is the unit matrix. It can be shown that for the PBE the matrix T is consistently ordered,11 i.e., it has the form
[
0 M1 M2 0
(13)
where t, q, and k are parameters. The iteration that arises when one tries to solve this equation in a Jacobi scheme is (n+1)
Figure 1. Schematic representation of an iterative procedure; the bold line represent iterations that converge, while the broken line shows iterations that diverges.
]
(16)
where M1 and M2 are two sub matrices. Thus, T can be thought of as the sum of a lower triangular matrix E and an upper triangular matrix F, the properties of which will be used below. Equation 16 can be rewritten in a form suitable for the Gauss-Seidel technique
Bφ ) Fφ + Q where B ) I - E.
(17)
6510 J. Phys. Chem. B, Vol. 105, No. 28, 2001
Rocchia et al.
Solving eq 17 with respect to φ, we obtain
φ ) B-1Fφ + B-1Q ) Jφ + B-1Q
(18)
A major advantage of the Gauss-Seidel method is the ability to condense two Jacobi iterations in one step thus, accelerating convergence. In addition, further optimization is possible.11 It can be shown that the speed of convergence depends on the spectral radius of J, which is the highest eigenvalue of J in modulus. A transformation on the system is carried out which tends to decrease the spectral radius of J preserving the asymptotic value for the succession φ(n) ) Jφ(n-1) + B-1Q. This transformation technique introduces the relaxation parameter ω
B)
I -E ω
(19)
where 0 < ω < 2. If ω < 1, the technique is called under-relaxation, whereas if ω > 1 it is called over-relaxation. The case of ω ) 1 corresponds to the original Gauss-Seidel form.29 This transformation does not modify the symmetry properties of the matrices that ensure the convergence of Gauss-Seidel and provides the optimal convergence rate for the relaxation parameter
ωop )
2
(20)
1 + x1 - F2
where F is the spectral radius of J before the transformation. As it can be seen, if F is real and less than 1, then ωop is greater than 1. Non linearity makes the system much more complicated. In our approach the nonlinearity is treated as a perturbation of the linear system, and the expression for the error at the n-th step is
∆(n) ) φ(n) - φ(∞) ) J∆(n-1) - B-1(u (n-1) - u (∞))
(21)
where φ(∞) is the final (correct) potential and u (n) is a vector whose j-th component is
κ2h2 κh + 2 2
∑i
(n) [sinh(φ(n) j ) - φj ]
(22)
i
Thus, the j-th component (u (n) - u (∞))j can be rewritten as follows
κ2h2 κh + 2 2
∆(n-1) × j
∑i i
[
]
(n-1) exp(∆(n-1) ) - 1 exp(φ(∞) ) j j ) - exp(φj ‚ - 1 (23) 2 ∆(n-1) j
If we assume that the values of the vector of errors ∆(n-1) are small enough, the eq 23 can be approximated as
κ2h2 κh + 2 2
∆(n-1) [sinh(φj) - 1] ) j
∑i i
〈
κ2h2 κh + 2 2
〉
[sinh(φj) - 1] ∆(n-1) ≡ (s∆(n-1))j (24) j
∑i i
Figure 2. Magnitude of the relaxation parameter “ω” for the nonlinear Poisson-Boltzmann solver as a function of net charge of the protein (file 1lyz.pdb43 from the Protein Data Bank44).
Equation 24 introduces a new quantity “s”, which is an average throughout the space of the expression between the large brackets. It can be seen that “s” depends on several parameters but mostly on the local potential through the hyperbolic sin function. Under these approximations, the error vector can be rewritten as
∆(n) ) (J - sB-1)∆(n-1)
(25)
Attempts to derive a formula providing the optimal relaxation parameter that minimizes the spectral radius of the new matrix system were unsuccessful because “s” depends on the value of the potential throughout the space, which is not known a priori. The “instability” resulting from the nonlinearity is greatest when the absolute value of “s” is high. This occurs when the electrostatic potential is large in solution, where the nonlinear equation has to be solved. As is expected, this occurs when the system is highly charged and when the charges are close to the molecular surface. For these reasons, we consider the convergence problems of a “difficult” system which has a high net charge and where the individual charges are close to the surface. The effect of the net charge on the value of the optimal relaxation parameter for such a case is shown in Figure 2 , a protein of about thousand heavy atoms and average radius of about 30 Å. The charges were assigned with the Parse set30 but the partial charge of the backbone nitrogens was artificially varied to test the relaxation at different net charges of the protein. In our attempt to study the behavior of the system, we need to find the spectral radius of J - sB-1 before the relaxation takes place, i.e., when ω ) 1. This means that we have to estimate the largest eigenvalue, in modulus, for J - s(I - E)-1. Given the structure of the matrices we are considering, we still are dealing with a system ruled by a consistently ordered matrix that can be transformed in the following way J -s(I + E ˜), where E ˜ is modified lower triangular matrix. A plausible estimate for the spectral radius of this matrix is F* ) |s(1 + xF) -F|, F is the spectral radius before relaxation. The above expression can be greater than one and in that case the iteration process diverges. Thus, will extend the theorem given in Ref.29 to the case of a transformed matrix with complex eigenvalues. According to control theory31 such an iterative process will converge if the real part of the eigenvalues is less than one. The descending qualitative behavior of ωop versus “s” can be
Nonlinear Poisson-Boltzmann Equation
Figure 3. Relaxation parameter ω as a function of the variable s.
seen in Figure 3. The value corresponding to s ) 0 is the value for ωop obtained in the linear case, which in Figure 3 has been chosen to be 1.58. To summarize this section, the larger the nonlinearity, the lower the relaxation parameter has to be in order to ensure convergence. However, a low relaxation parameter dramatically increases the number of steps required for convergence. The new version of Delphi has a heuristic algorithm which estimates the optimal parameter at run time, trying to keep it just below the threshold of divergence. The first nonlinear iteration is carried out using the relaxation parameter that is estimated from the linear equation. As described above, the value of the relaxation parameter for the linear equation is always larger (due to the overrelaxation procedure), than that appropriate to the nonlinear equation. Thus, the program makes its first trial using the relaxation parameter estimated with the linear algorithm. The maximum change in the rmsd of the potential calculated at the current iteration is then compared to the corresponding values of the three previous steps. If there is evidence that the system is leaving the “convergence basin” the procedure decreases the relaxation parameter trying to find a better one. Given the peculiarity of this procedure, and despite the fact that it has been tested successfully on many different systems, it is still possible that in some particular case it will be unable to react quickly enough to the system evolution to give the correct estimate. In this case, the user has the option of manually assigning a particular relaxation parameter so as to achieve improved convergence. 3. Results Test on an Idealized System. Calculations were carried out on a system consisting of a sphere of dielectric constant 1 and radius “a”, with a point dipole (P0 ) 4 D) placed at its center. The sphere was surrounded by a shell of thickness “b-a” with a dielectric constant 2. The rest of space was considered as an infinite medium with a dielectric constant 3 and Debye-Huckel parameter κ. There is an analytic solution to the PB equation for this problem (see the Appendix). The general formalism can be found in many textbooks on electrostatics26 and a similar example has been solved in ref 32. The numerical solution for the electrostatic potential was tested against the analytical solution. In all cases, the calculated potentials matched the analytical solution. The following three examples were used to test the numerical procedures. (a) 1 ) 4, a ) 4 Å; 2 ) 20, b ) 6 Å and 3 ) 80. At a point 14.1 Å away from the center of the dipole and a radius
J. Phys. Chem. B, Vol. 105, No. 28, 2001 6511 vector making a 45° angle with the direction of the dipole, the analytical solution is 0.1693 [kT/e]. The calculated potential is 0.1676 [kT/e] at grid resolution of 1.3 [grids/Å]. The error is 1%. (b)the above example, but 1 ) 1. At the same point, the analytical solution is 0.1784 [kT/e] and the calculated value is 0.1793 [kT/e], which corresponds to a 0.5% error at the same grid resolution. (c) same as example (a), but including ions in the aqueous phase (I ) 0.1 M). At the same point the analytical solution is 0.1059 [kT/e], and the numerical value is 0.104 [kT/e]. The corresponding error is 1.8% at a grid resolution of 1.3 [grids/ Å]. Using Multiple Dielectrics to Study Electron Transfer. Electron transfer in proteins is strongly affected by relaxation processes where nuclei rearrange to accommodate a charge distribution induced by the transfer of the electron. The effect of these motions can be described in terms of a dielectric constant, but given the known heterogeneity in protein flexibility, the appropriate value for might be quite different for different regions of a protein. As an example of how this issue might come into play, we discuss electron transfer in photosynthetic reaction center (RC) proteins. RCs transfers electron(s) from the extracellular to the intracellular side of a membrane through a cascade of cofactors. The last step in this sequence involves transfer from QA to QB, where both species, Q, correspond to ubiquinone 10. The simplest explanation for how a protein might drive an electron from QA to QB is that the electrostatic potential at QB is more positive than at QA. However, several attempts to compute the free energy of the reaction as well as 3D potential distributions33,34 have found that the reaction is uphill in free energy, that is, the potential is more positive at the QA site. These studies all assumed a static protein which had a uniform low dielectric constant. Recently, the reaction was studied with the MCCE method,6 that includes the effects of conformational changes in the calculation of proton equilibria and redox effects, and it was found that reaction can occur only if the QB pocket undergoes significant rearrangement upon electron transfer. It was demonstrated that many polar side chains and water molecules reorient as a result of the electron transfer and that most of them are in the QB pocket. Thus, it appears that the RC protein is designed so that the QA pocket is much more rigid then the QB pocket, resulting in different polarizabilities, i.e., different dielectric constants, for the two sites. To estimate the consequences of this result within the continuum approximation, DelPhi calculations were carried out with the QB pocket treated as a region of dielectric constant 20, whereas the rest of the protein was assigned a dielectric constant of 4. On the basis of the analysis of Alexov and Gunner,6 residues L201-L226, M232-M236, H119-H124, H228-H230, M14-M20 were use to define the high dielectric region. All titratable groups were assumed to be charged with exception of all Tyr and His, and GluL104 and GluL212.6,33,34 Charges and radii were taken from Parse data set.30 Hydrogen atom positions were generated so as to generate optimal hydrogen bond geometries.35,6 To demonstrate the effect of variable dielectric constant, the electrostatic potential distribution was first calculated with the entire protein assigned a dielectric constant of p ) 4. As can seen in Figure 4a, the potential at the QB site is far more negative that at QA position. The reason for this is that QB pocket is rich in acidic groups that generate a strong negative potential at the QB position33,34 thus opposing electron transfer from QA to QB.
6512 J. Phys. Chem. B, Vol. 105, No. 28, 2001
Rocchia et al.
Figure 5. Electrostatic energy of association of the barnase-barstar complex as a function of salt concentration. Experimental data taken from Ref.46 b- experimental data, 0 - linear PB and 2 - nonlinear PB.
added to polar heavy atoms in the complex and in the individual subunits, using the same technique that was mentioned above.35 Charges and radii were taken from the Parse parameter set.30 The electrostatic contribution to binding was calculated as
∆∆Gelec (binding) ) ∆Gelec (complex) ∆Gelec (barnase) - ∆Gelec (barstar) (26)
Figure 4. Electrostatic potentials through a plane in a RC protein. The plane is determined by the position of the Mg ion and the O2 atoms of the quinines. The peptide backbone is shown in green; QA and QB are black. Red corresponds to negative potentials, whereas blue corresponds to positive potentials. The figure is produced with GRASP.45 (a) Homogeneous protein with dielectric constant, ) 4. (b) Inhomogeneous protein. Dielecric constant of protein is 4 except for the QB pocket which has a dielectric constant, ) 20.
Figure 4b shows the potential distribution for the case where the QB pocked has a dielectric constant of 20, whereas the rest of protein has a dielectric constant of 4. It can be seen (Figure 4b) that the high dielectric effectively reduces the negative potential at QB site and in this way facilitates the reaction. The potential distribution away from QB pocket is essentially unaffected. Another effect that we have not considered here directly, but which may be of greater importance, is that the rearrangement of polar groups can stabilize the electron on the QB site through direct interactions with the quinone. This effect can also be modeled by assigning a higher dielectric constant to QB than to QA and can also be treated with the new version of DelPhi. Salt Efects on Protein-Protein Binding Free Energies. Salt effects on binding have been treated extensively with the PB equation and, for example, remarkable agreement with experiment has been obtained for the binding of proteins and ligands to DNA.36-38 In the following, we compare theory and experiment for the salt dependence of the binding of barstar to barnase,39 as system for which detailed experimental results are available.40 DelPhi calculations were carried out on the barnasebarstar complex; file 1b2s.pdb41 was used. Hydrogen atoms were
Three different runs were carried out for each ionic strength (I) using both the linear and nonlinear PB solver in DelPhi. As can be seen in Figure 5, the observed decrease in binding affinity at higher salt concentrations is reproduced by the calculations, although the calculated magnitude of the effect is somewhat smaller than observed experimentally. The calculated energies are adjusted at I ) 0 to match the experimental value. The linear and nonlinear equations yield very similar results, in part because of the relatively low net charge on the individual proteins. The calculations were carried out at on a 1293 grid without focusing. One run was carried out for each molecule resulting in grid resolutions of 1.0, 0.87, and 0.74 grids/Å for barstar, barnase and the complex, respectively. The proteins were automatically centered by DelPhi in the center of the cube and no attempt was made to keep them at the same grid position in different runs. To test the sensitivity of the algorithm to the position of the proteins, their centers within the grid were systematically changed by a non integer number. For example, by default DelPhi centers the objects at the center of the cube. We manipulated the program to center the objects at a 1.5 grid offset from the geometrical center along the x, y, and z directions. The resulting variation of the free energy was 0.01% demonstrating the robustness of the energy calculations. To test the effect of the grid resolution, the electrostatic free energy of barnase was calculated at different grid spacing. The calculations were carried out proceeding from a quite low resolution (0.43 grids/Å) with 653 lattice points to 2.16 grids/Å with 3193 lattice points (Figure 6). It can be seen that at grid resolutions of about 1grid/Å and higher, the energy reaches a constant value. Program Performance. In general, an expert user of a program such as DelPhi needs to explore the sensitivity of the various parameters in the context of the system being studied. In this section, we summarize the performance of the current
Nonlinear Poisson-Boltzmann Equation
Figure 6. Electrostatic free energy of barnase as a function of grid resolution.
version of DelPhi under a limited set of conditions. The linear PB solver in DelPhi is extremely fast due to the overrelaxation procedure it uses which also makes use of the symmetry of the linear equation (allowing the number of cells grids calculated at given iteration to be reduced by a factor of 2).11 The new nonlinear solver uses the output of the linear algorithm as a starting point and, in this way exploits its computational efficiency. In cases where the nonlinear solution is close to the linear one, this provides a significant advantage. Recently, Holst et al.16 have described a quasi-Newtonian algorithm to solve the nonlinear PB equation whose convergence properties are superior to the Gauss-Seidel method used in the work. We have not adopted this method so as to retain maximum consistency between our linear and nonlinear solvers. A hybrid approach which combines the advantages of the two, as suggested by a referee, will be considered in future work. The convergence algorithm was tested on an Octane 2 SGI computer, with a 400 MHz MIPS R12000 CPU and 2.2 GB RAM. At low grid dimensions 653 and an input file of several thousand atoms, convergence was achieved in 2.06 s. For a system of more than 60 000 atoms and a net charge of 130 electron units, convergence was achieved in 5.7 s, and the total grid independent energy was calculated in an additional 28.1 s. Convergence was defined as a change in the rmsd of the potential between two successive iterations of less than 0.001. The largest grid size currently allowable in DelPhi is 319.3 In this case, for the smaller system, convergence (accuracy < 0.001) was achieved in 593.1 s and the grid independent energy was calculated in an additional 55.2 s. This limit can easily be extended, it is the calculation of salt contributions to the free energy that is time-consuming. 4. Concluding Remarks We have described a number of features of a new version of the DelPhi program, that provides numerical solutions to the PB equation. These include the ability to assign different dielectric constants to different regions of space, free energy calculations for mixed salt solutions, and a new numerical approach to solving the nonlinear equation. The program is extremely fast and can handle finite difference lattices of extremely high dimension. There are many PB solvers currently available and the technology continues to improve.17,16,42 Our own goals have been to extend the capabilities of the DelPhi program so that it
J. Phys. Chem. B, Vol. 105, No. 28, 2001 6513 can handle increasingly more complex systems. From the perspective of physical chemistry, it is quite remarkable how a classical theory has seen a resurgence due to the application of numerical methods and the availability of fast computers. Although the PB equation had been used to analyze molecules for many years, its full power was only realized when numerical methods made it possible describe molecules in atomic detail. Because the equation has no analytical solution for irregularly shaped objects, it had been necessary to use simplified descriptions of molecules, e.g. spheres, to solve the PB equation. Failures of the approach were frequently attributed to the limitations of the continuum approximation rather than to the extreme geometric simplifications that were used. We now know that the opposite is frequently the case. The availability of programs such as DelPhi make it possible for researchers to explore the approximations inherent in the models used to treat electrostatic phenomena in quite different areas of science. Acknowledgment. This work was supported by Grant No. GM 41137-13 from NIH and Grant No. DBI - 9904841 from NSF. Appendix Electrostatic Potential of an Ideal Dipole Immersed into Two Concentric Spheres. Consider a system composed of an ideal dipole with a dipole moment P0 placed at the center of a sphere of radius “a”. The dielectric constant of the sphere is 1. The space taken up by the sphere is labeled as region I. The sphere is surrounded by a shell of thickness “b - a”, where “b” is the distance between the center of the sphere and the outer surface of the layer. The layer has a dielectric constant 2 and the space occupied by it is labeled as II. The system is then immersed into a bulk aqueous medium that has a dielectric constant 3 and salt concentration corresponding to the DebyeHuckel parameter κ. The water phase is called region III. Thus, one can write the potentials in each of the regions as26,32
P0 cos ϑ
(a) region I:
Ψ(r,ϑ) )
(b) region II:
Ψ(r,ϑ) ) (Cr + D/r2) cos ϑ
(c) region III:
Ψ(r,ϑ) ) F
1r2
+ Ar cos ϑ
x
exp(-κr) 2 (κr + 1) cos ϑ 2 πκ r2
These equations have unknown constants that can be determined from the boundary conditions. Thus, it can be shown that the electrostatic potential outside the sphere (in the aqueous phase) is given by the expression
Ψ(r,ϑ) )
3P02 b3 Ω cos ϑ (κr + 1) exp(-κ(r - b)) 2 3 (1 - 2) a Ξ r
where
b3 22 + 1 +2 a3 1 - 2 Ω)1+ b3 22 + 1 1- 3 a 1 - 2
6514 J. Phys. Chem. B, Vol. 105, No. 28, 2001
[
Ξ)
]
b3 22 + 1 a3 1 - 2 - 3[2 + 2κb + (κb)2]. b3 22 + 1 1- 3 a 1 - 2
2(κb + 1) 2 +
References and Notes (1) Gilson, M.; Honig, B. Nature 1987, 330, 84-84. (2) Honig, B.; Nicholls, A. Science 1995, 268, 1144-1149. (3) Ben-Tal, N.; Honig, B.; Miller, C.; McLaughlin, S. Biophys. J. 1997, 73, 1717-1727. (4) Murray, D.; Hermida-Matsumoto, L.; Buser, C.; Tsang, J.; Sigal, C.; Ben-Tal, N.; Honig, B.; Resh, M.; McLaughlin, S. Biochemistry 1998, 37, 2145-2159. (5) McCammon, J. A. Curr. Opin. Struct. Biol. 1998, 8, 245-249. (6) Alexov, E.; Gunner, M. Biochemistry 1999, 38, 8253-8270. (7) Sheinerman, F.; Norel, R.; Honig, B. Curr. Op. Struct. Biol. 2000, 10, 153-159. (8) Warwicker, J.; Watson, H. C. J. Mol. Biol. 1982, 157, 671. (9) Gilson, M.; Sharp, K. A.; Honig, B. J. Comput. Chem. 1987, 9, 327-335. (10) Davis, M. E.; McCammon, J. A. J. Comput. Chem. 1989, 10, 386. (11) Nicholls, A.; Honig, B. J. Comput. Chem. 1991, 12, 435-445. (12) Sharp, K. A. Biopolymers 1995, 36, 227-243. (13) Zauhar, R.; Morgan, R. J. J. Mol. Biol. 1985, 186, 815-820. (14) Bharadwaj, R.; Windemuth, A.; Sridharan, S.; Honig, B.; Nicholls, A. J. Comp. Chem. 1995, 16, 898-913. (15) You, T. J.; Harvey, S. C. J. Comput. Chem. 1993, 14, 484-501. (16) Holst, M.; Baker, N.; Wang, M. J. Comput. Chem. 2000, 21, 13191342. (17) Backer, N.; Holst, M.; Wang, F. J. Comput. Chem. 2000, 21, 13431352. (18) An Object-Oriented Programming Suite for Electrostatic Effects in Biological Molecules; Bashford, D., Ed.; Springer: Berlin, 1997; p 223. (19) Cortis, C.; Friesner, R. J. Comput. Chem. 1997, 18, 1591-1608. (20) Sharp, K. A.; Honig, B. J. Phys. Chem. 1990, 94, 7684-92. (21) Chen, S. W.; Honig, B. J. Phys. Chem. 1997, 101, 9113-9118. (22) Karshikoff, A.; Spassov, V.; Cowan, S.; Ladenstein, R.; Schirmen, T. J. Mol. Biol. 1994, 240, 372-384.
Rocchia et al. (23) Nicholls, A.; Sharp, K. A.; Honig, B. DelPhi; 3.0 ed.; Department of Biochemistry and Molecular Biophysics, Columbia University: New York, 1990. (24) Sharp, K.; Honig, B. J. Phys. Chem. 1990, 94, 7684-7692. (25) Sharp, K. Biopolymers 1995, 36, 227-243. (26) Jackson, J. D. Classical Electrodynamics; John Wiley and Sons: New York, 1962. (27) Rocchia, W.; Sridharan, S.; Nicholls, A.; Alexov, E.; Chiabrera, A.; Honig, B. J. Comput. Chem. 2001, submitted. (28) Pack, G.; Wong, L.; Lamm, G. Biopolymers 1999, 49, 575-590. (29) Stoer, J.; Burlirsch, R. Introduction to Numerical Analysis; SpringerVerlag: New York:, 1980. (30) Sitkoff, D.; Sharp, K. A.; Honig, B. J. Phys. Chem. 1994, 98, 19781988. (31) Marro, M. Teoria dei Sistemi e del Contollo; Zanichelli, Ed. Bologna, 1989. (32) Alexov, E.; Michonova, E. J. Macromol. Sci.-Phys. 1994, B33, 21-35. (33) Gunner, M.; Nicholls, A.; Honig, B. J. Phys. Chem. 1996, 100, 4277-4291. (34) Lancaster, R.; Michel, H.; Honig, B.; Gunner, M. R. Biophys. J. 1996, 70, 2469-2492. (35) Alexov, E.; Gunner, M. Biophys. J. 1997, 74, 2075-2093. (36) Misra, V. K.; Sharp, K. A.; Friedman, R. A.; Honig, B. J. Mol. Biol. 1994, 238, 245-63. (37) Misra, V. K.; Hecht, J. L.; Sharp, K. A.; Friedman, R. A.; Honig, B. J. Mol. Biol. 1994, 238, 264-280. (38) Misra, V.; Honig, B. Proc. Natl. Acad. Sci. U.S.A. 1995, 92, 46914695. (39) Hartley, R. W. Biochemistry 1993, 32, 5978-5984. (40) Schreiber, G.; Ferst, A. Biochemistry 1993, 32, 5145-5150. (41) Buckle, A.; Schereiber, S.; Ferst, A. Biochemistry 1994, 33, 88788883. (42) Grant, A.; Pickup, B.; Nicholls, A. J. Com. Chem. 2001, 22, 608640. (43) Diamond, R. J. Mol. Biol. 1974, 82, 371. (44) Bernstein, F. C.; Koetzle, T. F.; Williams, G. J.; Meyer, E. F.; Brice, M. D.; Rodgers, J. R.; Kennard, O.; Shimanouchi, T.; Tasumi, M. J. Mol. Biol. 1977, 112, 535-542. (45) Nicholls, A.; Sharp, K. A.; Honig, B. Proteins: Struct. Funct. Genet. 1991, 11, 281-296. (46) Schreiber, G.; Fersht, A. R. Biochemistry 1993, 32, 5145-5150.
