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.     

Adaptive finite element simulation of chronoamperometry at microdisc electrodes

In this paper we shall introduce a transient finite element algorithm by considering the simplest problem of numerical simulation of the chronoamperometric current for an E reaction mechanism at a microdisc electrode. Such a numerical simulation is made difficult by the presence of a boundary singularity where the electrode meets the insulator (often known as the ‘edge-effect’) and by a time-singularity caused by the impulsive start of the experiment. We attempt to overcome these problems by using an adaptive finite element algorithm in which we derive an a posteriori bound on the error in the computed current. This is used to guide mesh refinement and adaptive time-stepping, resulting in a fully automated algorithm which is both accurate and efficient.     

Adaptative finite element simulation of currents at microelectrodes to a guaranteed accuracy. First-order EC′ mechanism at inlaid and recessed discs

In this paper we extend the introductory work described in the accompanying papers on the use of adaptive finite element methods in electrochemical simulation (K. Harriman et al., Electrochem. Commun. 2 (2000) 150 and 157) to the case of a (pseudo) first-order EC′ reaction mechanism at both an inlaid and a recessed disc. The recessed disc is shown to be a particularly suitable example for illustrating the power of the technique in providing the simulated current to a guaranteed accuracy on near-optimal meshes. For both problems we demonstrate that we can obtain excellent accuracy across the spectrum of reaction rates using just a few seconds of CPU time. Our results also confirm the accuracy of some recently published analytical solutions to these problems.     

Adaptive finite element methods in electrochemistry

In this article, we review some of our previous work that considers the general problem of numerical simulation of the currents at microelectrodes using an adaptive finite element approach. Microelectrodes typically consist of an electrode embedded (or recessed) in an insulating material. For all such electrodes, numerical simulation is made difficult by the presence of a boundary singularity at the electrode edge (where the electrode meets the insulator), manifested by the large increase in the current density at this point, often referred to as the edge effect. Our approach to overcoming this problem has involved the derivation of an a posteriori bound on the error in the numerical approximation for the current that can be used to drive an adaptive mesh-generation algorithm, allowing calculation of the quantity of interest (the current) to within a prescribed tolerance. We illustrate the generic applicability of the approach by considering a broad range of steady-state applications of the technique.     

Computation of dynamic adsorption with adaptive integral,finite difference,and finite element methods

Analysis of diffusion-controlled adsorption and surface tension in one-dimensional planar coordinates with a finite diffusion length and a nonlinear isotherm, such as the Langmuir or Frumkin isotherm, requires numerical solution of the governing equations. This paper presents three numerical methods for solving this problem. First, the often-used integral (I) method with the trapezoidal rule approximation is improved by implementing a technique for error estimation and choosing time-step sizes adaptively. Next, an improved finite difference (FD) method and a new finite element (FE) method are developed. Both methods incorporate (a). an algorithm for generating spatially stretched grids and (b). a predictor-corrector method with adaptive time integration. The analytical solution of the problem for a linear dynamic isotherm (Henry isotherm) is used to validate the numerical solutions. Solutions for the Langmuir and Frumkin isotherms obtained using the I, FD, and FE methods are compared with regard to accuracy and efficiency. The results show that to attain the same accuracy, the FE method is the most efficient of the three methods used.     

Application of Lobatto polynomials in atomic physics

The approximation properties of Lobatto polynomials are analyzed and then applied to approximate the atomic Kohn–Sham eigenfunction. In the first part of this article the approximation algorithm based on the Galerkin finite element method is derived. To obtain the approximation of the function, based on the presented algorithm, the linear set of equations must be solved. The matrix of the equation set is very sparse and its elements can be evaluated analytically. In the second part of this article, the algorithm is applied to evaluate adaptive polynomial approximation of selected Kohn–Sham eigenstates of indium (In) atom. The proposed r‐adaptive algorithm evaluates the minimum number of subintervals needed to represent the eigenfunction with required accuracy. Based on the r‐adaptive algorithm, the approximations of 4d, 5s 5p In eigenfunctions were calculated. © 2009 Wiley Periodicals, Inc. Int J Quantum Chem, 2010     

