Michael Holst

Electrostatics of nanosystems: Application to microtubules and the ribosome

Evaluation of the electrostatic properties of biomolecules has become a standard practice in molecular biophysics. Foremost among the models used to elucidate the electrostatic potential is the Poisson-Boltzmann equation; however, existing methods for solving this equation have limited the scope of accurate electrostatic calculations to relatively small biomolecular systems. Here we present the application of numerical methods to enable the trivially parallel solution of the Poisson-Boltzmann equation for supramolecular structures that are orders of magnitude larger in size. As a demonstration of this methodology, electrostatic potentials have been calculated for large microtubule and ribosome structures. The results point to the likely role of electrostatics in a variety of activities of these structures.

Multigrid solution of the Poisson-Boltzmann equation

A multigrid method is presented for the numerical solution of the linearized Poisson–Boltzmann equation arising in molecular biophysics. The equation is discretized with the finite volume method, and the numerical solution of the discrete equations is accomplished with multiple grid techniques originally developed for twodimensional interface problems occurring in reactor physics. A detailed analysis of the resulting method is presented for several computer architectures, including comparisons to diagonally scaled CG, ICCG, vectorized ICCG and MICCG, and to SOR provided with an optimal relaxation parameter. Our results indicate that the multigrid method is superior to the preconditioned CG methods and SOR and that the advantage of multigrid grows with the problem size.

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.

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.

Three-dimensional electron microscopy reveals new details of membrane systems for Ca2+ signaling in the heart.

In the current study, the three-dimensional (3D) topologies of dyadic clefts and associated membrane organelles were mapped in mouse ventricular myocardium using electron tomography. The morphological details and the distribution of membrane systems, including transverse tubules (T-tubules), junctional sarcoplasmic reticulum (SR) and vicinal mitochondria, were determined and presumed to be crucial for controlling cardiac Ca2+ dynamics. The geometric complexity of T-tubules that varied in diameter with frequent branching was clarified. Dyadic clefts were intricately shaped and remarkably small (average 4.39×105 nm3, median 2.81×105 nm3). Although a dyadic cleft of average size could hold maximum 43 ryanodine receptor (RyR) tetramers, more than one-third of clefts were smaller than the size that is able to package as many as 15 RyR tetramers. The dyadic clefts were also adjacent to one another (average end-to-end distance to the nearest dyadic cleft, 19.9 nm) and were distributed irregularly along T-tubule branches. Electron-dense structures that linked membrane organelles were frequently observed between mitochondrial outer membranes and SR or T-tubules. We, thus, propose that the topology of dyadic clefts and the neighboring cellular micro-architecture are the major determinants of the local control of Ca2+ in the heart, including the establishment of the quantal nature of SR Ca2+ releases (e.g. Ca2+ sparks).

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.

Adaptive Numerical Treatment of Elliptic Systems on Manifolds

Adaptive multilevel finite element methods are developed and analyzed for certain elliptic systems arising in geometric analysis and general relativity. This class of nonlinear elliptic systems of tensor equations on manifolds is first reviewed, and then adaptive multilevel finite element methods for approximating solutions to this class of problems are considered in some detail. Two a posteriori error indicators are derived, based on local residuals and on global linearized adjoint or dual problems. The design of Manifold Code (MC) is then discussed; MC is an adaptive multilevel finite element software package for 2- and 3-manifolds developed over several years at Caltech and UC San Diego. It employs a posteriori error estimation, adaptive simplex subdivision, unstructured algebraic multilevel methods, global inexact Newton methods, and numerical continuation methods for the numerical solution of nonlinear covariant elliptic systems on 2- and 3-manifolds. Some of the more interesting features of MC are described in detail, including some new ideas for topology and geometry representation in simplex meshes, and an unusual partition of unity-based method for exploiting parallel computers. A short example is then given which involves the Hamiltonian and momentum constraints in the Einstein equations, a representative nonlinear 4-component covariant elliptic system on a Riemannian 3-manifold which arises in general relativity. A number of operator properties and solvability results recently established in [55] are first summarized, making possible two quasi-optimal a priori error estimates for Galerkin approximations which are then derived. These two results complete the theoretical framework for effective use of adaptive multilevel finite element methods. A sample calculation using the MC software is then presented; more detailed examples using MC for this application may be found in [26].

