The current-voltage relation of a pore and its asymptotic behavior in a Nernst-Planck model

In this article, we provide a simple explanation for the nonlinearity and the asymmetry of the current-voltage characteristic of a pore in a membrane, using a NernstPlanck (NP) model. The NP equation is difficult to solve analytically because it is nonlinear. The solution of the NP equation becomes easy to find assuming a linear profile for the ion's potential energy inside the pore, in this particular case a triangular asymmetric profile. The model we found for the electrical current through a pore is bolstered by similar features found experimentally by other investigators including current nonlinearity within two voltage domains of linear behavior (at low and high voltages) and current asymmetry (rectification) where ions flow more easily through a pore entrance with a mild energy gradient than through one with a steeper energy gradient.

The current-voltage relation of a pore is essential for studying ion channels [1, 2], or for designing synthetic nanopores or nanodevices [3, 4, 5, 6]. The current-voltage relation of a porous membrane is useful in understanding the way in which epithelia function, especially the skin's epidermal stratum corneum [7] or for describing and predicting electroporation [8, 9]. Thus, an improved understanding of the current-voltage relationship of a pore in a membrane is of use to various research fields of biology, material science, bioengineering, and nanoelectronics.

There are various approaches to modeling ion transport in membrane pores from molecular dynamics and Brownian dynamics to continuum theories [10]. A continuum model, like the one we are presenting herein, is an explicative model offering an increased understanding of the importance which various processes involved in the problem have, thus giving this type of model an increased predictive power. Such a model usually provides an analytical formula for the current-voltage relationship, a very useful feature for interpreting experimental data. The down side of the continuum model is related to the modeling accuracy of ion–channel interactions, because it works with macroscopic quantities, such as the dielectric constant, on a microscopic-mesoscopic scale [10, 11]. This drawback is to some degree compensated by the fact that the ion potential energy can be estimated from molecular or Brownian dynamics simulations.

We will first present the proposed model with its basic assumptions, followed by results for the current-voltage relations for the assumed profiles of ion energy inside the pore, and their behavior at low-level and high-level voltages. Finally, we will discuss the implications and present the related work that supports our results.

Our model assumes the existence of a cylindrical pore, with length h and radius R, in a plane membrane bathed by an electrolyte with only one monovalent ionic species. The generalization for multiple ionic species is straightforward [12, 13] but unnecessary here.

The transport equation states that the stationary flux, F, of ions through the pore in the membrane is given by the product of the ionic concentration, c, the mobility, b, and the gradient of the electrochemical potential, μ [13, 14, 15]:

The mobility, b, is related to the diffusion coefficient of the ion, D, by Einstein's relation, b=D/(kT), where k is the Boltzmann constant and T is absolute temperature. The electrochemical potential μ is:

where the first term is the concentration dependence (the ionic activity is assumed to be unity), the second is the contribution of the applied potential with elementary charge, q, the third is the work, W, required to transfer the ion from a distant point in the aqueous phase to point, x, inside the pore, and the last term is the standard chemical potential in the pore, independent of position.

The potential energy, W(x), of an ion inside the pore is the key element for pore conductance and will be treated here as an assumed parameter of the model. The potential energy has two distinct parts [10]: 1) the self-energy due to the induced charges on the dielectric boundary, surrounding liquid and pore wall (Born energy and reaction-field energy) proportional to q² which is always repulsive and 2) the Coulomb interaction with the charges on the pore's wall (if the pore's wall is charged), proportional to q, which can be either attractive or repulsive depending on the signs of interacting charges.

Substituting Eq. 2 into Eq. 1 we can represent the stationary current density, J, as:

where the dimensionless energy variables are:

In Eq. 3 the exponential factor in front of the derivative may be transferred to the left side of the equation to arrive at

Eq. 5 can then integrated with respect to x

to give the standard expression for stationary current density [10, 13]:

It is presupposed that the ionic concentration in the bulk solution and at the pore's extremities are equal: c(0)=c(h)=c.

The denominator can be integrated if the function relating energy to position is known. To this end, the profile of the potential energy inside the pore, W(x), is approximated here by an asymmetric triangle (Eq. 8), a symmetric triangle (Eq. 9), and a trapezium (Eq. 10) as in DeBruin and Krassowska (1999) [12,16], to produce the following formulations:

Here w is the dimensionless height of the energy barrier (the potential energy maximum in kT units) and d is the length of the pore entrance, or more specifically the length of the pore domain where the potential energy varies. It should be mentioned that the energy barrier height, w, and the pore radius, R, are tightly bound. The simplest relationship between them is w=5nm/R(nm) from Glaser et al. (1988) [16]. More complex formulas can be found in Kuyucak et al. (2001) [10]. This simple relation tells us that the model is appropriate for nanopores that are not too large (so that w→0) and not too narrow (to avoid ionic crowding or depletion, kinetic effects, steric effects, etc. [6,10, 11,14]).