Adaptive multilevel finite element solution of the Poisson–Boltzmann equation I. Algorithms and examples

This article is the first of two articles on the adaptive multilevel finite element treatment of the nonlinear Poisson–Boltzmann equation (PBE), a nonlinear eliptic equation arising in biomolecular modeling. Fast and accurate numerical solution of the PBE is usually difficult to accomplish, due to the presence of discontinuous coefficients, delta functions, three spatial dimensions, unbounded domain, and rapid (exponential) nonlinearity. In this first article, we explain how adaptive multilevel finite element methods can be used to obtain extremely accurate solutions to the PBE with very modest computational resources, and we present some illustrative examples using two well‐known test problems. The PBE is first discretized with piece‐wise linear finite elements over a very coarse simplex triangulation of the domain. The resulting nonlinear algebraic equations are solved with global inexact Newton methods, which we have described in an article appearing previously in this journal. A posteriori error estimates are then computed from this discrete solution, which then drives a simplex subdivision algorithm for performing adaptive mesh refinement. The discretize–solve–estimate–refine procedure is then repeated, until a nearly uniform solution quality is obtained. The sequence of unstructured meshes is used to apply multilevel methods in conjunction with global inexact Newton methods, so that the cost of solving the nonlinear algebraic equations at each step approaches optimal O(N) linear complexity. All of the numerical procedures are implemented in MANIFOLD CODE (MC), a computer program designed and built by the first author over several years at Caltech and UC San Diego. MC is designed to solve a very general class of nonlinear elliptic equations on complicated domains in two and three dimensions. We describe some of the key features of MC, and give a detailed analysis of its performance for two model PBE problems, with comparisons to the alternative methods. It is shown that the best available uniform mesh‐based finite difference or box‐method algorithms, including multilevel methods, require substantially more time to reach a target PBE solution accuracy than the adaptive multilevel methods in MC. In the second article, we develop an error estimator based on geometric solvent accessibility, and present a series of detailed numerical experiments for several complex biomolecules. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 1319–1342, 2000     

A new strategy for simulation of electrochemical mechanisms involving acute reaction fronts in solution: Application to model mechanisms

The new strategy for adaptive simulation of electrochemical reaction mechanisms described in a previous paper (C. Amatore, O. Klymenko, I. Svir. Electrochem. Commun., (2010), doi:10.1016/j.elecom.2010.06.009) [1] provides an efficient method of obtaining accurate concentration distributions and electrochemical currents. In this paper, this strategy is illustrated and tested more deeply upon simulating several representative classical electrochemical mechanisms involving fast homogeneous comproportionation or disproportionation reactions that pose severe difficulties when simulated by classical finite difference methods including those based on exponentially expanding grids.     

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.     

Chronoamperometric current at finite disk electrodes

The chronoamperometric current at a stationary finite disk electrode is studied using both analytical and digital simulation techniques. The exact long-time expansion of the current is obtained and its short-time behavior is considered. Digital simulation of the current using an explicit hopscotch algorithm is presented. In contrast to the usual explicit difference method, the ‘hopscotch’ algorithm is unconditionally stable, and thus, it is particularly suited for studying electrochemical problems at intermediate and long times. A simple analytic expressions for the current, which is accurate to 0.6% for all times, is presented.     

A parallel finite element simulator for ion transport through three‐dimensional ion channel systems

