Cubic-like Features of I–V Relations via Classical Poisson–Nernst–Planck Systems Under Relaxed Electroneutrality Boundary Conditions

Abstract

We focus on higher-order matched asymptotic expansions of a one-dimensional classical Poisson–Nernst–Planck system for ionic flow through membrane channels with two oppositely charged ion species under relaxed electroneutrality boundary conditions. Of particular interest are the current–voltage (I–V) relations, which are used to characterize the two most relevant biological properties of ion channels—permeation and selectivity—experimentally. Our result shows that, up to the second order in $\epsilon =\sqrt{\lambda /r}$, where $\lambda$ is the Debye length and r is the characteristic radius of the channel, the cubic I–V relation has either three distinct real roots or a unique real root with a multiplicity of three, which sensitively depends on the boundary layers because of the relaxation of the electroneutrality boundary conditions. This indicates more rich dynamics of ionic flows under our more realistic setups and provides a better understanding of the mechanism of ionic flows through membrane channels.

1. Introduction

Ion channels are large proteins embedded in membranes that provide a major pathway for cells to communicate with each other and with the outside to transform signals and to conduct group tasks [1,2,3,4]. Ion channel research encompasses two major, interlinked areas: structural characterization and ionic flow analysis. The physical configuration of ion channels is defined by their geometry and the spatial distribution of permanent and polarization charges. Very often, the shape of a typical ion channel has a cylindrical-like domain with a non-uniform cross-sectional area. Within a large class of ion channels, amino acid side chains are distributed mainly over a “short” and “narrow” portion of the channel [5,6,7,8]. The function of channel structures is to select the types of ions and to facilitate the diffusion of ions across cell membranes.

Currently, these permeation and selectivity properties of an ion channel are mainly characterized by the I–V relation that is measured experimentally under different ionic conditions [2,9]. However, all present experimental measurements about ionic flow are of input–output type [2]; that is, the internal dynamics within the channel cannot be measured with the current technology. Mathematical analysis equipped with physically rigorous models becomes essential for a comprehensive exploration of ion channel behavior. Poisson–Nernst–Planck (PNP) systems serve as foundational models for describing ionic motion in these channels. Recently, some nice results in mathematical analysis of the PNP system have been obtained [10,11,12,13,14,15,16,17,18,19,20,21,22].

1.1. Poisson–Nernst– Planck Model for Ionic Flows

Based on the structural characteristics, the basic continuum model for ionic flows is the PNP system that treats the aqueous medium as a dielectric continuum (see [23,24,25,26,27,28,29,30] and the reference therein). Under various conditions, the PNP system can be derived as a reduced model from molecular dynamics [31], from Boltzmann equations [32], and from variational principles [33,34,35]. Considering the key feature of the biological system, the PNP system represents an appropriate model for both analysis and numerical simulations of ionic flows.

In this work, we adopt the one-dimensional steady-state PNP model analyzed in [36] and first proposed by [37], which reads, for $k=1,2,\cdots ,n,$

$$\begin{array}{cc}\hfill & {\displaystyle \frac{{\epsilon }^2}{h\left(x\right)}}{\displaystyle \frac{d}{dx}}\left(h\left(x\right){\displaystyle \frac{d\varphi }{dx}}\right)=-\sum _{j=1}^n{\alpha }_j{c}_j\left(x\right)-Q\left(x\right),{\displaystyle \frac{d{\mathcal{J}}_k}{dx}}=0,h\left(x\right){\displaystyle \frac{d{c}_k}{dx}}+{\alpha }_k{c}_{k}h\left(x\right){\displaystyle \frac{d\varphi }{dx}}=-J_k,\hfill \end{array}$$

(1)

where $x\in [0,1]$ is the coordinate along the axis of the channel (here, the channel is normalized from $x=0$ to $x=1$); $h\left(x\right)$ is the area of the cross-section over the point x; ${\epsilon }^2={\displaystyle \frac{\lambda }{r}},$ where $\lambda$ is the Debye length and r is the characteristic radius of the channel; $\varphi \left(x\right)$ is the electric potential; e is the elementary charge; and for each $k,$ ${\alpha }_k$ is the valence, $c_k$ is the number density, $J_k$ is the flux density (scaled by the diffusion constant), and $Q\left(x\right)$ is the permanent charge density of the channel.