Fig. 1

The externally applied potential V(x) is assumed to vary linearly across the membrane and along the pore as:

where u is the dimensionless transmembrane potential:

We introduce two dimensionless quantities: the relative size of the entrance region of the pore r, Eq. 13, and a nondimensional current density j, Eq. 14:

From these, we derive the formulas of the dimensionless current densities for the asymmetric triangular (j₁), symmetric triangular (j₂), and trapezoidal (j₃) energy profiles:

It is worth noting that without an energy barrier (w=0) the current density behavior is ohmic for all three formulas, the same as that of the bulk solution current density, j₀:

For the symmetric profiles, triangular and trapezoidal, the current density formulas are antisymmetric, i.e. j(u)=−j(−u). In other words, for these energy profiles there is no rectification. Only the asymmetric energy barrier generates asymmetry in the j-u curve as seen in Fig. 2. The rule here is that current will flow easier from the side with a smaller energy gradient to the side with a larger gradient, than in the opposite direction. This finding tells us about the importance of energy profile at the pore entrance.

Fig. 2

At large applied voltages, when |u|>>1, the behavior of the current density is quasilinear for all three formulas, as shown in Fig. 2. It is not a pure ohmic behavior, because the asymptotic line is shifted relative to the bulk current density, j₀, but it has the same slope as that of the bulk current density. This linear behavior occurs analytically, using the approximation e–|^u| ≈e–^r|^u| ≈ 0 in Eq. 15 and 17 for the asymmetric triangular and trapezoidal profiles. Thus the current density formulas for large voltages are

Equations 19 and 20 also apply for the symmetric triangular profile simply by using r=1/2. The term w/r or w/(1–r) represents a threshold voltage, u_T, which marks the region where the current density slope switches from a low value below u_T (|u|<u_T) to a large value above u_T (|u|>u_T).

Below |u|<1, the current densities of the two triangular profiles practically coincide, as shown in the double logarithmic plot of current density versus voltage of Fig. 3. The logarithmic slope of the three current densities is similar to that of the bulk current density, j₀, i.e. there is a quasi-ohmic, linear dependence between current and voltage at low voltages. This is analytically confirmed using a Taylor series expansion for |u|<<1 of the current density formulas. Neglecting higher order terms, we may represent the current density formulas for low voltages as

Fig. 3

In the case of small applied voltages, the current density has a quasi-ohmic behavior with a gentle slope, expressing a low conductance, high resistance behavior imposed by the energy barrier height, w, through the exponential term e^w.

For the "inverse current" region, u>0, the current density for the asymmetric triangular profile, j₁, coincides with the current density for the symmetric triangular profile, j₂, at low voltages (u<1) and for large voltages, above u_T, j₁ it is similar with the current density of the trapezoidal profile, j₃, as shown in Figs. 2 and 3.

For the j-u characteristic plotted here, there are two distinct regions of quasi-linear behavior: near the origin (|u|<1), a slow rising quasi-ohmic current density described by Eq. 21 or 22, and a fast rising linear, non-ohmic current density described by Eq. 19 or 20, when the potential exceeds a certain threshold (|u|>u_T). The transition between these two regimes is made in a smooth nonlinear fashion.

This Nernst-Planck model for the current density in a membrane's pore shows that the current density as a function of voltage has three distinct regions: 1) a low voltage region with a quasi-ohmic, slow rising linear current density whose slope depends mainly on the height of the energy barrier, w, through the exponential term e^w; 2) an intermediate voltage region with a nonlinear current-voltage dependence; and 3) a high voltage region with a non-ohmic, fast rising linear current whose slope is equal to that of the bulk current density, but shifted relative to it by a threshold voltage, u_T, which is also the potential energy gradient at the pore entrance, w/r.

Here we consider the j-u relationship quasi-ohmic at low voltages because the ratio j/u is constant and thus follows the standard Ohm's law. For the linear j-u relationship at high voltages the ratio j/u is not a constant and for this reason we consider it non-ohmic.

For voltages u>>1, the current density for the asymmetric triangular profile, j₁, coincides with the current density for the trapezoidal profile, j₃, as seen in Fig. 2. This fact indicates that the potential energy gradient at the pore entrance, w/r, controls the current density behavior and that the potential energy gradient at the pore exit has no influence upon the current density, as shown by Eqs. 19 and 20.

