Protein Electrostatics: Rapid Multigrid-Based Newton Algorithm for Solution of the Full Nonlinear Poisson-Boltzmann Equation

Abstract

A new method for solving the full nonlinear Poisson-Boltzmann equation is outlined. This method is robust and efficient, and uses a combination of the multigrid and inexact Newton algorithms. The novelty of this approach lies in the appropriate combination of the two methods, neither of which by themselves are capable of solving the nonlinear problem accurately. Features of the Poisson-Boltzmann equation are fully exploited by each component of the hybrid algorithm to provide robustness and speed. The advantages inherent in this method increase with the size of the problem. The efficacy of the method is illustrated by calculations of the electrostatic potential around the enzyme Superoxide Dismutase. The CPU time required to solve the full nonlinear equation is less than half that needed for a conjugate gradient solution of the corresponding linearized Poisson-Boltzmann equation. The solutions reveal that the field around the active sites is significantly reduced as compared to that obtained by solving the corresponding linearized Poisson-Boltzmann equation. This new method for the nonlinear Poisson-Boltzmann equation will enable fast and accurate solutions of large protein electrostatics problems.     

Calculation of the Electrophoretic Mobility of a Particle Bearing Bound Polyelectrolyte Using the Nonlinear Poisson-Boltzmann Equation

A numerical method for determining the electrophoretic mobility of a polyelectrolyte-coated particle is presented. The particle surface is modeled as having a permeable layer of polyelectrolyte molecules anchored to its surface. Fluid flow within the polyelectrolyte layer is subject to Stokes drag arising from the polyelectrolyte segments. The method allows arbitrary distribution of polymer segments and charge density normal to the surface to be used. The hydrodynamic plane of shear may also be varied. The potential profile is determined by a numerical solution to the nonlinearized Poisson-Boltzmann equation. The potential profile is then used in a numerical solution to the Navier-Stokes equation to give the required mobility. The use of the nonlinearized Poisson-Boltzmann equation extends the results to higher charge density/lower ionic strength conditions than previous treatments. The surface potentials and mobilities for three limiting charge distributions are compared for both the linear and nonlinear treatments to delimit the range of validity of the linear treatment. The utility of the numerical, nonlinear treatment is demonstrated by an improved fit to the electrophoretic mobility of human erythrocytes as a function of ionic strength in the range 10 to 150 mM.     

Investigation of ionic stability criteria for ion-permeable charged membranes

The flocculation criteria in the DLVO theory of colloid stability are applied to ion-permeable membranes containing ionizable fixed groups. These groups are not restricted to the membrane surface but are uniformly distributed throughout a thick surface layer. The flocculation concentrations for such membranes are calculated by using a numerical method to solve the nonlinear Poisson-Boltzmann equation. Results are compared with calculations previously carried out for more restrictive models of biological membranes. The flocculation concentrations are shown to depend on the density of ionizable groups, the dissociation constant of these groups, and the pH of the bulk solution.     

Energetics of charge-charge interactions in proteins

Electrostatic interactions between pairs of atoms in proteins are calculated with a model based on the linearized Poisson-Boltzmann equation. The equation is solved accurately by a method that takes into account the detailed shape of the protein. This paper presents applications to several systems. Experimental data for the interaction of ionized residues with an active site histidine in subtilisin BPN' allow the model to be tested, using various assumptions for the electrical properties of the protein and solvent. The electrostatic stabilization of the active site thiolate of rhodanese is analyzed, with attention to the influence of alpha-helices. Finally, relationships between electrostatic potential and charge-charge distance are reported for large and small globular proteins. The above results are compared with those of simpler electrostatic models, including Coulomb's law with both a distance-dependent dielectric constant (epsilon = R) and a fixed dielectric constant (epsilon = 2), and Tanford-Kirkwood theory. The primary conclusions are as follows: 1) The Poisson-Boltzmann model agrees with the subtilisin data over a range of ionic strengths; 2) two alpha-helices generate a large potential in the active site of rhodanese; 3) epsilon = R overestimates weak electrostatic interactions but yields relatively good results for strong ones; 4) Tanford-Kirkwood theory is a useful approximation to detailed solutions of the linearized Poisson-Boltzmann equation in globular proteins; and 5) the modified Tanford-Kirkwood theory over-screens the measured electrostatic interactions in subtilisin.     

Fluid mechanical and electrostatic study on the osmotic flow through cylindrical pores

