Simulation of the diffusion-controlled reaction between superoxide and superoxide dismutase. II. Detailed models

A two-stage Brownian dynamics simulation method is used to study the diffusion-influenced bimolecular reaction between superoxide and superoxide dismutase (SOD). The crystal structure of the dimeric enzyme is used in constructing detailed topographical and electrostatic models. Several electrostatic models are considered. In the most realistic, the excluded volume of the protein, which is impermeable to penetration by mobile ions, is assigned a dielectric constant of 2 and the surrounding “solvent” is assigned a value of 78. A finite difference method is used to solve the linearized Poisson–Boltzmann equation. For native SOD, the simulations reproduce the pronounced salt dependence of the rate constant observed experimentally. This salt dependence is attributed to electrostatic interactions between enzyme and substrate that are inherently attractive and amplified by the low dielectric constant of the protein interior. The simulation method is also applied to a modified enzyme, acylated SOD.     

Finite element analysis of the Poisson–Boltzmann equation coupled with chemical equilibriums: redistribution and transport of protons in nanophase separated polymeric acid–base proton exchange membranes

Abstract

The finite element analysis is applied to the study of the redistribution and transport of protons in model nanophase separated polymeric acid–base composite membranes by the Poisson–Boltzmann equation coupled with the acid and base dissociation equilibriums for the first time. Space charge redistribution in terms of proton and hydroxide redistributions is observed at the interfaces of acidic and basic domains. The space charge redistribution causes internal electrostatic potential, and thus, promotes the macroscopic transport of protons in the acid–base composite membranes.     

Hydration effects on the electrostatic potential around tuftsin

The electrostatic potential and component dielectric constants from molecular dynamics (MD) trajectories of tuftsin, a tetrapeptide with the amino acid sequence Thr–Lys–Pro–Arg in water and in saline solution are presented. The results obtained from the analysis of the MD trajectories for the total electrostatic potential at points on a grid using the Ewald technique are compared with the solution to the Poisson–Boltzmann (PB) equation. The latter was solved using several sets of dielectric constant parameters. The effects of structural averaging on the PB results were also considered. Solute conformational mobility in simulations gives rise to an electrostatic potential map around the solute dominated by the solute monopole (or lowest order multipole). The detailed spatial variation of the electrostatic potential on the molecular surface brought about by the compounded effects of the distribution of water and ions close to the peptide, solvent mobility, and solute conformational mobility are not qualitatively reproducible from a reparametrization of the input solute and solvent dielectric constants to the PB equation for a single structure or for structurally averaged PB calculations. Nevertheless, by fitting the PB to the MD electrostatic potential surfaces with the dielectric constants as fitting parameters, we found that the values that give the best fit are the values calculated from the MD trajectories. Implications of using such field calculations on the design of tuftsin peptide analogues are discussed. © 1999 John Wiley & Sons, Inc. Biopoly 50: 133–143, 1999     

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.     

Electroosmotic Flow in Nanoscale Parallel-plate Channels: Molecular Simulation Study and Comparison with Classical Poisson–Boltzmann Theory

Molecular dynamics simulations have been carried out for simple electrolyte systems to study the electrokinetically driven osmotic flow in parallel-plate channels of widths ～10–120?nm. The results are compared with the classical theory predictions based on the solution to the Poisson–Boltzmann equation. We find that despite some of the limitations in the Poisson–Boltzmann equation, such as assumption of the Boltzmann distribution for the ions, the classical theory captures the general trend of the variations of the osmotic flow with channel width, as characterized by the mobility of the fluid in channels between ～10 and 120?nm at moderate to low ion concentration. At moderate concentration (corresponding to relatively low surface potential), the classical theory is almost quantitative. The theory and simulation show more disagreement at low concentration, primarily caused by the high surface potential where the assumption of Boltzmann distribution becomes inaccurate. We discuss the limitations of the Poisson–Boltzmann equation as applied to the nanoscale channels.     

Electrostatic potential and Born energy of charged molecules interacting with phospholipid membranes: calculation via 3-D numerical solution of the full Poisson equation.

