Boundary element solution of the linear Poisson-Boltzmann equation and a multipole method for the rapid calculation of forces on macromolecules in solution

The Poisson-Boltzmann equation is widely used to describe the electrostatic potential of molecules in an ionic solution that is treated as a continuous dielectric medium. The linearized form of this equation, applicable to many biologic macromolecules, may be solved using the boundary element method. A single-layer formulation of the boundary element method, which yields simpler integral equations than the direct formulations previously discussed in the literature, is given. It is shown that the electrostatic force and torque on a molecule may be calculated using its boundary element representation and also the polarization charge for two rigid molecules may be rapidly calculated using a noniterative scheme. An algorithm based on a fast adaptive multipole method is introduced to further increase the speed of the calculation. This method is particularly suited for Brownian dynamics or molecular dynamics simulations of large molecules, in which the electrostatic forces must be calculated for many different relative positions and orientations of the molecules. It has been implemented as a set of programs in C++, which are used to study the accuracy and speed of this method for two actin monomers.     

Iterative Solution Method for the Linearized Poisson-Boltzmann Equation: Indirect Boundary Integral Equation Approach

An iterative solution scheme is proposed for solving the electrical double-layer interactions governed by the linearized Poisson-Boltzmann equation. The method is based on the indirect integral equation formulation with the double-layer potential kernel of the linearized Poisson-Boltzmann equation. In contrast to the conventional direct integral equation approach that yields Fredholm integral equations of the first kind, the indirect integral equation approach yields well-posed Fredholm integral equations of the second kind. The eigenvalue analysis reveals that the spectral radius of the double-layer integral operator is always less than one. Thus, iterative solution schemes can be successfully implemented for solving the electrical double-layer interactions for very large and complex systems. The utility of the iterative indirect method is demonstrated for several examples which include spherical and spheroidal particles. Copyright 2001 Academic Press.     

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.     

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

A hybrid approach for solving the nonlinear Poisson-Boltzmann equation (PBE) is presented. Under this approach, the electrostatic potential is separated into (1) a linear component satisfying the linear PBE and solved using a fast boundary element method and (2) a correction term accounting for nonlinear effects and optionally, the presence of an ion-exclusion layer. Because the correction potential contains no singularities (in particular, it is smooth at charge sites) it can be accurately and efficiently solved using a finite difference method. The motivation for and formulation of such a decomposition are presented together with the numerical method for calculating the linear and correction potentials. For comparison, we also develop an integral equation representation of the solution to the nonlinear PBE. When implemented upon regular lattice grids, the hybrid scheme is found to outperform the integral equation method when treating nonlinear PBE problems. Results are presented for a spherical cavity containing a central charge, where the objective is to compare computed 1D nonlinear PBE solutions against ones obtained with alternate numerical solution methods. This is followed by examination of the electrostatic properties of nucleic acid structures.     

Numerical modeling of Joule heating-induced temperature gradient focusing in microfluidic channels

Temperature gradient focusing (TGF) is a recently developed technique for spatially focusing and separating ionic analytes in microchannels. The temperature gradient required for TGF can be generated either by an imposed temperature gradient or by Joule heating resulting from an applied electric field that also drives the flow. In this study, a comprehensive numerical model describing the Joule heating induced temperature development and TGF is developed. The model consists of a set of governing equations including the Poisson-Boltzmann equation, the Laplace equation, the Navier-Stokes equations, the energy equations and the mass transport equation. As the thermophysical and electrical properties including the liquid dielectric constant, viscosity, and electric conductivity are temperature-dependent, these governing equations are coupled, and therefore the coupled governing equations are solved numerically by using a CFD-based numerical method. The numerical simulations agree well with the experimental results, suggesting the valid mathematical model presented in this study.     

DRBEM solution of biomagnetic fluid flow and heat transfer in cavities-CMMSE2016

In this paper, we investigate the fully developed, laminar, forced convection flow of an electrically non-conducting, viscous, biomagnetic fluid in the 2D cross-section (cavity) of a long impermeable pipe. The fluid is under the influence of a point magnetic source placed below the cavity. The dual reciprocity boundary element method (DRBEM) with constant and linear elements is used for solving the governing equations resulting from the Navier–Stokes and energy equations together with magnetization and buoyancy forces. The fundamental solution of Laplace equation is made use of converting differential equations to boundary integral equations by taking all the terms other than Laplacian as inhomogeneity in the Poisson’s equations for the velocity components, pressure and the temperature of the fluid. The unknown pressure boundary conditions are approximated through momentum equations by using finite difference approximation for the pressure gradients and DRBEM coordinate matrix for the other terms. All the space derivatives are also calculated by DRBEM coordinate matrix which is one of the main advantages of DRBEM. Pipe axis velocity is also computed. The effects of magnetization and the buoyancy force on the fluid with or without viscous dissipation term in the energy equation are investigated in square and lid-driven cavities for several values of magnetic (Mn) and Rayleigh (Ra) numbers. It is observed that the flow and heat transfer are significantly affected with increasing values of Mn and Ra. DRBEM gives small sized linear systems due to its boundary only nature at a considerably low computational expense.     