When two solutions differing in solute concentration are separated by a porous membrane, the osmotic pressure will generate a net volume flux of the suspending fluid across the membrane; this is termed osmotic flow. We consider the osmotic flow across a membrane with circular cylindrical pores when the solute and the pore walls are electrically charged, and the suspending fluid is an electrolytic solution containing small cations and anions. Under the condition in which the radius of the pores and that of the solute molecules greatly exceed those of the solvent as well as the ions, a fluid mechanical and electrostatic theory is introduced to describe the osmotic flow in the presence of electric charge. The interaction energy, including the electrostatic interaction between the solute and the pore wall, plays a key role in determining the osmotic flow. We examine the electrostatic effect on the osmotic flow and discuss the difference in the interaction energy determined from the nonlinear Poisson-Boltzmann equation and from its linearized equation (the Debye-Hückel equation).     

Density functional theory for the nonspecific binding of salt to polyelectrolytes: thermodynamic properties

The thermodynamics of the nonspecific binding of salt to a polyelectrolyte molecule is studied using a density functional approach. The polyelectrolyte molecule is modeled as an infinite, inflexible, and impenetrable charged cylinder and the counterions and co-ions are modeled as charged hard spheres of equal diameter. The density functional theory is based on a hybrid approach where the hard-sphere contribution to the one-particle correlation function is evaluated nonperturbatively and the ionic contribution to the one-particle correlation function is evaluated perturbatively. The advantage of the approach is that analytical expressions are available for all the correlation functions. The calculated single ion preferential interaction coefficients, excess free energy, and activity coefficients show a nonmonotonic variation as a function of polyion charge in the presence of divalent ions. These properties display considerable departure from the predictions of the nonlinear Poisson-Boltzmann (NLPB) equation, with qualitative differences in some cases, which may be attributed to correlation effects neglected in the NLPB theory.     

The theoretical distributions and diffusivities of small ions in chondroitin sulphate and hyaluronate

The electrostatic interactions between polyionic glycosaminoglycans and small mobile ions are investigated using the Poisson-Boltzmann equation and a rod-in-cell model of the polyelectrolyte. Calculations are made for the range of polyelectrolyte concentrations and buffer compositions for which measurements of ion distributions and diffusivities are reported in a companion paper (Maroudas et al., Biophys. Chem. 32 (1988) 257). We conclude that the distribution of mobile ions is largely determined by the 'far-field' potential and is adequately described by the Poisson-Boltzmann theory and also by more approximate theories such as ideal Donnan or 'condensation' theory. The measured variations in cation diffusivities, particularly the increase in diffusivity with increasing matrix concentration at low ionic strengths, are predicted qualitatively using an approximate diffusion theory together with the calculated potential fields. However, the same theory applied to anion diffusion gives qualitatively wrong results.     

Average electrostatic potential between the filaments in striated muscle and its relation to a simple Donnan potential.

An approximate analytical solution to the Poisson-Boltzmann equation for a cylindrical particle was used to calculate the relationship between the charge on the filaments and the average electrostatic potential. Both thick and then filaments were considered in the muscle lattice with a filament charge ratio of 4 to 1. Comparing this with a similar relationship obtained using simple Donnan theory showed a discrepancy at high charge where the Poisson-Boltzmann equation leads to saturation of the average potential. However, using two separate experiments from the literature, we have shown that at pH 7.0 muscle must not be close to saturation and thus is in a region of the curve where the two approaches agree.     

Tension fluctuations in contracting myofibrils and their interpretation

Repulsive pressure has been measured as a function of lattice spacing in gels of tobacco mosaic virus (TMV) and in the filament lattice of vertebrate striated muscle. External pressures up to ten atm have been applied to these lattices by an osmotic stress method. Numerical solutions to the Poisson-Boltzmann equation in hexagonal lattices have been obtained and compared to the TMV and muscle data. The theoretical curves using values for k calculated from the ionic strength give a good fit to experimental data from TMV gels, and an approximate fit to that from the muscle lattice, provided that a charge radius for the muscle thick filaments of approximately 16 nm is assumed. Variations in ionic strength, sarcomere length and state of the muscle give results which agree qualitatively with the theory, though a good fit between experiment and theory in the muscle case will clearly require consideration of other types of forces. We conclude that Poisson-Boltzmann theory can provide a good first approximation to the long-range electrostatic forces operating in such biological gel systems.     

Donnan membrane equilibrium,sedimentation equilibrium,and coil expansion of DNA in salt solutions