The boundary conditions are, for $i=1,2,...,n$,

$$\varphi \left(0\right)=V,c_k\left(0\right)=L_k>0;\varphi \left(1\right)=0,c_k\left(1\right)=R_k>0.$$

(2)

For given $V,Q\left(x\right),L_k^{\prime }$s, and $R_k^{\prime }$s, if $(\varphi (x;\epsilon ),c_k(x;\epsilon ),J_k\left(\epsilon \right))$ is a solution of the boundary value problem in Equations (1) and (2), then the current–voltage (I–V) relation I is defined by

$$\begin{array}{c}\hfill I=\sum _{s=1}^n{\alpha }_s{J}_s\left(\epsilon \right),\end{array}$$

(3)

which means the dependence of the current on the voltage V for fixed boundary concentrations $L_k$ and $R_k$. Correspondingly, the flow rate of matter $T(V;\epsilon )$ is defined by

$$T(V,\epsilon )=\sum _{s=1}^n{J}_s\left(\epsilon \right).$$

1.2. Electroneutrality Boundary Conditions vs. Boundary Layers

To describe the behavior of channels or functional transistors accurately, it is essential to consider macroscopic reservoirs connected via ion channels [38,39,40,41]. These reservoirs require macroscopic boundary conditions, which introduce boundary layers characterized by concentration and charge distributions. When these boundary layers extend into the device region that operates under atomic-level control, they can significantly impact device behavior. In particular, charge boundary layers may cause long-distance artifacts, as the electric field generated by these layers can extend far beyond the immediate boundary area [12].

