Grid positioning independence and the reduction of self-energy in the solution of the Poisson—Boltzmann equation

A common problem in the solution of the Poisson–Boltzmann equation using finite difference methods is the self-energy of the system, also known as the grid energy. Because atoms are typically modeled as a point charge, the infinite self-energy of a point charge is likewise modeled. In this article, a simple, alternate treatment of atomic charge is described where each atom is represented as a sphere of uniform charge. Unlike the point charge model, this method converges as the grid spacing is reduced. The uniform charge model generates the same electrostatic field outside the atoms. In addition, the use of fine grids reduces the variations in the potential due to variations in the position of atoms relative to the grid. Calculations of Born ion solvation energies, small-molecule solvation energies, and the electrostatic field of superoxide dismutase are used to demonstrate that this method yields the same results as the point charge model. © John Wiley & Sons, Inc.     

Finite-difference solution of the Poisson–Boltzmann equation: Complete elimination of self-energy

A new procedure to solve the Poisson–Boltzmann equation is proposed and shown to be efficient. The electrostatic potential due to the reaction field is calculated directly. Self-interactions among the charges are completely eliminated. Therefore, the reference calculation to cancel out the self-energy is not needed. © 1996 by John Wiley & Sons, Inc.     

Multi-multigrid solution of modified Poisson–Boltzmann equation for arbitrarily shaped molecules

A new multi-multigrid method is presented for solving the modified Poisson–Boltzmann equation based on the Kirkwood Hierarchy of equations, with Loeb's closure, on a three-dimensional grid. The results are compared with standard Poisson–Boltzmann calculations, which are known to underestimate the local concentration of counterions near charged parts of molecules, mainly due to neglect of fluctuations in the ionic concentrations. In the present study, the Kirkwood hierarchy of equations is discretized with the finite volume method and solved using multigrid techniques. The new possibility for solution of the three-dimensional modified Poisson–Boltzmann equation, for the first time within a model including a dielectric discontinuity, and within reasonable computational time, enables the calculation of higher valence ion distributions around arbitrarily shaped biological macromolecules. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 893–901, 1998     

Solving the finite-difference non-linear Poisson–Boltzmann equation

The Poisson–Boltzmann equation can be used to calculate the electrostatic potential field of a molecule surrounded by a solvent containing mobile ions. The Poisson–Boltzmann equation is a non-linear partial differential equation. Finite-difference methods of solving this equation have been restricted to the linearized form of the equation or a finite number of non-linear terms. Here we introduce a method based on a variational formulation of the electrostatic potential and standard multi-dimensional maximization methods that can be used to solve the full non-linear equation. © 1992 by John Wiley & Sons, Inc.     

Numerical solution of the nonlinear Poisson–Boltzmann equation: Developing more robust and efficient methods

We present a robust and efficient numerical method for solution of the nonlinear Poisson-Boltzmann equation arising in molecular biophysics. The equation is discretized with the box method, and solution of the discrete equations is accomplished with a global inexact-Newton method, combined with linear multilevel techniques we have described in an article appearing previously in this journal. A detailed analysis of the resulting method is presented, with comparisons to other methods that have been proposed in the literature, including the classical nonlinear multigrid method, the nonlinear conjugate gradient method, and nonlinear relaxation methods such as successive overrelaxation. Both theoretical and numerical evidence suggests that this method will converge in the case of molecules for which many of the existing methods will not. In addition, for problems which the other methods are able to solve, numerical experiments show that the new method is substantially more efficient, and the superiority of this method grows with the problem size. The method is easy to implement once a linear multilevel solver is available and can also easily be used in conjunction with linear methods other than multigrid. © 1995 by John Wiley & Sons, Inc.     

Adaptive multilevel finite element solution of the Poisson–Boltzmann equation I. Algorithms and examples