A New Paradigm for Parallel Adaptive Meshing Algorithms

We present a new approach to the use of parallel computers with adaptive finite element methods. This approach addresses the load balancing problem in a new way, requiring far less communication than current approaches. It also allows existing sequential adaptive PDE codes such as PLTMG and MC to run in a parallel environment without a large investment in recoding. In this new approach, the load balancing problem is reduced to the numerical solution of a small elliptic problem on a single processor, using a sequential adaptive solver, without requiring any modifications to the sequential solver. The small elliptic problem is used to produce a posteriori error estimates to predict future element densities in the mesh, which are then used in a weighted recursive spectral bisection of the initial mesh. The bulk of the calculation then takes place independently on each processor, with no communication, using possibly the same sequential adaptive solver. Each processor adapts its region of the mesh independently, and a nearly load-balanced mesh distribution is usually obtained as a result of the initial weighted spectral bisection. Only the initial fan-out of the mesh decomposition to the processors requires communication. Two additional steps requiring boundary exchange communication may be employed after the individual processors reach an adapted solution, namely, the construction of a global conforming mesh from the independent subproblems, followed by a final smoothing phase using the subdomain solutions as an initial guess. We present a series of convincing numerical experiments which illustrate the effectiveness of this approach. The justification of the initial refinement prediction step, as well as the justification of skipping the two communication-intensive steps, is supported by some recent [J. Xu and A. Zhou, Math. Comp., to appear] and not so recent [J. A. Nitsche and A. H. Schatz, Math. Comp., 28 (1974), pp. 937--958; A. H. Schatz and L. B. Wahlbin, Math. Comp., 31 (1977), pp. 414--442; A. H. Schatz and L. B. Wahlbin, Math. Comp., 64 (1995), pp. 907--928] results on local a priori and a posteriori error estimation.

The adaptive multilevel finite element solution of the Poisson-Boltzmann equation on massively parallel computers

By using new methods for the parallel solution of elliptic partial differential equations, the teraflops computing power of massively parallel computers can be leveraged to perform electrostatic calculations on large biological systems. This paper describes the adaptive multilevel finite element solution of the Poisson-Boltzmann equation for a microtubule on the NPACI Blue Horizon--a massively parallel IBM RS/6000® SP with eight POWER3 SMP nodes. The microtubule system is 40 nm in length and 24 nm in diameter, consists of roughly 600000 atoms, and has a net charge of -1800 e. Poisson-Boltzmann calculations are performed for several processor configurations, and the algorithm used shows excellent parallel scaling.

The Finite Element Approximation of the Nonlinear Poisson-Boltzmann Equation

A widely used electrostatics model in the biomolecular modeling community, the nonlinear Poisson-Boltzmann equation, along with its finite element approximation, are analyzed in this paper. A regularized Poisson-Boltzmann equation is introduced as an auxiliary problem, making it possible to study the original nonlinear equation with delta distribution sources. A priori error estimates for the finite element approximation are obtained for the regularized Poisson-Boltzmann equation based on certain quasi-uniform grids in two and three dimensions. Adaptive finite element approximation through local refinement driven by an a posteriori error estimate is shown to converge. The Poisson-Boltzmann equation does not appear to have been previously studied in detail theoretically, and it is hoped that this paper will help provide molecular modelers with a better foundation for their analytical and computational work with the Poisson-Boltzmann equation. Note that this article apparently gives the first rigorous convergence result for a numerical discretization technique for the nonlinear Poisson-Boltzmann equation with delta distribution sources, and it also introduces the first provably convergent adaptive method for the equation. This last result is currently one of only a handful of existing convergence results of this type for nonlinear problems.