Decoupling of the nernst-planck and poisson equations. Application to a membrane system at overlimiting currents

This paper deals with one-dimensional stationary Nernst-Planck and Poisson (NPP) equations describing ion electrodiffusion in multicomponent solution/electrode or ion-conductive membrane systems. A general method for resolving ordinary and singularly perturbed problems with these equations is developed. This method is based on the decoupling of NPP equations that results in deduction of an equation containing only the terms with different powers of the electrical field and its derivatives. Then, the solution of this equation, analytical in several cases or numerical, is substituted into the Nernst-Planck equations for calculating the concentration profile for each ion present in the system. Different ionic species are grouped in valency classes that allows one to reduce the dimension of the original set of equations and leads to a relatively easy treatment of multi-ion systems. When applying the method developed, the main attention is paid to ion transfer at limiting and overlimiting currents, where a significant deviation from local electroneutrality occurs. The boundary conditions and different approximations are examined: the local electroneutrality (LEN) assumption and the original assumption of quasi-uniform distribution of the space charge density (QCD). The relations between the ion fluxes at limiting and overlimiting currents are discussed. In particular, attention is paid to the "exaltation" of counterion transfer toward an ion-exchange membrane by co-ion flux leaking through the membrane or generated at the membrane/solution interface. The structure of the multi-ion concentration field in a depleted diffusion boundary layer (DBL) near an ion-exchange membrane at overlimiting currents is analyzed. The presence of salt ions and hydrogen and hydroxyl ions generated in the course of the water "splitting" reaction is considered. The thickness of the DBL and its different zones, as functions of applied current density, are found by fitting experimental current-voltage curves.     

Simul 6: A fast dynamic simulator of electromigration

Simul 6 is a 1D dynamic simulator of electromigration based on the mathematical model of electromigration in free solutions. The model consists of continuity equations for the movement of electrolytes in a separation channel, acid–base equilibria of weak electrolytes, and the electroneutrality condition. It accounts for any number of multivalent electrolytes or ampholytes and provides a complete picture about dynamics of electromigration and diffusion in the separation channel. The equations are solved numerically using software means which allow for parallelization and multithreaded computation. Simul 6 has a user-friendly graphical interface. It is typically used for inspection of system peaks (zones) in electrophoresis, stacking and preconcentrating analytes, optimization of separation conditions, method development in either capillary zone electrophoresis, isotachophoresis, and isoelectric focusing. Simul 6 is the successor of Simul 5, and has been launched as a free software available for download at https://simul6.app/ .     

The electroneutrality approximation in electrochemistry

The electroneutrality approximation assumes that charge separation is impossible in electrolytic solutions. It has a long and successful history dating back to 1889 and may be justified because of the small absolute values for the permittivities of typical solvents. Dimensional analysis shows that the approximation becomes invalid only at nanosecond and nanometre scales. Recent work, however, has taken advantage of the capabilities of modern numerical simulation in order to relax this approximation, with concomitant advantages such as avoiding paradoxes and permitting a clear and consistent ??physical picture?? to describe charge dynamics in solution. These new theoretical techniques have been applied to liquid junction potentials and weakly supported voltammetry, with strong experimental corroboration for the latter. So long as dynamic processes are being studied, for which analytical solutions are unavailable in any case, numerical simulation is shown to render electroneutrality unnecessary as an a priori assumption.     

Mathematical model of electromigration allowing the deviation from electroneutrality

The structure of the double layer on the boundary between solid and liquid phases is described by various models, of which the Stern–Gouy–Chapman model is still commonly accepted. Generally, the solid phase is charged, which also causes the distribution of the electric charge in the adjacent diffuse layer in the liquid phase. We propose a new mathematical model of electromigration considering the high deviation from electroneutrality in the diffuse layer of the double layer when the liquid phase is composed of solution of weak multivalent electrolytes of any valence and of any complexity. The mathematical model joins together the Poisson equation, the continuity equation for electric charge, the mass continuity equations, and the modified G-function. The model is able to calculate the volume charge density, electric potential, and concentration profiles of all ionic forms of all electrolytes in the diffuse part of the double layer, which consequently enables to calculate conductivity, pH, and deviation from electroneutrality. The model can easily be implemented into the numerical simulation software such as Comsol. Its outcome is demonstrated by the numerical simulation of the double layer composed of a charged silica surface and an adjacent liquid solution composed of weak multivalent electrolytes. The validity of the model is not limited only to the diffuse part of the double layer but is valid for electromigration of electrolytes in general.     