This is a review of applications of the McMillan-Mayer-Hill virial theory and the ionic double-layer theory to dilute colloidal solutions, in particular, solutions of DNA. Interactions of highly charged colloidal rods are developed in terms of the second virial coefficients between two rods, and between one rod and one small co-ion. The relevant cluster integrals are evaluated with interaction potentials based on the Poisson-Boltzmann equation. The treatment is extended to the intrachain repulsion responsible for the statistical swelling of coiled DNA (excluded volume effect). The theory is compared with three sets of experimental data: The salt distribution in Donnan membrane equilibria of DNA-salt solutions, sedimentation equilibria of short DNA fragments at different ionic strengths, and the intrinsic viscosity of T7 DNA in NaCl solutions. In all cases the theory agrees well with the experiments. The agreement is not convincing for the sedimentation equilibrium at low ionic strength, because here the experimental DNA concentration is too high for the truncated dilute solution expansion of the DNA-salt repulsion.     

Comparison of Poisson-Boltzmann and condensation model expression for the colligative properties of cylindrical polyions

Closed-form expression have been derived for the polyelectrolyte contribution to the colligative properties of solutions containing rodlike polyions in the presence of excess added salt. The derivations are based on: the conventional Poisson-Boltzmann equation for cylindrical symmetry; the thermodynamics of the cell model developed by Marcus [J. Chem. Phys. 23 , 1057–1068 (1955)]; and an equation derived from the cylindrical Poisson-Boltzmann cell model by Anderson and Record [Biophys. Chem. 11 , 353–360 (1980)]. Subject to the inherent limitations of the Poisson-Boltzmann approximation [Fixman (1979) J. Chem. Phys. 70 , 4995–5005], the resulting expressions are nevertheless applicable outside the “limit of infinite dilution.” They conform over a range of salt concentrations to the limiting laws deduced by Manning from the hypothesis of counterion condensation [J. Chem. Phys. 51 , 924–933 (1969)]. This hypothesis is found to be compatible with the Poisson-Boltzmann cell model but is not required in the derivation of the thermodynamic coefficients presented here. It is demonstrated that the magnitude of the polyion axial charge density plays a critical role in determining the low-salt limiting forms of the colligative properties obtained from the Poisson-Boltzmann equation, in close analogy with Manning's model.     

The electrostatic potential of B-DNA

Electrostatic potentials around DNA are obtained by solving the nonlinear Poisson-Boltzmann (PB) equation. The detailed charge distribution of the DNA and the different polarizabilities of the macromolecule and solvent are included explicitly in the calculations. The PB equation is solved using extensions of a finite difference approach applied previously to proteins. Electrical potentials and ion concentrations are compared to those obtained with simpler models. It is found that the shape of the dielectric boundary between the macromolecule and solvent has significant effects on the calculated potentials near the surface, particularly in the grooves. Sequence-specific patterns are found, the most surprising result being the existence of positive regions of potential near the bases in both the major and minor grooves. The effect of solvent and ionic atmosphere screening of phosphate-phosphate repulsions is studied, and an effective dielectric function, appropriate for molecular mechanics simulations, is derived.     

The Poisson-Boltzmann equation for biomolecular electrostatics: a tool for structural biology

Electrostatics plays a fundamental role in virtually all processes involving biomolecules in solution. The Poisson-Boltzmann equation constitutes one of the most fundamental approaches to treat electrostatic effects in solution. The theoretical basis of the Poisson-Boltzmann equation is reviewed and a wide range of applications is presented, including the computation of the electrostatic potential at the solvent-accessible molecular surface, the computation of encounter rates between molecules in solution, the computation of the free energy of association and its salt dependence, the study of pKa shifts and the combination with classical molecular mechanics and dynamics. Theoretical results may be used for rationalizing or predicting experimental results, or for suggesting working hypotheses. An ever-increasing body of successful applications proves that the Poisson-Boltzmann equation is a useful tool for structural biology and complementary to other established experimental and theoretical methodologies.     

Thermodynamics of nucleosomal core particles

In this paper, a numerically detailed thermodynamic investigation of nucleosomal core particles is presented. The nonlinear Poisson-Boltzmann equation governs the electrostatic properties of both the DNA and histone protein. Brownian dynamics is used as the leading method, in combination with the analysis of the electrical features of the nucleosome. At elevated temperature, the structure of the nucleosome is destabilized by the decrease in electrical interactions of DNA-histone complexes, which can be explained with the EDL theory. Two obvious unwrapping transitions can be found, occurring within the temperature ranges 43-52 and 65-80 degrees C. The first transition is characterized by the melting of DNA terminal domains, and the feature of the second transition is the massive unwrapping of the DNA middle domain. It can be concluded that the nucleosomal DNA consists of two distinct structures, where the DNA terminal domains are less tightly bound to the histone than the DNA middle domain.     

