Electrical potential isosurfaces (as predicted by our Poisson-Boltzmann solver) for a poliovirus pentamer at -0.5 volts (blue) and +0.5 volts (red). 152 million grid points were used to compute the potential around 59,465 atoms.

Poisson-Boltzmann ^®: a 3-D numerical simulator

Electrostatic interactions play a crucial role in determining the structure and behavior of proteins and more complex structures (e.g. enzymes and viruses). For example, DNA has an overall negative charge due to phosphate groups. Biomolecules reside in an aqueous electrolyte, which affects their conformation and function due to screening and dielectric effects. Aqueous physiological media contain many mobile ions (e.g. Cl-, Na+, K+, Mg++, and Ca++), which redistribute and screen the Coulomb potential of the fixed charges on the macromolecules by creating counter-ions layers. Orientable or polarizable molecules of the host medium decrease the electrical forces that determine the structure and function of macromolecules and other nano-structures. Accurate calculation of the electrostatic potential can enhance our understanding of the behavior and structure of macromolecules. One can use such calculations to estimate the solvation energy, find pKa values, and titration curves for proteins. Also, one can calculate the electrostatic forces between biomolecules for use in molecular dynamics. In carrying out molecular mechanics computations, it is important to account for the channeling of the electric field along a macromolecule due to the dielectric constant contrast between the aqueous medium and the interior of a protein, an effect not accounted for in 1/r Coulomb computations as in the CHARMM force field.

The rigorous statistical mechanical theory of Coulomb systems is fraught with challenges due to the long range of the interparticle potential. The starting point is the determination of the quantum energy states of the electron­nucleus system. From the Born­Oppenheimer approximation, one arrives at an effective N-body potential for the nuclei when the electronic system can be assumed to be in the ground state. Thus, the problem is reduced to the statistical mechanics of many nuclei evolving in the Born­Oppenheimer potential. As a practical approach, phenomenological formulations based on heuristic arguments following from classical continuum electrostatics and equilibrium statistical mechanics are usually used to compute the electrostatic potential around a molecule immersed in an electrolyte. In such an approach, the host medium consists of water molecules and mobile ions accounted for by using macroscopic continuum densities.

The Poisson-Boltzmann (PB) equation has been traditionally used to find the electrostatic potential around a macromolecule. A classic PB model ignores the volume of ions in the medium. Therefore, the PB equation is only valid for dilute ionic solutions beyond several Debye lengths away from the fixed charges on the mesostructure of interest (i.e. concentration <0.15 M). This model has been applied to calculate properties of molecules in solution. Extensions of the PB model to account for the finite size of the mobile ions and correlations have been developed. The PB model accounts for solvent molecules implicitly via the dielectric constant is low within the molecule, and is assumed to increase gradually to an unperturbed, far-field value over several angstroms away from the molecule of interest. Solutions to simple problems with spherical, cylindrical, or planar symmetry are available for the linearized PB equation. A closed form formula for the solution of nonlinear PB equations does not seem to be possible except for the planar case and the infinite cylindrical symmetric case, where only counter-ions exist in the solution.

In attempting to use numerical methods to solve the PB equation with complicated molecular geometries and charge distributions, a number of challenges present themselves. The equation is highly nonlinear due to the exponential dependence of the mobile ion concentrations on the potential. The electrical potential varies on more than one length scale (the size of an atom, the Debye length, and the size of the structure of interest). In addition, the length required for the dielectric constant to reach its far field value has a considerable effect on the potential.

One of the most common approaches used to solve the linear and the nonlinear PB is the finite difference formulation where spatial derivatives are approximated using neighboring points. A successive over-relaxation method has been used to get rapid convergence in solving linear systems obtained from the finite difference discretization. Recently, the finite element method has been successfully used to solve the PB equation for larger systems. An adaptive multilevel approach based on tetrahedral elements to create a dense mesh to capture the dielectric discontinuity across the molecular surface has been used. Although irregular grids can be used in a finite element approach, it is inherently less accurate than regular grid methods due to neighboring large and small elements (unless a gradual increase in the element size is imposed). Furthermore, one might argue that in a physically relevant model, the dielectric discontinuity is not present; rather, gradually increases to its unperturbed bulk value with distance from the molecule­medium boundary. Also, regular grids lead to simpler and higher order numerical schemes and they allow implementation of multi-grid methods. The modified PB equation agrees with the results obtained from Monte Carlo simulations on 1:1 and 2:1 ionic solutions.

The development of our Poisson-Boltzmann solver was motivated by the need to have fine meshes to resolve the structure of, and the potential distribution around, large macromolecules and their complexes. A novel numerical approach is developed to reduce the memory requirement by approximately one order of magnitude in comparison to a Galerkin finite element approach. Our solution of the nonlinear PB equation is obtained from the time evolution of a diffusion­advection equation where there is no need to use Newton-Raphson-like methods. A detailed verification of the solution obtained from our approach was accomplished by comparison with the solution of a charged spherical particle. Results on the polio virus protomer and pentamer illustrated the viability of our approach for large systems. On this website you may solve the PB problem for a macro-molecule immersed in an electrolyte.