The stationary phase method in the uniform WKB theory of inelastic scattering

When the modified Magnus approximation to the equation of motion in the uniform WKB theory of inelastic scattering is combined with the method of stationary phase for certain integrals, excellent numerical agreement is obtained between the transition probabilities thus calculated and the corresponding transition probabilities by the numerical solutions of the equation of motion and equivalent Schrödinger equations. The computation time required for calculation of transition probabilities is an order of magnitude less than that for numerical solution of Schrödinger equation by the Runge-Kutta-Gill method.     

'Objective' determination of coupling constants and conformational equilibria by solvent variation: a re-evaluation

An 'objective' method for determining conformational equilibria in substituted ethanes, proposed by Abraham et al., has been evaluated by computational methods. Abraham's method involves measuring vicinal couplings, such as (3)J(H,H) and (3)J(H,F), between methine and methylene protons with methine, methylene protons and fluorine in a range of solvents, on the assumption that the underlying coupling constants of the individual conformers are constant, but the fractions of each conformer in each solvent are different and unknown. Abraham posited that this would produce an 'over-determined' data set with more equations than unknowns would. Abraham's procedure is re-evaluated, and it is demonstrated that the type of system being considered here, where there are more equations than unknowns, is not necessarily over-determined. A computer equation solver and Monte Carlo-type procedures were employed to demonstrate that multiple numerical solutions exist for a representative 'over-determined' data set provided by Abraham. A statistical method was also developed to determine precisely which parameter sets constitute plausible solutions.     

On the nature of liquid junction and membrane potentials

Whenever a spatially inhomogeneous electrolyte, composed of ions with different mobilities, is allowed to diffuse, charge separation and an electric potential difference is created. Such potential differences across very thin membranes (e.g. biomembranes) are often interpreted using the steady state Goldman equation, which is usually derived by assuming a spatially constant electric field. Through the fundamental Poisson equation of electrostatics, this implies the absence of free charge density that must provide the source of any such field. A similarly paradoxical situation is encountered for thick membranes (e.g. in ion-selective electrodes) for which the diffusion potential is normally interpreted using the Henderson equation. Standard derivations of the Henderson equation appeal to local electroneutrality, which is also incompatible with sources of electric fields, as these require separated charges. We analyse self-consistent solutions of the Nernst-Planck-Poisson equations for a 1 : 1-univalent electrolyte to show that the Goldman and Henderson steady-state membrane potentials are artefacts of extraneous charges created in the reservoirs of electrolyte solution on either side of the membrane, due to the unphysical nature of the usual (Dirichlet) boundary conditions assumed to apply at the membrane-electrolyte interfaces. We also show, with the aid of numerical simulations, that a transient electric potential difference develops in any confined, but initially non-uniform, electrolyte solution. This potential difference ultimately decays to zero in the real steady state of the electrolyte, which corresponds to thermodynamic equilibrium. We explain the surprising fact that such transient potential differences are well described by the Henderson equation by using a computer algebra system to extend previous steady-state singular perturbation theories to the time-dependent case. Our work therefore accounts for the success of the Henderson equation in analysing experimental liquid-junction potentials.     

Applying Tikhonov regularization to process pendant droplet tensiometry data

The problem of obtaining the first and second derivatives of the profile of a pendant droplet is formulated as an integral equation of the first kind. This equation is solved by Tikhonov regularization in which the method of general cross validation is used to guide the selection of the regularization parameter. These derivatives are converted into mean curvature as a function of droplet height. Surface tension is then obtained by regression computation between the mean curvature and two possible algebraic expressions suggested by the Laplace-Young equation. This way of obtaining surface tension is demonstrated by applying it to a number of published droplet profiles. Some of the problems encountered are discussed and solutions suggested.     

Validity of the electroneutrality and goldman constant-field assumptions in describing the diffusion potential for ternary electrolyte systems in simple,porous membranes

Three numerical algorithms capable of simulating transport processes through simple, porous membranes in the steady state have been employed in order to study the change in the diffusion potential with the membrane thickness and the ionic concentrations for the ternary systems NaClHClH20 and CaCI2NaC1H₂O. The first simulation procedure uses Poisson's equation, the two others replace this equation by the electroneutrality and Goldman constant-field approximations respectively. From the results presented here, conditions for the applicability of the electroneutrality and constantfield assumption to ternary electrolyte systems are given.     