Comparison of polyelectrolyte theories of the binding of cations to DNA.

Predictions of the binding of counterions to DNA made using the counterion condensation theory developed by Manning are compared with those made using the Poisson-Boltzmann equation, solved numerically by the Runge-Kutta procedure. Ions are defined as territorially or atmospherically bound if they fall within a given distance, defined by counterion condensation theory, from the DNA surface. Two types of experimental situations are considered. The first is the delocalized binding of a single type of counterion to DNA. In this case the Poisson-Boltzmann treatment predicts somewhat lower extents of binding TO DNA, modeled as a 10-A radius cylinder, than does Manning theory. The two theories converge as the radius decreases. The second type of experiment is the competition of ions of different valence for binding to DNA. The theories are compared with literature values of binding constants of divalent ions in the presence of monovalent ions, and of spermidine 3+ in the presence of Na+ or Mg2+. Both predict with fair accuracy the salt dependence of the equilibrium constants.     

Comparison of hard-cylinder and screened coulomb interactions in the modeling of supercoiled DNAs

A 1000 base pair (bp) model supercoiled DNA is simulated using spherical screened Coulomb interactions between subunits on one hand and equivalent hard-cylinder interactions on the other. The amplitudes, or effective charges, of the spherical screened Coulomb electrostatic potentials are chosen so that the electrostatic potential surrounding the middle of a linear array of 2001 subunits (31.8 Å diameter) closely matches the solution of the nonlinear Poisson-Boltzmann equation for a cylinder with 12 Å radius and the full linear charge density of DNA at all distances beyond the 24 Å hard-core diameter. This superposition of spherical screened Coulomb potentials is practically identical to the particular solution of the cylindrical linearized Poisson-Boltzmann equation that matches the solution of the nonlinear Poisson-Boltzmann equation at large distances. The interaction energy between subunits is reckoned from the effective charges according to the standard DLVO expression. The equivalent hard-cylinder diameter is chosen following Stigter's protocol for matching second virial coefficients, but for the full linear charge density of DNA. The electrostatic persistence length of the model with screened Coulomb interactions is extremely sensitive to the (arbitrarily) chosen subunit length at the higher salt concentrations. The persistence length of the hard-cylinder model is adjusted to match that of the screened Coulomb model for each ionic condition. Simulations for a superhelix density σ = -0.05 using a spherical screened Coulomb interaction plus a 24 Å hard-cylinder core (SCPHC) potential indicate that the radius of gyration of this 1000 bp DNA actually undergoes a slight increase as the NaCl concentration is raised from 0.01 to 1.0M. Thus, merely softening the potential from hard-cylinder to screened Coulomb form does not produce a large decrease in radius of gyration with increasing NaCl concentration for DNAs of this size. Radii of gyration, static structure factors, and diffusion coefficients obtained using the equivalent hard-cylinder (EHC) potential agree well with those obtained using the SCPHC potential in 1.0M NaCl, but in 0.1M NaCl the agreement is not as good, and in 0.01M NaCl the agreement is definitely unsatisfactory. These conclusions differ in significant respects from those obtained in previous studies. © 1997 John Wiley & Sons, Inc. Biopoly 42: 455–470, 1997     

A theory of the interaction between electrical double layers

A generalization of the Poisson-Boltzmann approach to the repulsive electrostatic force between similar electrical double layers is presented. It is based on the integral equation formalism of the statistical mechanical theories of fluids and it is shown that the Poisson-Boltzmann result follows from a well-defined and improvable set of approximations.     

Nucleic acid helix stability: effects of salt concentration, cation valence and size, and chain length

Metal ions play crucial roles in thermal stability and folding kinetics of nucleic acids. For ions (especially multivalent ions) in the close vicinity of nucleic acid surface, interion correlations and ion-binding mode fluctuations may be important. Poisson-Boltzmann theory ignores these effects whereas the recently developed tightly bound ion (TBI) theory explicitly accounts for these effects. Extensive experimental data demonstrate that the TBI theory gives improved predictions for multivalent ions (e.g., Mg2+) than the Poisson-Boltzmann theory. In this study, we use the TBI theory to investigate how the metal ions affect the folding stability of B-DNA helices. We quantitatively evaluate the effects of ion concentration, ion size and valence, and helix length on the helix stability. Moreover, we derive practically useful analytical formulas for the thermodynamic parameters as functions of finite helix length, ion type, and ion concentration. We find that the helix stability is additive for high ion concentration and long helix and nonadditive for low ion concentration and short helix. All these results are tested against and supported by extensive experimental data.