Next: PBE Data Structures, Previous: Poisson-Boltzmann Electrostatics, Up: Poisson-Boltzmann Electrostatics

One of the most successful implicit models for the treatment of electrostatic
effects is the Poisson-Boltzmann equation^{1}
which is given below

grad (epsilon(x) grad phi(x)) - kappa(x)^2*sinh(phi(x)e/kT) = -4 pi k_e rho(x) where

phi is the electrostatic potential, rho is the charge density, epsilon(x) is the dielectric constant in units where the vacuum dielectric is exactly 1, k_e is the electrostatic force constant, e is the charge on the electron, k is the Boltzmann constant, T is the temperature, and kappa is a dielectric independent Debye-Huckel parameter, kappa = sqrt(epsilon) * kappa

The equation can be solved over a volume enclosing a molecule of
interest using finite difference methods^{2} ^{3} ^{4}
appropriate
for the solution of boundary value problems. In these methods, a grid
is laid down over the volume, and each grid point is assigned a value
for the ionic strength and charge density. The lines between each grid point are
assigned values for dielectric constant. Boundary values for the potential are set
according to a variety of models and then,
interior values for the electrostatic potential are calculated.
Typically, the equation is linearized by substituting x for sinh(x) in
the above equation.

The implementation of the finite difference solution of the
Poisson-Boltzmann equation in CONGEN is written in C and uses dynamic
storage allocation throughout so any size grid can be accommodated. The
process of a potential calculation begins with the determination of the
grid position. This can be centered on the spatial origin, the center of a
rectangular box enclosing the molecule, or the geometric center of the
molecule. Next, the dielectric constant stored on all of the grid edges
is set to the solvent value, and the Debye-Huckel parameters are
likewise set. Next, the program loops over all atoms. Grid points of
the protein excluded space are assigned grid charges as described above,
their dielectric constants are set to interior values, and the
Debye-Huckel parameters are set to zero. The van der Waals, accessible,
or molecular surface
can be used to delimit the interior. Upon the direction of the
user, the program can set the boundary using one of three possible
rules. The first possibility is to set it to zero. This is very fast,
but is only appropriate when the boundary is very far from the atoms in
the system. The second possibility, which is denoted as the system
Debye rule, is to set it to the potential caused by a single charge
equal to the total charge of the system and whose radius is equal to
that of a sphere whose volume equals the solvent excluded volume of the
system. This method is also efficient, but it can capture some of the
electrostatic screening due to the solvent if there is a thick solvent
boundary around the system. The final boundary setting method, which is
denoted as the atomic Debye rule, sets the boundary points to the sum of
the potentials due to all of the atoms scaled by the Debye-Huckel
rule,^{5} The atomic Debye rule is the most
accurate.

phi = k_e q_i e^(-kappa r_s) / (\e_s r (1 + \kappa R))

where r is the distance from boundary point to the atom, R is the radius of the atom plus the water radius plus the length of the ion exclusion (Stern) layer and r_s is r-R.

This last option does a more detailed calculation of the boundary, but is more expensive computationally.

The potential is calculated using Gauss-Seidel iteration with odd-even
ordering of processing grid points.^{6} When the linearized form of the equation
is solved, overrelaxation is used. The default spectral radius used for the
overrelaxation parameter is given by

rho = 1/3 (cos(pi/n_x) +cos(pi/n_y) + cos(pi/n_z))

where n_x, n_y, and n_z are the grid dimensions in the x, y, and z directions.

[1] Stephen C. Harvey,
“Treatment of Electrostatic Effects in Macromolecular Modeling”,
*Proteins: Struct. Funct. Gen.* **5**, 78-92 (1989)

[2]
J. Warwicker and H. C. Watson, “Calculation of the Electric Potential in the Active
Site Cleft Due to alpha-Helix Dipoles”, *J. Mol. Biol.* **157**,
671-679, (1982).

[3]
Malcolm E. Davis and J. Andrew McCammon,
“Solving the Finite Difference Linearized
Poisson-Boltzmann Equation: A Comparison of
Relaxation and Conjugate Gradient Methods”,
*J. Comput. Chem*, **10**,
386-391, (1989).

[4]
Anthony Nicholls and Barry Honig,
“A Rapid Finite Difference Algorithm, Utilizing
Successive Over-Relaxation to Solve the
Poisson-Boltzmann Equation”,
*J. Comput. Chem* **12**,
435-445, (1991).

[5] Malcolm E. Davis and J. Andrew McCammon,
“Electrostatics in Biomolecular Structure and Dynamics”,
*Chem. Rev.* **90** 509-521 (1990).

[6]
W. H. Press and B. P. Flannery and S. A. Teukolsky and W. T. Vetterling,
*Numerical Recipes: The Art of Scientific Computing*,
Cambridge University Press, Cambridge, 1986,
pp. 652-659.