This article is the first of two articles on the adaptive multilevel finite element treatment of the nonlinear Poisson–Boltzmann equation (PBE), a nonlinear eliptic equation arising in biomolecular modeling. Fast and accurate numerical solution of the PBE is usually difficult to accomplish, due to the presence of discontinuous coefficients, delta functions, three spatial dimensions, unbounded domain, and rapid (exponential) nonlinearity. In this first article, we explain how adaptive multilevel finite element methods can be used to obtain extremely accurate solutions to the PBE with very modest computational resources, and we present some illustrative examples using two well‐known test problems. The PBE is first discretized with piece‐wise linear finite elements over a very coarse simplex triangulation of the domain. The resulting nonlinear algebraic equations are solved with global inexact Newton methods, which we have described in an article appearing previously in this journal. A posteriori error estimates are then computed from this discrete solution, which then drives a simplex subdivision algorithm for performing adaptive mesh refinement. The discretize–solve–estimate–refine procedure is then repeated, until a nearly uniform solution quality is obtained. The sequence of unstructured meshes is used to apply multilevel methods in conjunction with global inexact Newton methods, so that the cost of solving the nonlinear algebraic equations at each step approaches optimal O(N) linear complexity. All of the numerical procedures are implemented in MANIFOLD CODE (MC), a computer program designed and built by the first author over several years at Caltech and UC San Diego. MC is designed to solve a very general class of nonlinear elliptic equations on complicated domains in two and three dimensions. We describe some of the key features of MC, and give a detailed analysis of its performance for two model PBE problems, with comparisons to the alternative methods. It is shown that the best available uniform mesh‐based finite difference or box‐method algorithms, including multilevel methods, require substantially more time to reach a target PBE solution accuracy than the adaptive multilevel methods in MC. In the second article, we develop an error estimator based on geometric solvent accessibility, and present a series of detailed numerical experiments for several complex biomolecules. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 1319–1342, 2000     

pKA in proteins solving the Poisson–Boltzmann equation with finite elements

Knowledge on pK_A values is an eminent factor to understand the function of proteins in living systems. We present a novel approach demonstrating that the finite element (FE) method of solving the linearized Poisson–Boltzmann equation (lPBE) can successfully be used to compute pK_A values in proteins with high accuracy as a possible replacement to finite difference (FD) method. For this purpose, we implemented the software molecular Finite Element Solver (mFES) in the framework of the Karlsberg+ program to compute pK_A values. This work focuses on a comparison between pK_A computations obtained with the well‐established FD method and with the new developed FE method mFES, solving the lPBE using protein crystal structures without conformational changes. Accurate and coarse model systems are set up with mFES using a similar number of unknowns compared with the FD method. Our FE method delivers results for computations of pK_A values and interaction energies of titratable groups, which are comparable in accuracy. We introduce different thermodynamic cycles to evaluate pK_A values and we show for the FE method how different parameters influence the accuracy of computed pK_A values. © 2015 Wiley Periodicals, Inc.     

An automatic three-dimensional finite element mesh generation system for the Poisson–Boltzmann equation

We present an automatic three-dimensional mesh generation system for the solution of the Poisson–Boltzmann equation using a finite element discretization. The different algorithms presented allow the construction of a tetrahedral mesh using a predetermined spatial distribution of vertices adapted to the geometry of the dielectric continuum solvent model. A constrained mesh generation strategy, based on Bowyer's algorithm, is used to construct the tetrahedral elements incrementally and embed the Richards surface of the molecule into the mesh as a set of triangular faces. A direct mesh construction algorithm is then used to refine the existing mesh in the neighborhood of the dielectric interface. This will allow an accurate calculation of the induced polarization charge to be carried out while maintaining a sparse grid structure in the rest of the computational space. The inclusion of an ionic boundary at some finite distance from the dielectric interface can be automatically achieved as the grid point distribution outside the solute molecule is constructed using a set of surfaces topologically equivalent to this boundary. The meshes obtained by applying the algorithm to real molecular geometries are described. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1570–1590, 1997     

Solving the finite-difference,nonlinear, Poisson–Boltzmann equation under a linear approach

