Question

Theoretical Battery Science Question:

You have determined that ignoring convection is not valid for a new battery chemistry that you have invented. Write out how you would modify the Poisson Nernst Planck system of equations.

Answer 1

The conventional Poisson-Nernst-Planck equations do not account for the finite size of ions explicitly. This leads to solutions featuring unrealistically high ionic concentrations in the regions subject to external potentials, in particular, near highly charged surfaces. A modified form of the Poisson-Nernst-Planck equations accounts for steric effects and results in solutions with finite ion concentrations. Here, we evaluate numerical methods for solving the modified Poisson-Nernst-Planck equations by modeling electric field-driven transport of ions through a nanopore. We describe a novel, robust finite element solver that combines the applications of the Newton's method to the nonlinear Galerkin form of the equations, augmented with stabilization terms to appropriately handle the drift-diffusion processes.

To make direct comparison with particle-based simulations possible, our method is specifically designed to produce solutions under periodic boundary conditions and to conserve the number of ions in the solution domain. We test our finite element solver on a set of challenging numerical experiments that include calculations of the ion distribution in a volume confined between two charged plates, calculations of the ionic current though a nanopore subject to an external electric field, and modeling the effect of a DNA molecule on the ion concentration and nanopore current.

Beginning with the experiments that revealed the microscopic mechanisms of nerve cell excitation [1], measurements of ion currents through nanoscale channels and pores have become the basis of many experimental techniques in biology and biotechnology. In addition to permitting the study of the behavior of individual proteins that allow the passage of ions into and out of cells [2], ion current measurements through nanopores have been used to study the rupture of molecular bonds [3–5], to distinguish between similar molecules [6], and to determine the properties and sequences of nucleic acid molecules [7–11]. However, since direct experimental imaging of molecules within nano-pores is extremely difficult, computation plays an important role in associating current with nanoscale phenomenon [12–19] (see [20, 21] for recent reviews of the field).

Equilibrium and transport properties of ionic solutions can be simulated using explicit ion methods such as all-atom molecular dynamics [16, 20] or Brownian dynamics [22–24], or by using continuum models such as the Poisson-Boltzmann and Poisson-Nernst-Planck equations [25, 26]. While the explicit ion methods provide the most accurate description of the system's behavior, both in spatial and temporal domains, they are stochastic in nature and thus require long, computationally expensive simulations to obtain average properties. Furthermore, the application of an explicit ion method usually requires the system to be described with the same resolution over the entire simulation domain. Often, this leads to a situation where a majority of the computational effort is applied to simulate a nearly uniform solution where quantities of interest exhibit little variation. In contrast, continuum methods allow different regions of the same system to be described at varying levels of detail, and thus focus the computational effort on regions that require a more precise description. In addition to being more computationally efficient, continuum models more easily incorporate certain types of boundary conditions that arise in physical systems, such as boundaries of fixed concentration or electrostatic potential.

The traditional continuum approach to modeling ionic transport is based on the Poisson-Nernst-Planck equations (PNPE). Although the PNPE have been applied successfully to model the electro-diffusion phenomena [27, 28], the equations are not without drawbacks. Within the PNPE approach, ions are modeled as mathematical points of negligible physical dimension, thereby allowing for accumulation of ions at unrealistically high concentrations in certain regions of the system. A modified formulation of the PNPE, called the modified Poisson-Nernst-Planck equations (MPNPE) [29], explicitly takes the physical dimensions of ions into consideration, which limits the maximum concentration that attained in the system. The advantage of using MPNPE over PNPE becomes apparent in the systems that contain regions subject to strong attractive potentials, for example, near charged surfaces.

In this work, we explore the MPNPE approach for modeling equilibrium and transport properties of ionic solutions in realistic three-dimensional geometries subject to realistic applied potentials. The finite difference method has been widely used to solve the Nernst-Planck equations in one or three dimensions [27, 30–32]. Although the finite difference method is straightforward to implement, applying this method to systems that have curved boundaries and complicated geometries is challenging. In this respect, using a finite element method is more appropriate as it naturally handles complex geometries, such as the molecular surfaces of DNA molecules and ion channels. Finite element methods for solving the three-dimensional PNPE have already been described [33, 34]. However, numerical studies of the MPNPE have been limited to one-dimensional systems [29] or the three dimensional spherical case [35] and have not been applied to simulate ion flow through a solid-state nanopore, which is the main process considered in this work.

Here, we introduce a three-dimensional MPNPE solver for the simulation of ionic current through nanopores, which can handle the complex geometry of the system and the realistic microscopic potentials the ions are subject to. The nanopore system is illustrated in Fig. 2 and described in detail in Section 2. In contrast to the previous efforts, our finite element method conserves ion concentration, takes into account the sharp repulsive potentials present near the walls of an ion channel [13, 27], and is able to reproduce the results of explicit ions simulations. The presence of a sharp, repulsive potential at the interface of fluid and solid-state domains necessitates formulation of a new stable numerical method for finding the solution of the PNPE and the MPNPE. Specifically, we found that, when applied to our nanopore systems, standard finite element methods become unstable and produce spurious results such as negative concentrations. Fig. 1 gives an example of such behavior. The sharp repulsive potential near the walls of a nanopore causes instability of the Galerkin method, producing spurious negative concentration values (see Section 4, Experiment 3 for more details). Below, we describe a numerical procedure that stabilizes the finite element method in the presence of sharp repulsive potentials, which is one of the main results of this works.

