Hybrid boundary element and finite difference method for solving the nonlinear Poisson-Boltzmann equation

Electrodiffusion: a continuum modeling framework for biomolecular systems with realistic spatiotemporal resolution

A computational framework is presented for the continuum modeling of cellular biomolecular diffusion influenced by electrostatic driving forces. This framework is developed from a combination of state-of-the-art numerical methods, geometric meshing, and computer visualization tools. In particular, a hybrid of (adaptive) finite element and boundary element methods is adopted to solve the Smoluchowski equation (SE), the Poisson equation (PE), and the Poisson-Nernst-Planck equation (PNPE) in order to describe electrodiffusion processes. The finite element method is used because of its flexibility in modeling irregular geometries and complex boundary conditions. The boundary element method is used due to the convenience of treating the singularities in the source charge distribution and its accurate solution to electrostatic problems on molecular boundaries. Nonsteady-state diffusion can be studied using this framework, with the electric field computed using the densities of charged small molecules and mobile ions in the solvent. A solution for mesh generation for biomolecular systems is supplied, which is an essential component for the finite element and boundary element computations. The uncoupled Smoluchowski equation and Poisson-Boltzmann equation are considered as special cases of the PNPE in the numerical algorithm, and therefore can be solved in this framework as well. Two types of computations are reported in the results: stationary PNPE and time-dependent SE or Nernst-Planck equations solutions. A biological application of the first type is the ionic density distribution around a fragment of DNA determined by the equilibrium PNPE. The stationary PNPE with nonzero flux is also studied for a simple model system, and leads to an observation that the interference on electrostatic field of the substrate charges strongly affects the reaction rate coefficient. The second is a time-dependent diffusion process: the consumption of the neurotransmitter acetylcholine by acetylcholinesterase, determined by the SE and a single uncoupled solution of the Poisson-Boltzmann equation. The electrostatic effects, counterion compensation, spatiotemporal distribution, and diffusion-controlled reaction kinetics are analyzed and different methods are compared.

Treatment of charge singularities in implicit solvent models

This paper presents a novel method for solving the Poisson-Boltzmann (PB) equation based on a rigorous treatment of geometric singularities of the dielectric interface and a Green's function formulation of charge singularities. Geometric singularities, such as cusps and self-intersecting surfaces, in the dielectric interfaces are bottleneck in developing highly accurate PB solvers. Based on an advanced mathematical technique, the matched interface and boundary (MIB) method, we have recently developed a PB solver by rigorously enforcing the flux continuity conditions at the solvent-molecule interface where geometric singularities may occur. The resulting PB solver, denoted as MIBPB-II, is able to deliver second order accuracy for the molecular surfaces of proteins. However, when the mesh size approaches half of the van der Waals radius, the MIBPB-II cannot maintain its accuracy because the grid points that carry the interface information overlap with those that carry distributed singular charges. In the present Green's function formalism, the charge singularities are transformed into interface flux jump conditions, which are treated on an equal footing as the geometric singularities in our MIB framework. The resulting method, denoted as MIBPB-III, is able to provide highly accurate electrostatic potentials at a mesh as coarse as 1.2 A for proteins. Consequently, at a given level of accuracy, the MIBPB-III is about three times faster than the APBS, a recent multigrid PB solver. The MIBPB-III has been extensively validated by using analytically solvable problems, molecular surfaces of polyatomic systems, and 24 proteins. It provides reliable benchmark numerical solutions for the PB equation.

Treatment of geometric singularities in implicit solvent models