Numerical study of electroosmotic slip flow of fractional Oldroyd-B fluids at high zeta potentials

In this paper, an investigation of the electroosmotic flow of fractional Oldroyd-B fluids in a narrow circular tube with high zeta potential is presented. The Navier linear slip law at the walls is considered. The potential field is applied along the walls described by the nonlinear Poisson–Boltzmann equation. It's worth noting here that the linear Debye–Hückel approximation can't be used at the condition of high zeta potential and the exact solution of potential in cylindrical coordinates can't be obtained. Therefore, the Matlab bvp4c solver method and the finite difference method are employed to numerically solve the nonlinear Poisson–Boltzmann equation and the governing equations of the velocity distribution, respectively. To verify the validity of our numerical approach, a comparison has been made with the previous work in the case of low zeta potential and the excellent agreement between the solutions is clear. Then, in view of the obtained numerical solution for the velocity distribution, the numerical solutions of the flow rate and the shear stress are derived. Furthermore, based on numerical analysis, the influence of pertinent parameters on the potential distribution and the generation of flow is presented graphically.     

Finite difference method of simulation of non-steady-state ion transfer in electrochemical systems with allowance for migration

Finite difference methods of the second order of accuracy are elaborated for numerical calculation of non-steady-state ion transfer, which is caused by diffusion, migration, and convection in the unidimensional electrochemical systems. The methods of decoupling a set of coupled continuity equations of the electrolyte species are proposed, which ensures that the discrete equations are consistent with the initial differential equations and the electroneutrality condition is rigorously met. The methods of approximation of the boundary conditions of the second order temporal and spatial accuracy and the method of decoupling the transfer equations in the boundary nodes are elaborated. The explicit, fully implicit, and semi-implicit finite difference schemes are elaborated. For semi-implicit schemes, two versions of difference equation closure are proposed, which assure the unambiguity of determination of the distribution of electrical potential. Comparison analysis of the accuracy of elaborated finite difference methods of calculation of non-steady-state ion transfer is performed.     

Multigrid solution of the Poisson—Boltzmann equation

A multigrid method is presented for the numerical solution of the linearized Poisson–Boltzmann equation arising in molecular biophysics. The equation is discretized with the finite volume method, and the numerical solution of the discrete equations is accomplished with multiple grid techniques originally developed for twodimensional interface problems occurring in reactor physics. A detailed analysis of the resulting method is presented for several computer architectures, including comparisons to diagonally scaled CG, ICCG, vectorized ICCG and MICCG, and to SOR provided with an optimal relaxation parameter. Our results indicate that the multigrid method is superior to the preconditioned CG methods and SOR and that the advantage of multigrid grows with the problem size. © 1993 John Wiley & Sons, Inc.     

Ionic Ingress and Charge‐Neutralization Phenomena of Conducting‐Polymer Films

Ionic ingress and diffusion through a conducting‐polymer (CP) film containing embedded charges under potential and concentration gradients is studied. Electroneutrality, a common assumption in modeling of similar systems, is not justified in this case or similar diffusion‐limited processes, as the timescale of ionic diffusion in the solid matrix is quite large. Counter ions therefore cannot move instantaneously for effective neutralization of excess charges. Poisson–Nernst–Planck (PNP) equations have to be solved for such cases without any simplifying assumption. Analytical solution shows the existence of a charge boundary layer, which limits and strongly influences the ionic flux. A general numerical method for solution is also developed for the dynamic modeling, analysis, and design of these types of systems. The numerical results are validated by comparison with analytical solutions as well as with some experimental results available in the literature. With this modeling framework, the basic features of controlled release of molecules across a CP film under an applied electrical potential can be explained quantitatively.     

On the dilemma of the use of the electroneutrality constraint in electrochemical calculations

A simple treatment is developed to demonstrate that invoking the electroneutrality assumption to solve electrochemical transport problems is equivalent to invoking the limit of the exact (Nernst–Planck–Poisson) treatment when certain charge-density-related terms become negligibly small.     

Expressions for evaluating the possibility of slip at the liquid-solid interface in open tube capillary electrochromatography