Polarizability and Kerr constant of proteins by boundary element methods

A precise implementation of the boundary element method has been applied to the computation of the polarizability and the Kerr constant of eight soluble proteins. The method is demonstrated to be accurate and precise by comparison with analytical values for spheroids. Two different integral equations have been solved: (1) an exact equation with explicit dielectric constant dependence, and (2) an exact equation for a metallic body. The dielectric dependence for the metallic body case is built in with a generalization of the ellipsoid formula. Both methods agree quantitatively with each other for low relative dielectric constants. A full tensor expression for the Kerr constant yields perfect agreement with experiment for some proteins and badly under reports for the rest. The protein structure is obtained from a crystallographic database and is assumed rigid throughout the TEB measurement. Electrolyte effects are neglected. The Kerr constant is dominated by the protein dipole moment and is quite sensitive to the orientation of the dipole moment relative to the principal axes of the polarizability tensor. Several possible reasons for the large discrepancy between some computed and measured values are discussed.     

The interaction energy between two parallel plates with constant surface charge density

On the basis of Langmuir's suggestion we simplify the Poisson-Boltzmann equation and derive the relation of surface potential, potential midway, and the plate distance. Thus we obtain the interaction force and energy equations between two dissimilar plates in the case of constant surface charge density. Agreement with the exact numerical values of the interaction of dissimilar plates is good. This method may not only apply to the cases of high constant potential but to the case of high constant charge density.     

The dependence of electrostatic solvation energy on dielectric constants in Poisson-Boltzmann calculations

The Poisson-Boltzmann equation gives the electrostatic free energy of a solute molecule (with dielectric constant epsilon(l)) solvated in a continuum solvent (with dielectric constant epsilon(s)). Here a simple formula is presented that accurately predicts the electrostatic free energy for all combinations of epsilon(l) and epsilon(s) from the calculation on a single set of epsilon(l) and epsilon(s) values.     

Electric fields in an electrolyte solution near a strip of fixed potential

Electrostatic fields produced by flat electrodes are often used to manipulate particles in solution. To study the field produced by such an electrode, we consider the problem of an infinite strip of width 2a with imposed constant potential immersed in an electrolyte solution. Sufficiently close to the edge of the strip, the solution is determined by classical electrostatics and results in a field singularity. We examine two limiting cases, (a) when strip width a<1k, the Debye screening length, and (b) when strip width is much greater than the Debye screening length, a>1k. We present exact results for the two cases in the limit of small potentials where the Poisson-Boltzmann equation can be linearized. By drawing on an analogy with antiplane shear deformations of solids, and by employing the path-independent J integral of solid mechanics, we present a new method for determining the strength of the edge singularity. The strength of the singularity defines an exact near-field solution. In the far field the solution goes to that of a line of charge. The accuracy of the solution is demonstrated by comparison with the numerical solutions of the Poisson-Boltzmann equation using the finite element method.     

Electroosmotic Flow in Microchannels

The electroosmotic flow induced by an applied electrostatic potential field through microchannels between two parallel plates and a 90 degrees bend is analyzed in this work. A nonlinear, two-dimensional Poisson-Boltzmann equation governing the electrical double-layer field and the Laplace equation governing the electrostatic field distribution in microchannels are numerically solved using a finite-difference method. A body force caused by the interaction between the electrical double-layer field and the applied electrostatic field is included in the full Navier-Stokes equations. The effects of the electrical double-layer field and the applied electrostatic field on the fluid velocity distribution, pressure drop, and skin friction are discussed. A small pressure drop along the parallel plates is detected, although it is always neglected in the literature. Pressure is not a constant across the channel height. The axial velocity profile is no longer flat across the channel height when the Reynolds number is large. A separation bubble is detected near the 90 degrees junction when the Reynolds number is large. Copyright 2001 Academic Press.     

Solution of the linearized Poisson-Boltzmann equation