Geometric singularities, such as cusps and self-intersecting surfaces, are major obstacles to the accuracy, convergence, and stability of the numerical solution of the Poisson-Boltzmann (PB) equation. In earlier work, an interface technique based PB solver was developed using the matched interface and boundary (MIB) method, which explicitly enforces the flux jump condition at the solvent-solute interfaces and leads to highly accurate biomolecular electrostatics in continuum electric environments. However, such a PB solver, denoted as MIBPB-I, cannot maintain the designed second order convergence whenever there are geometric singularities, such as cusps and self-intersecting surfaces. Moreover, the matrix of the MIBPB-I is not optimally symmetrical, resulting in the convergence difficulty. The present work presents a new interface method based PB solver, denoted as MIBPB-II, to address the aforementioned problems. The present MIBPB-II solver is systematical and robust in treating geometric singularities and delivers second order convergence for arbitrarily complex molecular surfaces of proteins. A new procedure is introduced to make the MIBPB-II matrix optimally symmetrical and diagonally dominant. The MIBPB-II solver is extensively validated by the molecular surfaces of few-atom systems and a set of 24 proteins. Converged electrostatic potentials and solvation free energies are obtained at a coarse grid spacing of 0.5 A and are considerably more accurate than those obtained by the PBEQ and the APBS at finer grid spacings.

Monte Carlo and modified Tanford-Kirkwood results for macromolecular electrostatics calculations

The understanding of electrostatic interactions is an essential aspect of the complex correlation between structure and function of biological macromolecules. It is also important in protein engineering and design. Theoretical studies of such interactions are predominantly done within the framework of Debye-Hückel theory. A classical example is the Tanford-Kirkwood (TK) model. Besides other limitations, this model assumes an infinitesimally small macromolecule concentration. By comparison to Monte Carlo (MC) simulations, it is shown that TK predictions for the shifts in ion binding constants upon addition of salt become less reliable even at moderately macromolecular concentrations. A simple modification based on colloidal literature is suggested to the TK scheme. The modified TK models suggested here satisfactorily predict MC and experimental shifts in the calcium binding constant as a function of protein concentration for the calbindin D(9k) mutant and calmodulin.

A new outer boundary formulation and energy corrections for the nonlinear Poisson-Boltzmann equation

The nonlinear Poisson-Boltzmann equation (PBE) has been successfully used for the prediction of numerous electrostatic properties of highly charged biopolyelectrolytes immersed in aqueous salt solutions. While numerous numerical solvers for the 3D PBE have been developed, the formulation of the outer boundary treatments used in these methods has only been loosely addressed, especially in the nonlinear case. The de facto standard in current nonlinear PBE implementations is to either set the potential at the outer boundaries to zero or estimate it using the (linear) Debye-Hückel (DH) approximation. However, an assessment of how these outer boundary treatments affect the overall solution accuracy does not appear to have been previously made. As will be demonstrated here, both approximations can, under certain conditions, produce completely erroneous estimates of the potential and energy salt dependencies. A related concern for calculations carried out on grids of finite extent (e.g., all current finite difference and finite element implementations) is the contribution to the energy and salt dependence from the exterior region outside the computational grid. This too is shown to be significant, especially at low salt concentration where essentially all of the contributions to the excess osmotic pressure and ion stress energies originate from this exterior region. In this paper the authors introduce a new outer boundary treatment that is valid for both the linear and nonlinear PBE. The authors also formulate energy corrections to account for contributions from outside the computational domain. Finally, the authors also consider the effects of general ion exclusion layers upon biomolecular electrostatics. It is shown that while these layers tend to increase the surface electrostatic potential, under physiological salt conditions and high net charges their effect on the excess osmotic pressure term, which is a measure of the salt dependence of the total electrostatic free energy, is weak. To facilitate presentation and allow very fine resolutions and/or large computational domains to be considered, attention is restricted to the 1D spherically symmetric nonlinear PBE. Though geometrically limited, the modeling principles nevertheless extend to general PBE solvers as discussed in the Appendix. The 1D model can also be used to benchmark and validate the salt effect prediction capabilities of existing PBE solvers.

Recognition of the spliceosomal branch site RNA helix on the basis of surface and electrostatic features