Such a behavior of the current is found in porated bilayer lipid membranes [17] and in protegrin transmembrane pores [18], for example. There are other experimental situations where these two linear regimes are clearly visible. In some lipid pores the I–V curve is linear from approximately −150 to +150 mV, and the I–V relation becomes non-linear when |V|>150 mV [19]. At intermediate pH values of the bathing solution, in polymeric-synthetic nanopores (double-conical pore) the I– V curve deviates from the linear behavior in a similar fashion with our model [4]. In a synthetic conical nanopore, the current-voltage curve shows a mild linear slope for "inverse current" and a steeper slope for "direct current", the milder slope for the "direct current" being less visible for the chosen voltage scale, but more visible for 20 mM KCl [20]. The I–V relation in cellular ion channels, usually asymmetric, has a similar aspect with ours [2, 21].

The energy profiles used here simplify calculations and are also good approximations for many practical problems. For example, a triangular profile for the energy barrier arises from a more elaborate calculation of the ion self-energy inside the pore [22, 23], and the asymmetrical triangular profile explains the ratchet mechanism in Brownian motors [5].

The effects of asymmetry are now used in nanofluidic diodes [3, 4], for example, and in nanoslits [24]. The energy asymmetry can be obtained from the geometrical asymmetry of the pore, like in synthetic conical pores [25] or intentional misalignment in nanoslits [24], and from surface charges due to manufacturing processes [3,20,26]. Similarly, in biological channels, this asymmetry can arise from the inherent geometrical asymmetry of the protein pore [1,18] and/or from the charged state of the protein pore or the phospholipids heads of the lipid membrane [18,27]. In addition, the charged state of the pore's wall can be trimmed by the surrounding solution pH [2] or by changing the surface chemical composition [20].

The asymmetry of the potential energy profile at the pore extremities is very important for the rectification effect to occur. Other researchers have demonstrated this fact with somewhat different theoretical approaches, using the Smoluchowski equation [28] and Poisson-Nernst-Planck equations [29]. Our approach has the advantage of providing a direct link between cause and effect, embedding the complicated part of the pore physics into the potential energy profile. Molecular or Brownian dynamics can estimate the potential energy profile inside the pore by computer simulation [10,26].

Finer details and a more rigorous analytical approach that includes kinetic effects, steric effects, large voltages, etc. can be achieved using the method employed by Bazant [11,14]. Different I-V curves that markedly differ from the one presented here are obtained in low concentration electrolytes or at large voltages [24], but such cases that imply ionic crowding or depletion are beyond the scope of this paper.

This simple model tells us that an uncharged conical pore will rectify equally well the ionic current for anions and cations because the self-energy is proportional to q². If the pore has a charged wall, the ionic current for anions and cations will differ markedly, because the Coulomb interaction is proportional to q. Therefore, the pore will present charge selectivity because the anions and the cations will sense a different potential energy inside the pore, a fact already found experimentally [20].

Another important practical problem is the current-voltage characteristic of the skin. The manner in which human skin responds to an externally applied electric field is of importance for several medical applications with either a diagnostic or treatment purpose [30]. The skin and particularly the outer stratum corneum is a porous membrane with a strong nonlinear response to applied voltages or currents [31, 32], which also has a measurable asymmetry [33, 34, 35].

Usually, Bode plots, Nyquist representations, Cole diagrams/formulas, or equivalent electrical circuits descriptively present the skin's electrical response. In order to understand the large changes in skin resistance at high voltages, Weaver and his collaborators have developed an explicative theoretical model for the skin electrical behavior [8], starting from electroporation studies in bilayer lipid membranes [16]. The electroporation phenomenon drastically changes the electrical behavior of the skin. The applied voltage modifies the number of pores in the skin, with a rate of pore creation that depends exponentially on the squared voltage [36].

Electroporation theory consistently describes skin behavior at large voltages, but for low voltages there is not yet an accepted theory for skin electrical nonlinearity and asymmetry. Current models for skin at low voltages have taken into account the appendageal macropores [36] and sweat pores [33], but for low voltage across the skin (<5V) where electroporation is supposed to be absent, however, the current-voltage characteristic of the skin has a strong similarity to our model [36]. In order to develop medical applications based on the electrical properties of human skin such as low voltage iontophoresis for drug delivery [37, 38], an in-depth understanding of these properties is required. Complex modeling of the skin [7,8,36] can be improved by taking into account the single pore effects described by our model.

The Nernst-Planck model presented in the current paper has only two assumed parameters with a clear physical meaning, the energy barrier height, w, and the relative size of the entrance region of the pore, r. This is an important advantage for fitting and interpreting experimental data. The low voltage domain of linear current shown by this model has not yet experimentally explored, despite the fact that it offers direct information about the energy barrier inside a membrane's pore. Additionally, the high voltage domain of linearity validates the values for w and r, when it is accessible. For all these arguments, we believe this simple model for the current-voltage nonlinearity and asymmetry of a pore is a good starting point for explaining the electrical behavior of pores, with wide applications to many research fields.