Electrostatic interactions are among the key factors in determining the structure and function of biomolecules. Simulating such interactions involves solving the Poisson equation and the Poisson-Boltzmann (P-B) equation in the molecular interior and exterior region, respectively. The P-B equation is a nonlinear partial differential equation. The central processing unit (CPU) time for solving the full nonlinear P-B equation has been severalfold greater than the equivalent linear case. Here a simple method is proposed to solve the full nonlinear P-B equation under a linear approach, which has been tested both on a spherical case and on small molecules. Results show that our method converges rapidly even under highly charged cases. With this method, the CPU time for solving the full nonlinear P-B equation is somewhat less than the equivalent linear case in our calculations. © 1995 by John Wiley & Sons, Inc.     

Ion strength limit of computed excess functions based on the linearized Poisson–Boltzmann equation

The linearized Poisson–Boltzmann (L‐PB) equation is examined for its κ‐range of validity (κ, Debye reciprocal length). This is done for the Debye–Hückel (DH) theory, i.e., using a single ion size, and for the SiS treatment (D. Fraenkel, Mol. Phys. 2010 , 108, 1435), which extends the DH theory to the case of ion‐size dissimilarity (therefore dubbed DH–SiS). The linearization of the PB equation has been claimed responsible for the DH theory's failure to fit with experiment at > 0.1 m; but DH–SiS fits with data of the mean ionic activity coefficient, γ_± (molal), against m, even at m > 1 (κ > 0.33 Å^?1). The SiS expressions combine the overall extra‐electrostatic potential energy of the smaller ion, as central ion—Ψ_a>b(κ), with that of the larger ion, as central ion—Ψ_b>a(κ); a and b are, respectively, the counterion and co‐ion distances of closest approach. Ψ_a>b and Ψ_b>a are derived from the L‐PB equation, which appears to conflict with their being effective up to moderate electrolyte concentrations (≈1 m). However, the L‐PB equation can be valid up to κ ≥ 1.3 Å^?1 if one abandons the 1/κ criterion for its effectiveness and, instead, use, as criterion, the mean‐field electrostatic interaction potential of the central ion with its ion cloud, at a radial distance dividing the cloud charge into two equal parts. The DH theory's failure is, thus, not because of using the L‐PB equation; the lethal approximation is assigning a single size to the positive and negative ions. © 2015 Wiley Periodicals, Inc.     

Adaptive multilevel finite element solution of the Poisson–Boltzmann equation II. Refinement at solvent‐accessible surfaces in biomolecular systems

We apply the adaptive multilevel finite element techniques (Holst, Baker, and Wang 21 ) to the nonlinear Poisson–Boltzmann equation (PBE) in the context of biomolecules. Fast and accurate numerical solution of the PBE in this setting is usually difficult to accomplish due to presence of discontinuous coefficients, delta functions, three spatial dimensions, unbounded domains, and rapid (exponential) nonlinearity. However, these adaptive techniques have shown substantial improvement in solution time over conventional uniform‐mesh finite difference methods. One important aspect of the adaptive multilevel finite element method is the robust a posteriori error estimators necessary to drive the adaptive refinement routines. This article discusses the choice of solvent accessibility for a posteriori error estimation of PBE solutions and the implementation of such routines in the “Adaptive Poisson–Boltzmann Solver” (APBS) software package based on the “Manifold Code” (MC) libraries. Results are shown for the application of this method to several biomolecular systems. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 1343–1352, 2000     

A rapid finite difference algorithm,utilizing successive over-relaxation to solve the Poisson–Boltzmann equation

An efficient algorithm is presented for the numerical solution of the Poisson–Boltzmann equation by the finite difference method of successive over-relaxation. Improvements include the rapid estimation of the optimum relaxation parameter, reduction in number of operations per iteration, and vector-oriented array mapping. The algorithm has been incorporated into the electrostatic program DelPhi, reducing the required computing time by between one and two orders of magnitude. As a result the estimation of electrostatic effects such as solvent screening, ion distributions, and solvation energies of small solutes and biological macromolecules in solution, can be performed rapidly, and with minimal computing facilities.     