Understanding the physicochemical basis of the interaction of molecules with lipid bilayers is fundamental to membrane biology. In this study, a new, three-dimensional numerical solution of the full Poisson equation including local dielectric variation is developed using finite difference techniques in order to model electrostatic interactions of charged molecules with a non-uniform dielectric. This solution is used to describe the electric field and electrostatic potential profile of a charged molecule interacting with a phospholipid bilayer in a manner consistent with the known composition and structure of the membrane. Furthermore, the Born interaction energy is then calculated by appropriate integration of the electric field over whole space. Numerical computations indicate that the electrostatic potential profile surrounding a charge molecule and its resultant Born interaction energy are a function of molecular position within the membrane and change most significantly within the polar region of the bilayer. The maximum interaction energy is observed when the charge is placed at the center of the hydrophobic core of the membrane and is strongly dependent on the size of the charge and on the thickness of the hydrocarbon core of the bilayer. The numerical results of this continuum model are compared with various analytical approximations for the Born energy including models established for discontinuous slab dielectrics. The calculated energies agree with the well-known Born analytical expression only when the charge is located near the center of a hydrocarbon core of greater than 60 A in thickness. The Born-image model shows excellent agreement with the numerical results only when modified to include an appropriate effective thickness of the low dielectric region. In addition, a newly derived approximation which considers the local mean dielectric provides a simple and continuous solution that also agrees well with the numerical results.     

Applicability of a continuum solvation model to the octanol water transfer: CFF91 based model for amino acids

A continuum hydration model based upon the atomic charges provided with the CFF91 force field [A. B. Schmidt and R. M. Fine (1994) Molecular Simulation, 13. 347–365] has been extended to the octanol–water transfer. The electrostatic component of the transfer free energy is calculated using the finite-difference solution to the Poisson–Boltzmann equation while the nonpolar contributions are assumed to be proportional to the solute-excluded volume in water. All atomic charges and radii besides the aromatic carbon radius are equal in both solvents. The octanol dielectric constant and the probe radius are the main fitting parameters defining the octanol phase. The model has been tested for 38 organic molecules related to the amino acid residues and generally provides a high accuracy. In particular, the mean unsigned error for N-acetyl amino acid amides is 0.5 kcal/mol. © 1995 John Wiley & Sons, Inc.     

Modeling the electrophoresis of rigid polyions: application to lysozyme.

An algorithm is developed to determine the electrophoretic mobility of a rigid polyion modeled as a low dielectric volume element of arbitrary shape containing an arbitrary charge distribution. The solvent is modeled as a high dielectric continuum with salt distributed according to the linearized Poisson Boltzmann equation. Account is also taken of a Stern layer that separates the molecular surface and the surface of hydrodynamic shear, or Stern surface. Relaxation of the ion atmosphere because of the presence of the external field is ignored. The electrostatic and hydrodynamic problems are both solved by boundary element methods. The procedure is first applied to spherical polyions containing monopolar, dipolar, and quadrupolar charge distributions, and calculated mobilities are found to be in excellent agreement with the theory of Yoon and Kim. It is then applied to lysozyme by using models that account for the detailed shape and charge distribution of the enzyme. For reasonable choices of the molecular and Stern surfaces, calculated and experimental mobilities are found to be in fair agreement with each other. However, if a pH independent Stern layer (or, equivalently, translational diffusion constant, Dt) is assumed, the calculated mobilities exhibit a stronger pH dependence than is observed experimentally. A small increase in Dt with increasing pH could correct this discrepancy.     

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     

Monte Carlo and Poisson–Boltzmann calculations of the fraction of counterions bound to DNA