Therefore, careful consideration of boundary layer effects is crucial in the study of such systems, especially for ion channel problems. Nonetheless, in many studies focused on the qualitative properties of ionic flows through membrane channels, electroneutrality boundary conditions are imposed at both ends of the channel (see, e.g., [5,6,7,18,36,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66], given by

$$\begin{array}{c}\hfill \sum _{k=1}^n{z}_k{L}_k=\sum _{k=1}^n{z}_k{R}_k=0.\end{array}$$

(4)

Under the condition in Equation (4), the complexity of analyzing the dynamics of ionic flows is significantly reduced. However, this simplification prevents observation of the effects of boundary layers, which emerge from the relaxation of neutrality and contain much richer information about ionic behavior.

To gain a deeper understanding of the mechanisms governing ionic flows through membrane channels, it is crucial to incorporate boundary layer effects into the analysis. Given the sensitivity of electric potentials to these boundary layers, an initial but natural approach is to examine states that are nearly, though not entirely, neutral. This approach, while more challenging, provides a more realistic framework for studying ionic flow dynamics. Building on this concept, the author in [67] examined the cPNP system for one cation and one anion without permanent charges, focusing on the qualitative properties of ionic flows. For this, we introduce parameters $\sigma$ and $\rho$ with $(\sigma ,\rho )\to (1,1)$ as follows:

$$\begin{array}{c}\hfill z_1{L}_1=\sigma (-z_2{L}_2)\mathrm{a}\mathrm{n}\mathrm{d}{z}_1{R}_1=\rho (-z_2{R}_2).\end{array}$$

(5)

Note that $\sigma =1=\rho$ in Equation (5) implies the neutral state. Richer behavior of ionic flows was observed under this new setup. Later, the authors in [20,68] further studied the PNP system under different setups.

All the works show that the boundary layer plays an important role in the analysis of ionic flow properties. In the current work, we study the cubic-like feature of the I–V relations under the assumption of Equation (5) and analyze the boundary layer impacts on it.

2. Mathematical Methods

The boundary value problem in Equations (1) and (2) is treated as a singularly perturbed problem, from which the outer and inner systems for the asymptotic expansions are derived. More precisely, the outer system deals with the internal dynamics of ionic flows within the channel, that is, the characterization of the system Equation (1) for $0<x<1$, while the two inner systems handle the boundary layers occurring at $x=0$ and $x=1$, respectively. The matching principle provides the connection between the boundary layers and the internal dynamics. The outer system and the inner system have different time scales, and one needs to piece them together via asymptotic matching principles [69], from which one is able to derive the I–V relations $I(V;\epsilon )$ and the flow rate of matter $T(V;\epsilon )$ of the following forms:

$$\begin{array}{cc}\hfill I(V;\epsilon )=& I_0\left(V\right)+\epsilon I_1\left(V\right)+{\epsilon }^2{I}_2\left(V\right)+o\left({\epsilon }^2\right),T(V;\epsilon )=T_0\left(V\right)+\epsilon T_1\left(V\right)+{\epsilon }^2{T}_2\left(V\right)+o\left({\epsilon }^2\right).\hfill \end{array}$$

2.1. Previous Results and Assumptions

We recall some results from [36], which are the starting point of our analysis. In [36], the authors obtained $I_0\left(V\right)$ and $I_1\left(V\right)$ for general boundary conditions, and they are given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill T_0=& 2\sqrt{L_1{L}_2}-2\sqrt{R_1{R}_2},I_0={\displaystyle \frac{2(\sqrt{L_1{L}_2}-\sqrt{R_1{R}_2})\left(2V+ln\left(L_1{R}_2\right)-ln\left(L_2{R}_1\right)\right)}{ln\left(L_1{L}_2\right)-ln\left(R_1{R}_2\right)}},\hfill \\ \hfill T_1=& {\displaystyle \frac{4r(I_0+r{T}_0)}{\sqrt{N}(1+r)(1-r)}}-{\displaystyle \frac{4l(I_0+l{T}_0)}{\sqrt{M}(1+l)(1-l)}},\hfill \\ \hfill I_1=& T_0^{-1}{I}_0{T}_1-{\displaystyle \frac{T_0^{-1}{I}_0(a_1{T}_0-a_0{T}_1)}{ln\left(R_1{R}_2\right)-ln\left(L_1{L}_2\right)}}\left({\displaystyle \frac{1}{\sqrt{R_1{R}_2}}}-{\displaystyle \frac{1}{\sqrt{L_1{L}_2}}}\right)-{\displaystyle \frac{4r{T}_0(T_0+r{I}_0)}{\sqrt{R_1{R}_{2}N}(1+r)(1-r)}}\hfill \\ & +{\displaystyle \frac{4l{T}_0(T_0+l{I}_0)}{\sqrt{L_1{L}_{2}M}(1+l)(1-l)}},\hfill \end{array}\end{array}$$

(6)

where $a_0=2\sqrt{L_1{L}_2}=M,N=2\sqrt{R_1{R}_2},a_1=-{\displaystyle \frac{4l(I_0+l{T}_0)}{\sqrt{M}(1+l)(1-l)}},l={\displaystyle \frac{L_2^{\frac{1}{4}}-L_1^{\frac{1}{4}}}{L_2^{\frac{1}{4}}+L_1^{\frac{1}{4}}}},r={\displaystyle \frac{R_2^{\frac{1}{4}}-R_1^{\frac{1}{4}}}{R_2^{\frac{1}{4}}+R_1^{\frac{1}{4}}}}.$

In the current work, we focus on the study of the cubic-like feature of the I–V relations under relaxed boundary concentration conditions. More precisely, we assume, with $z_1=-z_2=1$,

$$\begin{array}{c}\hfill L_1=\sigma L_2=\sigma L,R_1=\rho R_2=\rho R,\end{array}$$

(7)

where both $\sigma$ and $\rho$ are some positive parameters with $(\sigma ,\rho )\to (1,1)$ but not equal to 1 simultaniously.

We now consider $T_0,I_0,T_1$, and $I_1$ to be functions of $(V;\sigma ,\rho )$, expand them at $(\sigma ,\rho )=(1,1)$, and neglect higher-order terms. For convenience, we introduce $f(L,R)={\displaystyle \frac{L-R}{lnL-lnR}}.$ It can be verified directly that $f(L,R)>0$ if $L\ne R$ and $f(L,R)\to R$ as $L\to R$. Careful calculation gives, for $k=0,1,$

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill T_k(\sigma ,\rho )& =T_k(1,1)+T_{k\sigma }(1,1)(\sigma -1)+T_{k\rho }(1,1)(\rho -1),\hfill \\ \hfill I_k(\sigma ,\rho )& =I_k(1,1)+I_{k\sigma }(1,1)(\sigma -1)+I_{k\rho }(1,1)(\rho -1),\hfill \end{array}\end{array}$$

(8)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & T_0(1,1)=2(L-R),T_{0\sigma }(1,1)=L,T_{0\rho }(1,1)=R,I_0(1,1)=2f(L,R)V,\hfill \\ \hfill & I_{0\sigma }(1,1)=f(L,R)-{\displaystyle \frac{L-f(L,R)}{lnL-lnR}}V,I_{0\rho }(1,1)={\displaystyle \frac{f(L,R)-R}{lnL-lnR}}V-f(L,R),T_1(1,1)=0,\hfill \\ \hfill & T_{1\sigma }(1,1)={\displaystyle \frac{f(L,R)V}{\sqrt{2L}}},I_1(1,1)=0,I_{1\sigma }(1,1)={\displaystyle \frac{T_{1\sigma }(1,1)V}{lnL-lnR}}-{\displaystyle \frac{f^3(L,R)V^2-4{(L-R)}^3}{(L-R)L\sqrt{2L}}},\hfill \\ \hfill & T_{1\rho }(1,1)=-{\displaystyle \frac{f(L,R)V}{\sqrt{2R}}},I_{1\rho }(1,1)={\displaystyle \frac{T_{1\rho }(1,1)V}{lnL-lnR}}+{\displaystyle \frac{f^3(L,R)V^2+4{(L-R)}^3}{(L-R)R\sqrt{2R}}}.\hfill \end{array}\end{array}$$

(9)

2.2. Approximation of $I_2(V;\sigma ,\rho )$ as $(\sigma ,\rho )\to (1,1)$

In [36], the authors derived $I_2$ under electroneutrality conditions. To consider the boundary layer effects, we re-derive $I_2$ with $(\sigma ,\rho )\to (1,1)$ (see Appendix A for details). We just present the related result for $T_2(\sigma ,\rho )$ and $I_2(\sigma ,\rho )$. Note that both $T_2$ and $I_2$ also depend on the potential V. For convenience, we did not include it in the expressions. To be specific, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill T_2(\sigma ,\rho )=& T_2(1,1)+T_{2\sigma }(1,1)(\sigma -1)+T_{2\rho }(1,1)(\rho -1),\hfill \\ \hfill I_2(\sigma ,\rho )=& I_2(1,1)+I_{2\sigma }(1,1)(\sigma -1)+I_{2\rho }(1,1)(\rho -1),\hfill \end{array}\end{array}$$

(10)

where

$$\begin{array}{cc}\hfill & T_2(1,1)={\displaystyle \frac{f^2(L,R)(L^2-R^2)V^2}{2L^2{R}^2}},T_{2\sigma }(1,1)={\displaystyle \frac{f^2(L,R)V^2}{2L}}\left({\displaystyle \frac{L+R}{R^2}}-{\displaystyle \frac{f(L,R)(L+R)}{L{R}^2}}+{\displaystyle \frac{1}{L}}\right),\hfill \\ \hfill & T_{2\rho }(1,1)={\displaystyle \frac{f^2(L,R)V^2}{2R}}\left(-{\displaystyle \frac{L+R}{L^2}}+{\displaystyle \frac{f(L,R)(L+R)}{R{L}^2}}-{\displaystyle \frac{1}{R}}\right),\hfill \\ \hfill & I_2(1,1)={\displaystyle \frac{f^2(L,R)}{L^2{R}^2}}\left[{\displaystyle \frac{(L-R)(L^3-R^3)}{3LR}}-f(L,R)\left({\displaystyle \frac{f(L,R)(L^2+LR+R^2)}{3LR}}-{\displaystyle \frac{L+R}{2}}\right)V^2\right]V,\hfill \\ \hfill & I_{2\sigma }(1,1)=V{f}^2(L,R)[{\displaystyle \frac{2f(L,R)(L^3-R^3)V^2}{3L^3{R}^3{(lnL-lnR)}^2}}-{\displaystyle \frac{(L^3-R^3)V^2}{2R^3{L}^2{(lnL-lnR)}^2}}-{\displaystyle \frac{f(L,R)V^2}{2L^3(lnL-lnR)}}\hfill \\ & -{\displaystyle \frac{f(L,R)(L+R)V^2}{R^2{L}^2(lnL-lnR)}}+{\displaystyle \frac{V^2}{2L^2(lnL-lnR)}}+{\displaystyle \frac{(L+R)V^2}{2R^{2}L(lnL-lnR)}}+{\displaystyle \frac{f(L,R)(L+R)V^2}{4R^2{L}^2}}\hfill \\ & -{\displaystyle \frac{f(L,R)(L^3-R^3)}{3R^3{L}^3}}+{\displaystyle \frac{2(L^3-R^3)}{3R^3{L}^2}}+{\displaystyle \frac{L-R}{2L^3}}],\hfill \end{array}$$

$$\begin{array}{cc}\hfill I_{2\rho }(1,1)=& V{f}^2(L,R)[-{\displaystyle \frac{2f(L,R)(L^3-R^3)V^2}{3L^3{R}^3{(lnL-lnR)}^2}}+{\displaystyle \frac{(L^3-R^3)V^2}{2R^2{L}^3{(lnL-lnR)}^2}}+{\displaystyle \frac{f(L,R)V^2}{2R^3(lnL-lnR)}}\hfill \\ & +{\displaystyle \frac{f(L,R)(L+R)V^2}{R^2{L}^2(lnL-lnR)}}-{\displaystyle \frac{V^2}{2R^2(lnL-lnR)}}-{\displaystyle \frac{(L+R)V^2}{2R{L}^2(lnL-lnR)}}\hfill \\ & -{\displaystyle \frac{f(L,R)(L+R)V^2}{4R^2{L}^2}}+{\displaystyle \frac{f(L,R)(L^3-R^3)}{3R^3{L}^3}}-{\displaystyle \frac{2(L^3-R^3)}{3R^3{L}^2}}-{\displaystyle \frac{L-R}{2R^3}}].\hfill \end{array}$$

For convenience in our following discussion, we rewrite $I_2(V;\sigma ,\rho )$ as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & I_2(V;\sigma ,\rho )\hfill \\ \hfill & ={\displaystyle \frac{f^2(L,R)}{L^2{R}^2}}[{\displaystyle \frac{(L+R)f(L,R)}{2}}-{\displaystyle \frac{(L^3-R^3)f(L,R)}{3RL(lnL-lnR)}}+{\displaystyle \frac{2(L^3-R^3)f(L,R)}{3RL{(lnL-lnR)}^2}}(\sigma -\rho )\hfill \\ \hfill & -{\displaystyle \frac{(L^3-R^3)\left(L(\sigma -1)-R(\rho -1)\right)}{2RL{(lnL-lnR)}^2}}-{\displaystyle \frac{f(L,R)\left(R^3(\sigma -1)-L^3(\rho -1)\right)}{2RL(lnL-lnR)}}\hfill \\ & +f(L,R)(L+R)\left({\displaystyle \frac{1}{4}}-{\displaystyle \frac{1}{lnL-lnR}}\right)(\sigma -\rho )+{\displaystyle \frac{(L^2+LR+R^2)(\sigma -\rho )}{2(lnL-lnR)}}]V^3\hfill \\ \hfill & +{\displaystyle \frac{f^2(L,R)(L^3-R^3)}{L^3{R}^3}}[{\displaystyle \frac{L-R}{3}}-{\displaystyle \frac{f(L,R)}{3}}(\sigma -\rho )+{\displaystyle \frac{2\left(L(\sigma -1)-R(\rho -1)\right)}{3}}\hfill \\ & +{\displaystyle \frac{R^3(\sigma -1)-L^3(\rho -1)}{2(L^2+LR+R^2)}}]V.\hfill \end{array}\end{array}$$

(11)

3. Results

Under the assumption of Equation (7), for the I–V relations, up to order o$\left({ϵ}^2\right)$, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill I(V;\sigma ,\rho )& =I_0(V;\sigma ,\rho )+ϵI_1(V;\sigma ,\rho )+{ϵ}^2{I}_2(V;\sigma ,\rho )\hfill \\ \hfill & ={ϵ}^2{f}^2(L,R)M_1{V}^3+ϵM_2{V}^2+({ϵ}^2{M}_3+M_4)V+M_5+ϵM_6,\hfill \end{array}\end{array}$$

(12)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & M_1={\displaystyle \frac{f(L,R)(L+R)}{2R^2{L}^2}}-{\displaystyle \frac{f(L,R)(L^3-R^3)}{3R^3{L}^3(lnL-lnR)}}+[{\displaystyle \frac{2f(L,R)(L^3-R^3)}{3R^3{L}^3{(lnL-lnR)}^2}}-{\displaystyle \frac{f(L,R)}{2L^3(lnL-lnR)}}\hfill \\ & -{\displaystyle \frac{L^3-R^3}{2R^3{L}^2{(lnL-lnR)}^2}}-{\displaystyle \frac{f(L,R)(L+R)}{R^2{L}^2(lnL-lnR)}}+{\displaystyle \frac{1}{2L^2(lnL-lnR)}}+{\displaystyle \frac{L+R}{2R^{2}L(lnL-lnR)}}\hfill \\ & +{\displaystyle \frac{f(L,R)(L+R)}{4R^2{L}^2}}]\sigma +[-{\displaystyle \frac{2f(L,R)(L^3-R^3)}{3R^3{L}^3{(lnL-lnR)}^2}}+{\displaystyle \frac{L^3-R^3}{2R^2{L}^3{(lnL-lnR)}^2}}\hfill \\ & +{\displaystyle \frac{f(L,R)}{2R^3(lnL-lnR)}}+{\displaystyle \frac{f(L,R)(L+R)}{R^2{L}^2(lnL-lnR)}}-{\displaystyle \frac{L^2-(L+R)R}{2R^2(lnL-lnR)}}-{\displaystyle \frac{f(L,R)(L+R)}{4R^2{L}^2}}]\rho ,\hfill \\ \hfill & M_2={\displaystyle \frac{f^2(L,R)(L-R)(L^3-R^3)}{6R^3{L}^3}}+f^2(L,R)\left[{\displaystyle \frac{(2L-f(L,R))(L^3-R^3)}{3R^3{L}^3}}+{\displaystyle \frac{L-R}{2L^3}}\right]\sigma \hfill \\ & +f^2(L,R)\left[{\displaystyle \frac{(f(L,R)-2R)(L^3-R^3)}{3R^3{L}^3}}-{\displaystyle \frac{L-R}{2R^3}}\right]\rho ,\hfill \\ \hfill & M_3={\displaystyle \frac{f(L,R)}{lnL-lnR}}\left({\displaystyle \frac{\sigma -1}{\sqrt{2L}}}-{\displaystyle \frac{\rho -1}{\sqrt{2R}}}\right)-{\displaystyle \frac{f^2(L,R)}{lnL-lnR}}\left({\displaystyle \frac{\sigma -1}{L\sqrt{2L}}}-{\displaystyle \frac{\rho -1}{R\sqrt{2R}}}\right),\hfill \\ \hfill & M_4={\displaystyle \frac{L(\sigma +1)-R(\rho +1)}{lnL-lnR}}-{\displaystyle \frac{f(L,R)(\sigma -\rho )}{lnL-lnR}},\hfill \\ \hfill & M_5=-4{(L-R)}^2\left({\displaystyle \frac{\sigma -1}{L\sqrt{2L}}}-{\displaystyle \frac{\rho -1}{R\sqrt{2R}}}\right),M_6=f(L,R)(\sigma -\rho ).\hfill \end{array}\end{array}$$

(13)

We are ready to state our main result.

Theorem 1. 

Under the assumption of Equation (7) and $L\ne R$, up to the second-order in ϵ, as $(\sigma ,\rho )\to (1,1)$, for the I–V relation $I=I\left(V\right)$ defined in (12), one has the following:

(i)

$I\left(V\right)=0$ has three distinct real roots if one of the following conditions holds:

(i1)

$\sigma =\rho \ne 1$;

(i2)

$\sigma =1$ and $\rho <1$;

(i3)

$\sigma >1$ and $\rho =1$;

(i4)

$\rho <\sigma <1$;

(i5)

$\rho <1<\sigma$;

(i6)

$\rho <\sigma <1$.

(ii)

$I\left(V\right)=0$ has a unique real root with multiplicity 3 if one the the following conditions holds:

(ii1)

$\sigma =1$ and $\rho >1$;

(ii2)

$\sigma <1$ and $\rho =1$;

(ii3)

$\sigma <\rho <1$;

(ii4)

$\sigma <1<\rho$;

(ii5)

$1<\sigma <\rho$.

Proof. 

The proofs for different cases listed in the theorem are similar. Here, we will just provide a detailed discussion of case (i1).

Direct calculation gives

$$I^{\prime }\left(V\right)=3{ϵ}^2{f}^2(L,R)M_1{V}^2+2ϵM_{2}V+{ϵ}^2{M}_3+M_4,$$

and a quadratic function in the potential V with the discriminant $\Delta$ given by

$$\Delta =4{ϵ}^2(M_2^2-3f^2(L,R)M_1{M}_4-3{ϵ}^2{f}^2(L,R)M_1{M}_3)=4{ϵ}^2\left(M_2^2-3f^2(L,R)M_1{M}_4+o\left({ϵ}^2\right)\right).$$

With $ϵ>0$ small, the sign of $\Delta$ is dominated by the term $M_2^2-3f^2(L,R)M_1{M}_4$.

For the case with $\sigma =\rho \ne 1$, from Equation (13), one has

$$\begin{array}{cc}\hfill M_1=& {\displaystyle \frac{f(L,R)(L+R)}{2L^2{R}^2}}-{\displaystyle \frac{f(L,R)(L^3-R^3)}{3L^3{R}^3(lnL-lnR)}},M_2={\displaystyle \frac{f^2(L,R)(L-R)(L^3-R^3)}{6L^3{R}^3}}(\sigma +1),\hfill \\ \hfill M_4=& f(L,R)(\sigma +1).\hfill \end{array}$$

It follows that

$$\begin{array}{cc}\hfill M_2^2-3f^2(L,R)M_1{M}_4=& M_2^2-{\displaystyle \frac{f^4(L,R)(\sigma +1)}{2R^3{L}^3(lnL-lnR)}}Z(L,R),\hfill \end{array}$$

where $Z(L,R)=3LR(L+R)(lnL-lnR)-2(L^3-R^3).$ We now claim that $Z(L,R)<0$ for $L>R$. To start, we rewrite $Z(L,R)$ as

$$Z(L,R)=R^3{Z}_0(L,R),$$

where

$$Z_0(L,R)=3{\displaystyle \frac{L}{R}}\left({\displaystyle \frac{L}{R}}+1\right)ln{\displaystyle \frac{L}{R}}-2\left({\displaystyle \frac{L^3}{R^3}}-1\right).$$

Upon introducing $x={\displaystyle \frac{L}{R}}$>1, one has

$$Z_0\left(x\right)=3(x^2+x)lnx-2(x^3-1).$$

Note that $Z_0\left(1\right)=Z_0^{\prime }\left(1\right)=Z_0^{″}\left(1\right)=0$ and

$$Z_0^{‴}\left(x\right)=-{\displaystyle \frac{12}{x^2}}\left({\left(x-{\displaystyle \frac{1}{4}}\right)}^2+{\displaystyle \frac{3}{16}}\right)<0\mathrm{f}\mathrm{o}\mathrm{r} \mathrm{a}\mathrm{l}\mathrm{l}x>1.$$

Therefore, $Z_0\left(x\right)<0$ for all $x>1$, and hence, $Z(L,R)<0$ for $L>R$. It follows that $-3f^2{M}_1{M}_4>0$. It follows that $M_2^2-3f^2(L,R)M_1{M}_4>0$. Therefore, $I^{\prime }\left(V\right)=0$ has two distinct real roots, which are given by

$$V_1={\displaystyle \frac{-2ϵM_2+\sqrt{\Delta }}{6{ϵ}^2{f}^2(L,R)M_1}}\mathrm{a}\mathrm{n}\mathrm{d}{V}_2={\displaystyle \frac{-2ϵM_2-\sqrt{\Delta }}{6{ϵ}^2{f}^2(L,R)M_1}}.$$

Direct calculation gives

$$I\left(V_1\right)I\left(V_2\right)={\displaystyle \frac{M_4^2}{27{ϵ}^2{f}^4(L,R)M_1^2}}(7f^2(L,R)M_1{M}_4-M_2^2+o\left(ϵ\right)),$$

where

$$\begin{array}{cc}\hfill 7f^2(L,R)M_1{M}_4-M_2^2=& {\displaystyle \frac{7f^4(L,R)(\sigma +1)}{6R^3{L}^3(lnL-lnR)}}Z(L,R)-M_2^2<0.\hfill \end{array}$$

This implies that $I\left(V\right)=0$ has three distinct real roots. We complete the proof. □

4. Discussion

Electrodiffusion—the diffusion and migration of charged particles—plays a critical role in the understanding of nature and in the inventions of modern electronic devices. Ionic flows through ion channels comprise an important particular process of electrodiffusion that shows extremely rich dynamics depending on complicated nonlinear interactions of many physical parameters, such as channel structures, boundary concentrations, electric potential differences, diffusion coefficients, ion sizes, etc. This makes the study of flow properties of interest extremely challenging. Mathematical analysis equipped with physically rigorous models becomes essential for a comprehensive exploration of ion channel behavior.

Currently, the I–V relation is the main tool used experimentally to characterize the two most relevant properties of an ion channel: permeation and selectivity. However, current experimental methods provide only external input–output data and fall short of revealing the detailed internal ionic dynamics within channels. This gap presents a significant challenge to fully understand the properties of ion channels. What we hope for is to rigorously analyze a mesoscopic model of an ion channel in sufficient detail to elucidate its salient features and their impact on conduction and selection. In this work, we employ both an asymptotic matching principle and regular perturbation analysis to a 1D PNP model for ionic flows through membrane channels. Of particular interest is the cubic-like feature of the I–V relations and its sensitive dependence on the boundary concentrations under more realistic setups. To be specific, under relaxed electroneutrality boundary conditions in Equation (7), up to the second order in $ϵ,$ the I–V relations still show the cubic-like feature for some setups as stated in Theorem 1. Meanwhile, the ionic flow properties are sensitive to the boundary layers because of the relaxation of electroneutrality boundary conditions. Under conditions (ii1)–(ii5) stated in Theorem 1, the three distinct real roots of $I(V;\sigma ,\rho )=0$ degenerate to a unique real root of multiplicity 3, and the cubic-like feature of the I–V relations vanishes.

The work in this paper is our first step in trying to better understand the mechanism of ionic flows through membrane channels, particularly the sensitive dependence of the I–V relations on the boundary concentration conditions, which show some degenerate phenomena (three distinct real roots degenerate into a unique real root with multiplicity three) as stated in (ii) of the Theorem 1. This leads to richer dynamics of ionic flows through membrane channels. The study in the current work provides efficient ways for one to adjust/control boundary conditions to observe distinct phenomena of ionic flows both numerically and experimentally. The cubic-like feature characterized in this work is consistent with the one adopted in the FitzHuge–Nagumo simplification of the famous Hodgkin–Huxley systems, which describe the propagation of action potential of an ensemble of channels in a biological membrane. This is a reasonable and necessary step before we consider more complicated PNP models in the future, such as the one with nonzero permanent charges and ion-specific finite ion sizes. Finally, we point out that without permanent charges considered in the model, the channel geometry will not affect our discussion of the ionic flow properties (see Remark 2.2 in [43] for more detailed discussion).