A first‐order system least‐squares finite element method for the Poisson‐Boltzmann equation

The Poisson‐Boltzmann equation is an important tool in modeling solvent in biomolecular systems. In this article, we focus on numerical approximations to the electrostatic potential expressed in the regularized linear Poisson‐Boltzmann equation. We expose the flux directly through a first‐order system form of the equation. Using this formulation, we propose a system that yields a tractable least‐squares finite element formulation and establish theory to support this approach. The least‐squares finite element approximation naturally provides an a posteriori error estimator and we present numerical evidence in support of the method. The computational results highlight optimality in the case of adaptive mesh refinement for a variety of molecular configurations. In particular, we show promising performance for the Born ion, Fasciculin 1, methanol, and a dipole, which highlights robustness of our approach. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Overcharging in Colloids: Beyond the Poisson–Boltzmann Approach

A broad range of manufactured products and biological fluids are colloids. The ability to understand and control the processes (of scientific, technological and industrial interest) in which such colloids are involved relies upon a precise knowledge of the electrical double layer. The traditional approach to describing this ion cloud around colloidal particles has been the Gouy–Chapman model, developed on the basis of the Poisson–Boltzmann equation. Since the early 1980s, however, more sophisticated theoretical treatments have revealed both quantitative and qualitative deficiencies in the Poisson–Boltzmann theory, particularly at high ionic strengths and/or high surface charge densities. This review deals with these novel approaches, which are mostly computer simulations and approximate integral equation theories based on the so‐called primitive model. Special attention is paid to phenomena that cannot be accounted for by the classic theory as a result of neglecting ion size correlations, such as overcharging, namely, the counterion concentration in the immediate neighborhood of the surface is so large that the particle surface is overcompensated. Other illustrative examples are the nonmonotonic behavior of the electrostatic potential and attractive interactions between equally charged surfaces. These predictions are certainly remarkable and, on paper, they can have an effect on experimentally measurable quantities (for instance, electrophoretic mobility). Even so, these new approaches have scarcely been applied in practice. Thus a critical survey on the relevance of ion size correlation in real systems is also included. Overcharging of macroions can also be brought about by adsorption of oppositely charged polyelectrolytes. Noteworthy examples and theoretical approaches for them are also briefly reviewed.     

An assessment of the accuracy of the RRIGS hydration potential: Comparison to solutions of the Poisson–Boltzmann equation

A rapid, pairwise hydration potential, the reduced radius independent Gaussian sphere (RRIGS) approximation, has been presented recently. Because experimental values of the conformational dependence of the hydration free energy are unavailable, this hydration potential is testable by comparison to a presumably more accurate (yet more computationally intensive) model. One such method is the electrostatic hydration approach, which models the protein as a collection of point charges in a low-dielectric medium and the solvent as a high-dielectric continuum. The electrostatic free energy can be determined by solving the Poisson–Boltzmann equation, which is carried out with the program DelPhi. The total free energy of hydration is calculated by adding a free energy of cavity formation term to this electrostatic term. Comparison is made for many conformations of two proteins, bovine pancreatic trypsin inhibitor (BPTI) and the carboxy-terminal fragment of the L7/L12 ribosomal protein (CTF). Thirty-nine near-native structures of BPTI, previously generated by Ripoll and coworkers, and 150 conformations of CTF, generated by a threading algorithm to cover a wide range of conformational space, were used in these comparisons. It is shown that, for the neutral forms of these proteins, the RRIGS hydration potential correlates very well with the electrostatic model hydration free energy, although the correlation is better for the CTF conformations than for the near-native BPTI conformations. For charged forms, the correlation is much poorer. These results serve as evidence that solvent-exposure models of hydration, which leave out cooperative effects between different groups, may be appropriate for modeling neutral or slightly charged species, because these cooperative effects are likely to be small. However, for highly charged species where cooperative effects are surely large, such an approach will be less accurate. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 :1072–1078, 1997     