The counterion density and the condensation region around DNA have been examined as functions of both ion size and added-salt concentration using Metropolis Monte Carlo (MC) and Poisson–Boltzmann (PB) methods. Two different definitions of the “bound” and “free” components of the electrolyte ion atmosphere were used to compare these approaches. First, calculation of the ion density in different spatial regions around the polyelectrolyte molecule indicates, in agreement with previous work, that the PB equation does not predict an invariance of the surface concentration of counterions as electrolyte is added to the system. Further, the PB equation underestimates the counterion concentration at the DNA surface, compared to the MC results, the difference being greatest in the grooves, where ionic concentrations are highest. If counterions within a fixed radius of the helical axis are considered to be bound, then the fraction of polyelectrolyte charge neutralized by counterions would be predicted to increase as the bulk electrolyte concentration increases. A second categorization—one in which monovalent cations in regions where the average electrostatic potential is ledd than ?kT are considered to be bound—provides an informative basis for comparison of MC and PB with each other and with counterion-condensation theory. By this criterion, PB calculations on the B from of DNA indicate that the amount of bound counterion charge per phosphate group is about .67 and is independent of salt concentration. A particularly provocative observatiob is that when this binding criterion is used, MC calculations quantitatively reproduce the bound fraction predicated by counterion-condensation theory for all-atom models of B-DNA and A-DNA as well as for charged cylindera of varying lineat charge densities. For example, for B-DNA and A-DNA, the fractions of phosphate groups neutralized by 2 Å hard sphere counterions are 0.768 and .817, respectively. For theoretical studies, the rediys enclosing the region in which the electrostatic potential is calculated studies, the radius enclosing the region in which the electrostatic potential is calculated to be less than ?kT is advocated s a more suitable binding or condensation radius that enclosing the fraction of counterions given by (1 – ξ^?1). A comparsion of radii calculated using both of these definitions is presented. © 1994 John Wiley & Sons, Inc.     

Focusing of electric fields in the active site of Cu-Zn superoxide dismutase: effects of ionic strength and amino-acid modification

In this paper we report the implementation of a finite-difference algorithm which solves the linearized Poisson-Boltzmann equation for molecules of arbitrary shape and charge distribution and which includes the screening effects of electrolytes. The microcoding of the algorithm on an ST-100 array processor allows us to obtain electrostatic potential maps in and around a protein, including the effects of ionic strength, in about 30 minutes. We have applied the algorithm to a dimer of the protein Cu-Zn superoxide dismutase (SOD) and compared our results to those obtained from uniform dielectric models based on coulombic potentials. We find that both the shape of the protein-solvent boundary and the ionic strength of the solvent have a profound effect on the potentials in the solvent. For the case of SOD, the cluster of positive charge at the bottom of the active site channel produces a strongly enhanced positive potential due to the focusing of field lines in the channel-a result that cannot be obtained with any uniform dielectric model. The remainder of the protein is surrounded by a weak negative potential. The electrostatic potential of the enzyme seems designed to provide a large cross-sectional area for productive collisions. Based on the ionic strength dependence of the size of the positive potential region emanating from the active site and the repulsive negative potential barrier surrounding the protein, we are able to suggest an explanation for the ionic strength dependence of the activity of the native and chemically modified forms of the enzyme.     

Equilibrium surface charge distribution in phospholipid vesicles. I. Method of calculation

This paper presents a method of calculation of the surface charge equilibrium distribution between the two surfaces of a spherically closed phospholipid bilayer suspended in aqueous electrolyte solution. The net surface charge is supposed to be provided by the ionized polar groups of the phospholipid molecules. Its equilibrium distribution is found by minimization of the free electrostatic energy. The procedure of minimization utilizes the solution of the Poisson-Boltzmann equation which describes the double electric layers of the membrane and an expression for the membrane potential derived under the assumption of absence of charges in the membrane phase. An analytical solution of the problem in the range of validity of the linearized Poisson-Boltzman equation is obtained. It is shown that in this case an equilibrium transmembrane potential exists, and the surface charge density is greater at the outer surface of the vesicle.     

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.     

A precise analytical method for calculating the electrostatic energy of macromolecules in aqueous solution