In this work, expressions are constructed and solved that describe the velocity field of electroosmotic flow (EOF) in open tube capillary electrochromatography (CEC) systems when the possibility of having unequal tangential velocities at the liquid-solid interface is considered and a slip condition is employed as a boundary condition for the velocity of the EOF at the capillary wall. The coupled equations of hydrodynamics (momentum balance equation) and electrostatics (Poisson equation) are solved numerically in order to obtain the distribution of the velocity field as well as the value of the volumetric flow rate in the open tube. Also, expressions for the velocity field and the volumetric flow rate of the EOF are presented that are valid for certain electrolytic systems and for certain parameter values for which analytical solutions to the momentum balance and Poisson equations could be obtained. The results presented in this work indicate that having slip in the velocity of the EOF at the wall of the capillary could (i) substantially increase the electroosmotic velocity in the plug-flow region of the radial domain of the open capillary tube and (ii) increase the portion of the radial domain of the open capillary tube where the velocity of the EOF has a plug-flow profile, which in turn could increase the average velocity and volumetric flow rate of the EOF in the open capillary tube. Furthermore, the modeling approach and the results presented in this work indicate a method for experimentally evaluating the possibility of having slip in the velocity of the EOF at the capillary wall.     

Numerical Modeling of Non-Steady-State Ion Transfer in Electrochemical Systems with Allowance for Migration

Numerical methods that are used for modeling non-steady-state ion transfer in electrochemical systems and account for the diffusion, migration, convection, and homogeneous chemical reactions are analyzed. It is shown that the violation of the electroneutrality condition (ENC) in the process of numerical solution is due to the difference equations being inconsistent with the initial differential equations. Difference schemes for numerical calculation of transfer processes, which make it possible to split a set of coupled equations, are designed and conditions for their stability are determined. The explicit difference scheme is self-consistent, i.e. it ensures that ENC is rigorously met. In the implicit difference scheme, ENC is probably violated when splitting the set of equations. To restore electroneutrality of the medium, it is proposed to use a physically substantiated analytical relation for the space charge relaxation under the action of a strong electric field.     

A finite-difference method for numerical solution of the steady-state nernst—planck equations with non-zero convection and electric current density

A computer algorithm has been developed for digital simulation of ionic transport through membranes obeying the Nernst—Planck and Poisson equations. The method of computation is quite general and allows the treatment of steady-state electrodiffusion equations for multiionic environments, the ionic species having arbitrary valences and mobilities, when convection and electric current are involved. The procedure provides a great flexibility in the choice of suitable boundary conditions and avoids numerical instabilities which are so frequent in numerical methods. Numerical results for concentration and electric potential gradient profiles are presented in the particular case of the ternary system NaClHClH₂O.     

Accurate solutions for the spiked oscillators

The one‐dimensional harmonic oscillator potential plus a term of the form λ / x^α is known as the spiked oscillator potential; it constitutes a very interesting system because of its difficulties to accept perturbative and variational solutions for certain regimes of the α parameter. By the use of a numerical method, we obtain accurate energy eigenvalues and eigenfunctions for a wide range of λ values and a few α values. The accuracy of the present results is by much higher than previously published results. © 2008 Wiley Periodicals, Inc. Int J Quantum Chem, 2009     

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.     

Stochastic theory of large-scale enzyme-reaction networks: finite copy number corrections to rate equation models

Chemical reactions inside cells occur in compartment volumes in the range of atto- to femtoliters. Physiological concentrations realized in such small volumes imply low copy numbers of interacting molecules with the consequence of considerable fluctuations in the concentrations. In contrast, rate equation models are based on the implicit assumption of infinitely large numbers of interacting molecules, or equivalently, that reactions occur in infinite volumes at constant macroscopic concentrations. In this article we compute the finite-volume corrections (or equivalently the finite copy number corrections) to the solutions of the rate equations for chemical reaction networks composed of arbitrarily large numbers of enzyme-catalyzed reactions which are confined inside a small subcellular compartment. This is achieved by applying a mesoscopic version of the quasisteady-state assumption to the exact Fokker-Planck equation associated with the Poisson representation of the chemical master equation. The procedure yields impressively simple and compact expressions for the finite-volume corrections. We prove that the predictions of the rate equations will always underestimate the actual steady-state substrate concentrations for an enzyme-reaction network confined in a small volume. In particular we show that the finite-volume corrections increase with decreasing subcellular volume, decreasing Michaelis-Menten constants, and increasing enzyme saturation. The magnitude of the corrections depends sensitively on the topology of the network. The predictions of the theory are shown to be in excellent agreement with stochastic simulations for two types of networks typically associated with protein methylation and metabolism.