Improved methods are formulated for solution of the linearized Poisson-Boltzmann equation, to be used in conjunction with electronic structure calculation on a solute together with dielectric continuum representation of the salt-containing solvent. Volume polarization effects due to quantum mechanical penetration of solute charge density outside the cavity that excludes solvent are treated by exact and by approximate methods analogous to those previously developed for the salt-free case. With boundary element approaches, exact solutions lead to coupled equations for a pair of cavity surface distributions that mimic the polarization of the solvent dielectric and the ionic atmosphere. A novel means is found to effectively decouple these equations, yielding more efficient practical methods for their numerical solution. Detailed comparisons are given to related boundary element formulations previously reported in the literature, which neglect volume polarization, and analogous decoupling is also found for the pair of surface distributions invoked there. Illustrative results are provided for a simple spherical example.     

Analysis of rotation-driven electrokinetic flow in microscale gap regions of rotating disk systems

In the present study, a novel theoretical model is developed for the analysis of rotating thermal-fluid flow characteristics in the presence of electrokinetic effects in the microscale gap region between two parallel disks under specified electrostatic, rotational, and thermal boundary conditions. The major flow configuration considered is a rotor-stator disk system. Axisymmetric Navier-Stokes equations with consideration of electric body force stemming from streaming potential are employed in the momentum balance. Variations of the fluid viscosity and permittivity with the local fluid temperature are considered. Between two disks, the axial distribution of the electric potential is determined by the Poisson equation with the concentration distributions of positive and negative ions obtained from Nernst-Planck equations for convection-diffusion of the ions in the flow field. Effects of disk rotation and electrostatic and thermal conditions on the electrokinetic flow and thermal characteristics are investigated. The electrohydrodynamic mechanisms are addressed with an interpretation of the coupling nature of the electric and flow fields. Finally, solutions with electric potential determined by employing nonlinear or linearized Poisson-Boltzmann equation and/or invoking assumptions of constant properties are compared with the predictions of the present model for justification of various levels of approximation in solution of the electrothermal flow behaviors in rotating microfluidic systems.     

Solution of the nonlinear Poisson-Boltzmann equation using pseudo-transient continuation and the finite element method

The nonlinear Poisson-Boltzmann (PB) equation is solved using Newton-Krylov iterations coupled with pseudo-transient continuation. The PB potential is used to compute the electrostatic energy and evaluate the force on a user-specified contour. The PB solver is embedded in a existing, 3D, massively parallel, unstructured-grid, finite element code. Either Dirichlet or mixed boundary conditions are allowed. The latter specifies surface charges, approximates far-field conditions, or linearizes conditions "regulating" the surface charge. Stability and robustness are proved using results for backward Euler differencing of diffusion equations. Potentials and energies of charged spheres and plates are computed and results compared to analysis. An approximation to the potential of the nonlinear, spherical charge is derived by combining two analytic formulae. The potential and force due to a conical probe interacting with a flat plate are computed for two types of boundary conditions: constant potential and constant charge. The second case is compared with direct force measurements by chemical force microscopy. The problem is highly nonlinear-surface potentials of the linear and nonlinear PB equations differ by over an order of magnitude. Comparison of the simulated and experimentally measured forces shows that approximately half of the surface carboxylic acid groups, of density 1/(0.2 nm2), ionize in the electrolyte implying surface charges of 0.4 C/m2, surface potentials of 0.27 V, and a force of 0.6 nN when the probe and plate are 8.7 nm apart.     

Radiative and nonradiative decay rates of a molecule close to a metal particle of complex shape

We present a model to evaluate the radiative and nonradiative lifetimes of electronic excited states of a molecule close to a metal particle of complex shape and, possibly, in the presence of a solvent. The molecule is treated quantum mechanically at Hartree-Fock (HF) or density-functional theory (DFT) level. The metal/solvent is considered as a continuous body, characterized by its frequency dependent local dielectric constant. For simple metal shapes (planar infinite surface and spherical particle) a version of the polarizable continuum model based on the integral equation formalism has been used, while an alternative methodology has been implemented to treat metal particles of arbitrary shape. In both cases, equations have been numerically solved using a boundary element method. Excitation energies and nonradiative decay rates due to the energy transfer from the molecule to the metal are evaluated exploiting the linear response theory (TDHF or TDDFT where TD--time dependent). The radiative decay rate of the whole system (molecule + metal/solvent) is calculated, still using a continuum model, in terms of the response of the surrounding to the molecular transition. The model presented has been applied to the study of the radiative and nonradiative lifetimes of a lissamine molecule in solution (water) and close to gold spherical nanoparticles of different radius. In addition, the influence of the metal shape has been analyzed by performing calculations on a system composed by a coumarin-type molecule close to silver aggregates of complex shape.     

Highly accurate biomolecular electrostatics in continuum dielectric environments