A new method for calculating the total electrostatic free energy of a macromolecule in solution is presented. It is applicable to molecules of arbitrary shape and size, including membranes or macromolecular assemblies with substrate molecules and ions. The method is derived from integrating the energy density of the electrostatic field and is termed the field energy method. It is based on the dielectric model, in which the solute and the surrounding water are regarded as different continuous dielectrics. The field energy method yields both the interaction energy between all charge pairs and the self energy of single charges, effectively accounting for the interaction with water. First, the dielectric boundary and mirror charges are determined for all charges of the solute. The energy is then given as a simple function of the interatomic distances, and the standard atomic partial charges and volumes. The interaction and self energy are shown to result from three-body and pairwise interactions. Both energy terms explicitly involve apolar atoms, revealing that apolar groups are also subject to electrostatic forces. We applied the field energy method to a spherical model protein. Comparison with the Kirkwood solution shows that errors are within a small percentage. As a further test, the field energy method was used to calculate the electrostatic potential of the protein superoxide dismutase. We obtained good agreement with the result from a program that implements the numerical finite difference algorithm. The field energy method provides a basis for energy minimization and dynamics programs that account for the solvent and screening effect of water at little computational expense.     

On the theory of diffusion of electrolytes

In the theory of diffusion of electrolytes the following assumptions are frequently made: (i) the electrolytic solution is electrically neutral everywhere, (ii) the ionic concentrations and the electric potential all depend on a single Cartesian coordinate as the only space variable. Often the electric potential of the solution is determined on the basis of the Poisson equation alone, disregarding any other relation between this potential and the ionic concentrations. Since the Poisson equation only represents a condition which the potential fulfills, the use of this equation alone may lead to error unless the explicit relation for the potential involving a space integration of ionic concentrations is also taken into account. But if this relation is used the Poisson equation becomes redundant and, more important, assumptions (i) and (ii) appear unacceptable, the former because it leads to a zero electric potential everywhere, the latter because it is mathematically incorrect. The present paper is based on general equations of diffusion of ions, excluding the Poisson equation. These equations form a system of nonlinear integrodifferential quations whose number equals the number of ionic species present in the solution. It appears that when all ions are distributed symmetrically around a point all functions related to the above system of equations can be made dependent on a single space coordinate: the distance from the center of symmetry. Two methods of successive approximations are given for the solution of the equations in the case of spherical symmetry with limitation to the steady state. These methods are then applied to the study of the distribution of ionic concentrations and electrical potentials inside a cell of spherical shape in equilibrium with its surroundings. These methods are rapidly convergent; exact theoretical values of the electric potential are calculable on the boundary of the cell. It appears that the potential at the center of the cell is not more than ∼50% higher than at its boundary and that variation of concentration inside the cell is not very large. For instance, with 100 mV on the boundary the ionic concentration there is about four times higher than at the center. Calculations show that extremely small amounts of electricity are sufficient to account for the electric potentials currently observed. In a cell of 100 micra diameter an average concentration of only 10⁻¹⁴ mole/cm³ of a monovalent ion would be sufficient to give 1 millivolt on the boundary. This concentration is directly proportional to the voltage and inversely proportional to the square of the cell diameter. Most of the numerical results given above are obtained by considering only those ions whose electrical charge is not compensated for by ions of an opposite sign. The total concentrations may be much higher than those quoted. The theory does not take into account possible effects of structural heterogeneities which may exist in the cell, particularly of various phase boundaries. An incidental result shows that the Boltzmann distribution function in the form employed in modern theory of electrolytes is fundamentally a consequence of the mathematical theory of diffusion alone. It is pointed out, however, that Boltzmann distribution is not always compatible with the definition of the electric potential.     

Modeling the free solution electrophoretic mobility of short DNA fragments

Boundary element methods are used to model the free solution electrophoretic mobility of short DNA fragments. The Stern surfaces of the DNA fragments are modeled as plated cylinders that reproduce translational and rotational diffusion constants. The solvent-accessible and ion-accessible surfaces are taken to be coincident with the Stern surface. The mobilities are computed by solving simultaneously the coupled Navier–Stokes, Poisson, and ion-transport equations. The equilibrium electrostatics are treated at the level of the full Poisson–Boltzmann equation and ion relaxation is included. For polyions as highly charged as short DNA fragments, ion relaxation is substantial. At .11 M KCl, the simulated mobilities of a 20 base pair DNA fragment are in excellent agreement with experiment. At .04 M Tris acetate, pH = 8.0, the simulated mobilities are about 10–15% higher than experimental values and this discrepancy is attributed to the relatively large size of the Tris counterion. The length dependence of the mobility at .11 M KCl is also investigated. Earlier mobility studies on lysozyme are reexamined in view of the present findings. In addition to electrophoretic mobilities, the effective polyion charge measured in steady state electrophoresis and its relationship to the preferential interaction parameter γgG is briefly considered. © 1998 John Wiley & Sons, Inc. Biopoly 46: 359–373, 1998     