We consider the Poisson-Nernst-Planck Equations (PNPE) for a 1:1 electrolyte solution (referred to as solvent) described over a computational domain, denoted by Ω = Ωs ∪ Ωm, which includes both the solvent region, represented as Ωs, as well as a molecular or membrane region, Ωm, which is void of solvent. The time dependent PNPE are given as [27]

∂c±∂t=D±∇⋅[∇c±+1kBT[±e(c±∇ϕ)+(c±∇U)]]inΩs,

(2.1)

−∇ ⋅ ∊∇ϕ = e(c+ − c−) in Ω,

(2.2)

where ϕ is the electrostatic potential and U is the potential due to other interactions (such as van der Waals and solvation forces), which is assumed to be the same for both ionic species. Hereafter, we will refer to potential U as a non-electrostatic potential, to differentiate it from the explicit electrostatic potential ϕ. In the Nernst-Planck equation, (2.1), the concentration of positive and negative ions are c+ and c–, respectively, kB is the Boltzmann's constant, T is temperature, e is the charge on an electron and D± are the diffusivities of the positive and negative ions, respectively. In the Poisson equation, (2.2), we assume a piecewise constant dielectric coefficient ε that is defined in the two sub-domains, Ωs and Ωm. For simplicity, we write the total potential energy experienced by an ion as V± =±eϕ+U.

The modified form of the PNPE (MPNPE) adds a nonlinear term to each of the two Nernst-Planck equations in (2.1) to model the steric repulsion. The Poisson equation remains unchanged, however the modified Nernst Planck equations are [29],

∂c±∂t=D±∇⋅[∇c±+1kBTc±∇V±+a3(c±∇(c++c−)1−c+a3−c−a3)]inΩs,

(2.3)

Here a is the size of the ion (assumed to be the same for both species). As a result, in this model the maximum permitted concentration is bounded by 1/a3, which we refer to as the steric limit. To simplify the presentation of the material that follows, we write the PNPE and the MPNPE as

∂c±∂t=D±∇⋅[∇c±+1kBTc±∇V±+Nα(c±)]inΩs,

(2.4)

where

Nα(c±)=α(a3c±∇(c++c−)1−c+a3−c−a3),α={0,1,forPNPE,forMPNPE.

(2.5)

A primary focus of this paper is the application of MPNPE solver to nanopores, wherein we compute the ionic current through a pore in a solid-state membrane. The domain we consider is depicted in Fig. 2(a). Here, solution reservoirs above and below the membrane are connected through a nanopore, allowing positive and negative ions to pass from one side of the membrane to the other. We also consider a system where a DNA molecule is present inside the pore. Thus, the membrane (and the DNA, if present) comprise the domain Ωm, whereas the ionic solution, which consists of the solution reservoirs above and below the membranes and the nanopore, comprise the domain Ωs.

With the concentration profiles of (2.4), one important quantity is the ionic current J through a surface G with normal n (see Fig. 2(b)), which is defined as

J=∑±±∫ΓeD±n⋅[∇c±+1kBT[±e(c±∇ϕ)+c±∇V]]ds.

(2.6)

For example, in the 2D cross-section of the problem domain shown in Fig. 2(b), we measure the ionic current through the plane in the middle of the pore, denoted by a dotted line.

The Poisson portion of the PNPE in (2.2) is solved with Dirichlet boundary conditions specified by ϕt and ϕb at the top and the bottom of the domain, and periodic boundary conditions along the other four sides. Further, the (unmodified and modified) Nernst-Planck equations in (2.1) and (2.3) use blocking boundary conditions on the interface of the membrane and the ionic solution, which is denoted ∂Ωs,n and is displayed with a dotted line in Fig. 2(b), while periodic boundary conditions are set at the remaining boundaries. Specifically, we consider blocking boundary conditions of the form

[∇c±+1kBTc±∇V±+Nα(c±)]⋅n=0on∂Ωs,n,

(2.7)

where n is the unit normal on the surface ∂Ωs,n. One consequence of the blocking boundary conditions is that the integral of the concentration remains constant. That is, the total number of ions of each ion species is conserved:

∂∂t∫Ωsc±dx=∫ΩsD±∇⋅[∇c±+1kBTc±∇V±+Nα(c±)]dx=∫∂Ωs,n[∇c±+1kBTc±∇V±+Nα(c±)]⋅nds=0.

(2.8)

Additionally, the number of ions in the domain is also conserved in the case of a partially periodic boundary. Moreover, we note that the numerical methods we develop are equally applicable for other boundary conditions, like Dirichlet boundary conditions for the concentrations.