A parallel finite element simulator, ichannel, is developed for ion transport through three‐dimensional ion channel systems that consist of protein and membrane. The coordinates of heavy atoms of the protein are taken from the Protein Data Bank and the membrane is represented as a slab. The simulator contains two components: a parallel adaptive finite element solver for a set of Poisson–Nernst–Planck (PNP) equations that describe the electrodiffusion process of ion transport, and a mesh generation tool chain for ion channel systems, which is an essential component for the finite element computations. The finite element method has advantages in modeling irregular geometries and complex boundary conditions. We have built a tool chain to get the surface and volume mesh for ion channel systems, which consists of a set of mesh generation tools. The adaptive finite element solver in our simulator is implemented using the parallel adaptive finite element package Parallel Hierarchical Grid (PHG) developed by one of the authors, which provides the capability of doing large scale parallel computations with high parallel efficiency and the flexibility of choosing high order elements to achieve high order accuracy. The simulator is applied to a real transmembrane protein, the gramicidin A (gA) channel protein, to calculate the electrostatic potential, ion concentrations and I – V curve, with which both primitive and transformed PNP equations are studied and their numerical performances are compared. To further validate the method, we also apply the simulator to two other ion channel systems, the voltage dependent anion channel (VDAC) and α‐Hemolysin (α‐HL). The simulation results agree well with Brownian dynamics (BD) simulation results and experimental results. Moreover, because ionic finite size effects can be included in PNP model now, we also perform simulations using a size‐modified PNP (SMPNP) model on VDAC and α‐HL. It is shown that the size effects in SMPNP can effectively lead to reduced current in the channel, and the results are closer to BD simulation results. © 2013 Wiley Periodicals, Inc.     

Numerical solutions for isoelectric focusing and isotachophoresis problems

This study combines an adaptive mesh redistribution (AMR) method and the space–time conservation element and solution element (CESE) method to construct a high-resolution scheme for the solution of electrophoresis pre-concentration and separation problems. In the proposed AMR–CESE scheme, the fine mesh points are moved toward the regions of discontinuity within the solution domain in accordance with the equidistribution principle. To reduce the numerical dissipation within the regions of the solution domain with a large spatial mesh, the spatial component of the CESE scheme is treated using a Courant–Friedrichs–Lewy (CFL) number insensitive scheme. The validity of the proposed approach is confirmed by comparing the results obtained for typical isoelectric focusing (IEF) and isotachophoresis (ITP) problems with those obtained from the conventional CESE scheme and the finite volume method (FVM), respectively. It is shown that the AMR–CESE scheme yields a better accuracy than uniform fixed-mesh solvers with no more than a minor increase in the computational cost.     

Adaptative finite element simulation of currents at microelectrodes to a guaranteed accuracy. Application to a simple model problem

In this series of papers we consider the general problem of numerical simulation of the currents at microelectrodes using an adaptive finite element approach. Microelectrodes typically consist of an electrode embedded (or recessed) in an insulating material. For all such electrodes, numerical simulation is made difficult by the presence of a boundary singularity at the electrode edge (where the electrode meets the insulator), manifested by the large increase in the current density at this point, often referred to as the ‘edge-effect’. Our approach to overcoming this problem involves the derivation of an a posteriori bound on the error in the numerical approximation for the current that can be used to drive an adaptive mesh-generation algorithm. This allows us to calculate the current to within a prescribed tolerance. We begin by demonstrating the power of the method for a simple model problem — an E reaction mechanism at a microdisc electrode — for which the analytical solution is known. In this paper we give the background to the problem, and show how an a posteriori error bound can be used to drive an adaptive mesh-generation algorithm. We then use the algorithm to solve our model problem and obtain very accurate results on comparatively coarse meshes in minimal computing time. We give the technical details of the background theory and the derivation of the error bound in the accompanying paper.     

Adaptive multilevel finite element solution of the Poisson–Boltzmann equation II. Refinement at solvent‐accessible surfaces in biomolecular systems

We apply the adaptive multilevel finite element techniques (Holst, Baker, and Wang 21 ) to the nonlinear Poisson–Boltzmann equation (PBE) in the context of biomolecules. Fast and accurate numerical solution of the PBE in this setting is usually difficult to accomplish due to presence of discontinuous coefficients, delta functions, three spatial dimensions, unbounded domains, and rapid (exponential) nonlinearity. However, these adaptive techniques have shown substantial improvement in solution time over conventional uniform‐mesh finite difference methods. One important aspect of the adaptive multilevel finite element method is the robust a posteriori error estimators necessary to drive the adaptive refinement routines. This article discusses the choice of solvent accessibility for a posteriori error estimation of PBE solutions and the implementation of such routines in the “Adaptive Poisson–Boltzmann Solver” (APBS) software package based on the “Manifold Code” (MC) libraries. Results are shown for the application of this method to several biomolecular systems. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 1343–1352, 2000     