Implicit solvent models based on the Poisson-Boltzmann (PB) equation are frequently used to describe the interactions of a biomolecule with its dielectric continuum environment. A novel, highly accurate Poisson-Boltzmann solver is developed based on the matched interface and boundary (MIB) method, which rigorously enforces the continuity conditions of both the electrostatic potential and its flux at the molecular surface. The MIB based PB solver attains much better convergence rates as a function of mesh size compared to conventional finite difference and finite element based PB solvers. Consequently, highly accurate electrostatic potentials and solvation energies are obtained at coarse mesh sizes. In the context of biomolecular electrostatic calculations it is demonstrated that the MIB method generates substantially more accurate solutions of the PB equation than other established methods, thus providing a new level of reference values for such models. Initial results also indicate that the MIB method can significantly improve the quality of electrostatic surface potentials of biomolecules that are frequently used in the study of biomolecular interactions based on experimental structures.     

A quantitative evaluation of the electrostatic theory for ion pair chromatography

Summary A quantitative model for ion pair chromatography based on the electrostatic theory is described. The model is based on the solution of the linearised Poisson-Boltzmann equation in a cylinder. The obtained equations are compared with experimental data from a number of different systems. The agreement between theory and experiments is satisfactorily. Systematic deviations due to the use of the linearised equation and ion correlation effects are discussed.     

A Model for Calculating Electrostatic Interactions between Colloidal Particles of Arbitrary Surface Topology

A numerical model for calculating the electrostatic interaction between two particles of arbitrary shape and topology is described. A key feature of the model is a generalized discretization program, capable of simulating any desired analytical shape as a set of flat, triangular elements. The relative sizes of the elements are adjusted using a density function to better match the desired shape and the spatial variation of the electrical surface properties on each particle. The distribution of either surface potential or surface charge density is then calculated using a boundary element approach to solve the linearized Poisson-Boltzmann equation. Example interaction energy profiles are calculated for three different types of roughness-bumps, pits, and surface waves. It is found that the interaction energy between rough particles remains different from that between two equivalent smooth spheres at all separations, even for gap widths much larger than either the solution Debye length or the characteristic roughness size. This behavior at large gap widths arises from the nature of the decay of the electric potential away from each particle. In addition, the magnitude of the roughness effect is found to depend greatly on the size and shape of the nonuniformity as well as the electrostatic boundary conditions. For example, for a sphere containing asperities of height equal to 0.2 times the particle radius, the interaction energy can be as much as 50% greater than that between two equivalent spheres under the condition of constant surface potential. At constant surface charge density, the ratio of the interaction energies between rough and smooth spheres was found to either diverge or become zero as contact between the two particles is approached, depending on the nature of the roughness. Changes of this magnitude could clearly have a substantial impact on the stability behavior of a dispersion of such particles. Copyright 2001 Academic Press.     

Computation of electrostatic forces between solvated molecules determined by the Poisson-Boltzmann equation using a boundary element method

A rigorous approach is proposed to calculate the electrostatic forces among an arbitrary number of solvated molecules in ionic solution determined by the linearized Poisson-Boltzmann equation. The variational principle is used and implemented in the frame of a boundary element method (BEM). This approach does not require the calculation of the Maxwell stress tensor on the molecular surface, therefore it totally avoids the hypersingularity problem in the direct BEM whenever one needs to calculate the gradient of the surface potential or the stress tensor. This method provides an accurate and efficient way to calculate the full intermolecular electrostatic interaction energy and force, which could potentially be used in Brownian dynamics simulation of biomolecular association. The method has been tested on some simple cases to demonstrate its reliability and efficiency, and parts of the results are compared with analytical results and with those obtained by some known methods such as adaptive Poisson-Boltzmann solver.     

Calculation of the Maxwell stress tensor and the Poisson-Boltzmann force on a solvated molecular surface using hypersingular boundary integrals

The electrostatic interaction among molecules solvated in ionic solution is governed by the Poisson-Boltzmann equation (PBE). Here the hypersingular integral technique is used in a boundary element method (BEM) for the three-dimensional (3D) linear PBE to calculate the Maxwell stress tensor on the solvated molecular surface, and then the PB forces and torques can be obtained from the stress tensor. Compared with the variational method (also in a BEM frame) that we proposed recently, this method provides an even more efficient way to calculate the full intermolecular electrostatic interaction force, especially for macromolecular systems. Thus, it may be more suitable for the application of Brownian dynamics methods to study the dynamics of protein/protein docking as well as the assembly of large 3D architectures involving many diffusing subunits. The method has been tested on two simple cases to demonstrate its reliability and efficiency, and also compared with our previous variational method used in BEM.