Assessment of linear finite‐difference Poisson–Boltzmann solvers

CPU time and memory usage are two vital issues that any numerical solvers for the Poisson–Boltzmann equation have to face in biomolecular applications. In this study, we systematically analyzed the CPU time and memory usage of five commonly used finite‐difference solvers with a large and diversified set of biomolecular structures. Our comparative analysis shows that modified incomplete Cholesky conjugate gradient and geometric multigrid are the most efficient in the diversified test set. For the two efficient solvers, our test shows that their CPU times increase approximately linearly with the numbers of grids. Their CPU times also increase almost linearly with the negative logarithm of the convergence criterion at very similar rate. Our comparison further shows that geometric multigrid performs better in the large set of tested biomolecules. However, modified incomplete Cholesky conjugate gradient is superior to geometric multigrid in molecular dynamics simulations of tested molecules. We also investigated other significant components in numerical solutions of the Poisson–Boltzmann equation. It turns out that the time‐limiting step is the free boundary condition setup for the linear systems for the selected proteins if the electrostatic focusing is not used. Thus, development of future numerical solvers for the Poisson–Boltzmann equation should balance all aspects of the numerical procedures in realistic biomolecular applications. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

A smooth permittivity function for Poisson–Boltzmann solvation methods

This work introduces a continuous smooth permittivity function into Poisson–Boltzmann techniques for continuum approaches to modeling the solvation of small molecules and proteins. The permittivity function is derived using a Gaussian method to describe volume exclusion. The new method allows a rigorous determination of solvent forces within a grid‐based technology. The generality of approach is demonstrated by considering a range of applications for small molecules and macromolecules. We also present a very complete statistical analysis of grid errors, and show that the accuracy of our Gaussian‐based method is improved over standard techniques. The method has been implemented in a new code called ZAP, which is freely available to academic institutions. 1 © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 608–640, 2001     

SMPBS: Web server for computing biomolecular electrostatics using finite element solvers of size modified Poisson‐Boltzmann equation

SMPBS (Size Modified Poisson‐Boltzmann Solvers) is a web server for computing biomolecular electrostatics using finite element solvers of the size modified Poisson‐Boltzmann equation (SMPBE). SMPBE not only reflects ionic size effects but also includes the classic Poisson‐Boltzmann equation (PBE) as a special case. Thus, its web server is expected to have a broader range of applications than a PBE web server. SMPBS is designed with a dynamic, mobile‐friendly user interface, and features easily accessible help text, asynchronous data submission, and an interactive, hardware‐accelerated molecular visualization viewer based on the 3Dmol.js library. In particular, the viewer allows computed electrostatics to be directly mapped onto an irregular triangular mesh of a molecular surface. Due to this functionality and the fast SMPBE finite element solvers, the web server is very efficient in the calculation and visualization of electrostatics. In addition, SMPBE is reconstructed using a new objective electrostatic free energy, clearly showing that the electrostatics and ionic concentrations predicted by SMPBE are optimal in the sense of minimizing the objective electrostatic free energy. SMPBS is available at the URL: smpbs.math.uwm.edu © 2017 Wiley Periodicals, Inc.     

A numerical solution of the linear Boltzmann equation using cubic B-splines

A numerical method using cubic B-splines is presented for solving the linear Boltzmann equation. The collision kernel for the system is chosen as the Wigner-Wilkins kernel. A total of three different representations for the distribution function are presented. Eigenvalues and eigenfunctions of the collision matrix are obtained for various mass ratios and compared with known values. Distribution functions, along with first and second moments, are evaluated for different mass and temperature ratios. Overall it is shown that the method is accurate and well behaved. In particular, moments can be predicted with very few points if the representation is chosen well. This method produces sparse matrices, can be easily generalized to higher dimensions, and can be cast into efficient parallel algorithms.