High-resolution modeling of isotachophoresis and zone electrophoresis

The space-time conservation element and solution element (CESE) method is applied to simulate the ITP and zone electrophoresis (ZE) separation phenomena. The CESE method expresses the governing equation in the integral form of the conservation law, and has a second-order accuracy in both space and time. The current results show that the CESE solutions for the ITP and ZE phenomena are more accurate than those obtained using conventional numerical schemes, which are characterized by serious numerical diffusion and oscillation. Furthermore, the CESE method suppresses the numerical oscillations or peaks observed in the results obtained using traditional second-order finite difference schemes. Finally, the results reveal that the CESE method accurately models the sharp boundaries between adjacent ITP samples under steady-state conditions. Overall, the results presented in this study demonstrate the numerical accuracy of the CESE method and confirm its applicability to the modeling of a range of electrophoretic phenomena.     

Filling factor of a paramagnetic sample in a rectangular cavity: theory and application

A computational method is presented for calculating the filling factor of an electron paramagnetic resonance (EPR) tube in a rectangular TE102 cavity. The algorithm employs the conventional finite element method. In addition to the filling factor, the algorithm allows to calculate the quality factor and the reflection coefficient of the loaded cavity. This method allows calculating very accurately the EPR signal intensities from which the spin concentration of paramagnetic samples can be determined. A comparison between the predicted EPR signal intensities to several experimental results was found to be satisfactory. The method also allows optimizing the EPR tube dimensions and its glass quality to improve measurement sensitivity.     

Electrochemical impedance spectroscopy investigations of a microelectrode behavior in a thin-layer cell: Experimental and theoretical studies

Electrochemical impedance spectroscopy experiments were performed on a microdisk electrode in a thin-layer cell using a scanning electrochemical microscope for controlling the cell geometry. Experimental data showed that when the thin-layer thickness diminished, an additional low-frequency response appeared. It was ascribed to the radial diffusion of the electroactive species and was strongly dependent on the thin-layer dimensions (both thickness and diameter). Moreover, the numerical simulation of the impedance diagrams by finite element method calculations confirmed this behavior. An equivalent circuit based on a Randles-type circuit was proposed. Thus, the diffusion was described by introducing two electrical elements: one for the spherical diffusion and the other for the radial contribution. A nonlinear Simplex algorithm was used, and this circuit was shown to fit the impedance diagrams with a good accuracy.     

Computational Electrochemistry: A Model to Studying Ohmic Distortion of Voltammetry in Multiple Working Electrode,Microfluidic Devices,an Adaptive FEM Approach

Characteristics of mass transport and potential distribution applicable to microfluidic electrochemical flow cell devices has been modelled using the finite element method. A flexible, automatic grid generation algorithm has been combined with an a‐posteriori error indication technique presented by Nann and Heinze to allow irregular cell geometries to be modelled. The code has been applied to the problem of steady state generator – detector linear sweep voltammetry in a channel flow cell showing the effects of IR drop on the voltammetric response of each electrode.     

A Brownian dynamics-finite element method for simulating DNA electrophoresis in nonhomogeneous electric fields

The objective of this work is to develop a numerical method to simulate DNA electrophoresis in complicated geometries. The proposed numerical scheme is composed of three parts: (1) a bead-spring Brownian dynamics (BD) simulation, (2) an iterative solver-enhanced finite element method (FEM) for the electric field, and (3) the connection algorithm between FEM and BD. A target-induced searching algorithm is developed to quickly address the electric field in the complex geometry which is discretized into unstructured finite element meshes. We also develop a method to use the hard-sphere interaction algorithm proposed by Heyes and Melrose [J. Non-Newtonian Fluid Mech. 46, 1 (1993)] in FEM. To verify the accuracy of our numerical schemes, our method is applied to the problem of lambda-DNA deformation around an isolated cylindrical obstacle for which the analytical solution of the electric field is available and experimental data exist. We compare our schemes with an analytical approach and there is a good agreement between the two. We expect that the present numerical method will be useful for the design of future microfluidic devices to stretch and/or separate DNA.