Reduced point charge models of proteins: assessment based on molecular dynamics simulations

A reduced point charge distribution is used to model Ubiquitin and two complexes, Vps27 UIM-1–Ubiquitin and Barnase–Barstar. It is designed from local extrema in charge density distributions obtained from the Poisson equation applied to smoothed molecular electrostatic potentials. A variant distribution is built by locating point charges on atoms. Various charge fitting conditions are selected, i.e. from either electrostatic Amber99 (Assisted Model Building with Energy Refinement) Coulomb potential or forces, considering reference grid points located within various distances from the protein atoms, with or without separate treatment of main and side chain charges. The program GROMACS (Groningen Machine for Chemical Simulations) is used to generate Amber99SB molecular dynamics (MD) trajectories of the solvated proteins modelled using the various reduced point charge models (RPCMs) so obtained. Point charges that are not located on atoms are considered as virtual sites. Some RPCMs lead to stable MD trajectories. They, however, involve a partial loss in the protein secondary structure and lead to a less-structured solute solvation shell. The model built by fitting charges on Coulomb forces calculated at grid points ranging between 1.4 and 2.0 times the van der Waals radius of the atoms, with a separate treatment of main chain and side chain charges, appears to best approximate all-atom MD trajectories.     

Modeling the electrophoresis of lysozyme. II. Inclusion of ion relaxation.

In this work, boundary element methods are used to model the electrophoretic mobility of lysozyme over the pH range 2-6. The model treats the protein as a rigid body of arbitrary shape and charge distribution derived from the crystal structure. Extending earlier studies, the present work treats the equilibrium electrostatic potential at the level of the full Poisson-Boltzmann (PB) equation and accounts for ion relaxation. This is achieved by solving simultaneously the Poisson, ion transport, and Navier-Stokes equations by an iterative boundary element procedure. Treating the equilibrium electrostatics at the level of the full rather than the linear PB equation, but leaving relaxation out, does improve agreement between experimental and simulated mobilities, including ion relaxation improves it even more. The effects of nonlinear electrostatics and ion relaxation are greatest at low pH, where the net charge on lysozyme is greatest. In the absence of relaxation, a linear dependence of mobility and average polyion surface potential, (lambda zero)s, is observed, and the mobility is well described by the equation [formula: see text] where epsilon 0 is the dielectric constant of the solvent, and eta is the solvent viscosity. This breaks down, however, when ion relaxation is included and the mobility is less than predicted by the above equation. Whether or not ion relaxation is included, the mobility is found to be fairly insensitive to the charge distribution within the lysozyme model or the internal dielectric constant.     

Influence of the membrane potential on the free energy of an intrinsic protein

A modified Poisson-Boltzmann equation is developed from statistical mechanical considerations to describe the influence of the transmembrane potential on macromolecular systems. Using a Green's function formalism, the electrostatic free energy of a protein associated with the membrane is expressed as the sum of three terms: a contribution from the energy required to charge the system's capacitance, a contribution corresponding to the interaction of the protein charges with the membrane potential, and a contribution corresponding to a voltage-independent reaction field free energy. The membrane potential, which is due to the polarization interface, is calculated in the absence of the protein charges, whereas the reaction field is calculated in the absence of transmembrane potential. Variations in the capacitive energy associated with typical molecular processes are negligible under physiological conditions. The formulation of the theory is closely related to standard algorithms used to solve the Poisson-Boltzmann equation and only small modifications to current source codes are required for its implementation. The theory is illustrated by examining the voltage-dependent membrane insertion of a simple polyalanine alpha-helix and by computing the electrostatic potential across a 60-A-diameter sphere meant to represent a large intrinsic protein.