We have investigated electrostatic and surface features of an essential region of the catalytic core of the spliceosome, the eukaryotic precursor messenger (pre-m)RNA splicing apparatus. The nucleophile for the first of two splicing reactions is the 2′-hydroxyl (OH) of the ribose of a specific adenosine within the intron. During assembly of the spliceosome's catalytic core, this adenosine is positioned by pairing with a short region of the U2 small nuclear (sn)RNA to form the pre-mRNA branch site helix. The solution structure of the spliceosomal pre-mRNA branch site [Newby,M.I. and Greenbaum,N.L. (2002) Nature Struct. Biol., 9, 958–965] showed that a phylogenetically conserved pseudouridine (ψ) residue in the segment of U2 snRNA that pairs with the intron induces a markedly different structure compared with that of its unmodified counterpart. In order to achieve a more detailed understanding of the factors that contribute to recognition of the spliceosome's branch site helix and activation of the nucleophile for the first step of pre-mRNA splicing, we have calculated surface areas and electrostatic potentials of ψ-modified and unmodified branch site duplexes. There was no significant difference between the total accessible area or ratio of total polar:nonpolar groups between modified and unmodified duplexes. However, there was substantially greater exposure of nonpolar area of the adenine base, and less exposure of the 2′-OH, in the ψ-modified structure. Electrostatic potentials computed using a hybrid boundary element and finite difference nonlinear Poisson–Boltzmann approach [Boschitsch, A.H. and Fenley, M.O. (2004) J. Comput. Chem., 25, 935–955] revealed a region of exceptionally negative potential in the major groove surrounding the 2′-OH of the branch site adenosine. These surface and electrostatic features may contribute to the overall recognition of the pre-mRNA branch site region by other components of the splicing reaction.

Structural and Thermodynamic Evidence for a Stabilizing Role of Nop5p in S-Adenosyl-L-methionine Binding to Fibrillarin

In Archaea, fibrillarin and Nop5p form the core complex of box C/D small ribonucleoprotein particles, which are responsible for site-specific 2'-hydroxyl methylation of ribosomal and transfer RNAs. Fibrillarin has a conserved methyltransferase fold and employs S-adenosyl-l-methionine (AdoMet) as the cofactor in methyl transfer reactions. Comparison between recently determined crystal structures of free fibrillarin and fibrillarin-Nop5p-AdoMet tertiary complex revealed large conformational differences at the cofactor-binding site in fibrillarin. To identify the structural elements responsible for these large conformational differences, we refined a crystal structure of Archaeoglobus fulgidus fibrillarin-Nop5p binary complex at 3.5 A. This structure exhibited a pre-formed backbone geometry at the cofactor binding site similar to that when the cofactor is bound, suggesting that binding of Nop5p alone to fibrillarin is sufficient to stabilize the AdoMet-binding pocket. Calorimetry studies of cofactor binding to fibrillarin alone and to fibrillarin-Nop5p binary complex provided further support for this role of Nop5p. Mutagenesis and thermodynamic data showed that a cation-pi bridge formed between Tyr-89 of fibrillarin and Arg-169 of Nop5p, although dispensable for in vitro methylation activity, could partially account for the enhanced binding of cofactor to fibrillarin by Nop5p. Finally, assessment of cofactor-binding thermodynamics and catalytic activities of enzyme mutants identified three additional fibrillarin residues (Thr-70, Glu-88, and Asp-133) to be important for cofactor binding and for catalysis.

Molecular Basis of Box C/D RNA-Protein Interactions

We have determined and refined a crystal structure of the initial assembly complex of archaeal box C/D sRNPs comprising the Archaeoglobus fulgidus (AF) L7Ae protein and a box C/D RNA. The box C/D RNA forms a classical kink-turn (K-turn) structure and the resulting protein-RNA complex serves as a distinct platform for recruitment of the fibrillarin-Nop5p complex. The cocrystal structure confirms previously proposed secondary structure of the box C/D RNA that includes a protruded U, a UU mismatch, and two sheared tandem GA base pairs. Detailed structural comparisons of the AF L7Ae-box C/D RNA complex with previously determined crystal structures of L7Ae homologs in complex with functionally distinct K-turn RNAs revealed a set of remarkably conserved principles in protein-RNA interactions. These analyses provide a structural basis for interpreting the functional roles of the box C/D sequences in directing specific assembly of box C/D sRNPs.