Mathematical Analysis on Current–Voltage Relations via Classical Poisson–Nernst–Planck Systems with Nonzero Permanent Charges under Relaxed Electroneutrality Boundary Conditions

We focus on a quasi-one-dimensional Poisson–Nernst–Planck model with small permanent charges for ionic flows of two oppositely charged ion species through an ion channel. Of particular interest is to examine the dynamics of ionic flows in terms of I–V (current–voltage) relations with boundary layers due to the relaxation of neutral conditions on boundary concentrations. This is achieved by employing the regular perturbation analysis on the solutions established through geometric singular perturbation analysis. Rich dynamics are observed, particularly, the nonlinear interplays among different physical parameters are characterized. Critical potentials are identified, which play critical roles in the study of ionic flows and can be estimated experimentally. Numerical simulations are performed to further illustrate and provide more intuitive understandings of our analytical results.


Introduction
Ion channels are pore-forming membrane proteins allowing charged particles to pass through the channel pore. Ion channels are embedded in cell membranes, which provide a major medium for cells to communicate with each other and with outside environment ( [1][2][3]). In this way, ion channels control a wide range of biological functions, in particular, many varied functions are necessary for life (see [2] for more discussion). Clinically, malfunctioning channels cause cystic fibrosis, cholera, neuronal disorders, and many other diseases ( [4]). Therefore, it is significant to explore the mechanism of ion channels.
The study of ion channels generally consists of two related major topics: structures of ion channels and ionic flow properties. The physical structure of ion channels is defined by the channel shape and the spacial distribution of permanent charges and the polarity of these charges. For open channels with given structures, the main interest is in the study of its electrodiffusion property. The most challenging part in examining the properties of ionic flows through membrane channels is the characterization of the nonlinear interplays among specific physical parameters involved in the system, particularly, the boundary concentrations and membrane potential, permanent charge distribution within the channel, channel geometry and diffusion coefficients. On the other hand, all present experimental measurements about ionic flow are of input-output type ( [1]); that is, the internal dynamics within the channel cannot be measured with the current technology. Therefore, it is extremely difficult to extract coherent properties or to formulate specific characteristic quantities from the experimental measurements.
Mathematical analysis plays important and unique roles for generalizing and understanding the principles that allow control of electrodiffusion, explaining the mechanics of observed biological phenomena and for discovering new ones, under the assumption that a more or less explicit solution of the associated mathematical model can be obtained. Recently, there have been some successes in mathematical analysis of Poisson-Nernst-Planck (PNP) models for ionic flows through membrane channels (  etc.). Particularly, for those that were studied under the dynamical system framework of geometric singular perturbation analysis, interesting phenomena of ionic flows were observed for relatively simple setups. However, for all the works mentioned here, the so-called electroneutrality boundary conditions are reinforced, and the boundary layers, which play critical roles in the study of ion channel problems, disappear.
To further understand the correlations/interactions among ions, we consider a PNP model with nonzero but small permanent charges under relaxed electroneutrality boundary concentration conditions. Of particular interest is the effect on the I-V relations from boundary layers, which can be mathematically extracted from solutions of the PNP system. Our study will take great advantage of the work conducted in [11]. To be specific, in [11], the authors treat the nonzero permanent charge as a small parameter, and employ the regular perturbation analysis to the solutions, from which they obtain the zerothorder and first-order individual fluxes under electroneutrality boundary conditions. This will be our starting point. More precisely, the individual fluxes obtained in [11] allow us to define the current-voltage (I-V) relations directly (see (15) in Section 2.2). To examine the boundary layer effects on the I-V relations, we relax the electroneutrality boundary conditions by introducing two positive parameters, σ and ρ close to but not simultaneously equal to 1 (see our discussion in Section 1.3). We next employ the regular perturbation analysis to these two parameters, in other words, we expand our I-V relations at (σ, ρ) = (1, 1) up to the first order and neglect higher order terms. Our main interest is then in the first-order terms, the leading terms that contain boundary layer effects.
Based on the fact that the channel is narrow and one can effectively view it as a onedimensional channel [0, l], where l (with unit nm) is the length of the channel together with the baths that the channel links. A quasi-one-dimensional steady-state PNP model for a mixtures of n ion species though a single channel reads (first proposed in [41]) where x ∈ [0, 1] is the coordinate along the axis of the channel that is normalized to [0, 1], h(x) is the area of cross-section of the channel over the location x.
For system (1), we have the following boundary conditions (see [8] for a reasoning), for k = 1, 2, . . . , n, where • e ≈ 1.60 × 10 −19 (C = coulomb) is the elementary charge; T is the absolute temperature (unit K (kelvin)), it is T = 273.16 (K); • Φ(x) is the electric potential with the unit V = Volt = JC −1 ; • Q(x) is the permanent charge density of the channel (with unit 1/m 3 ); h(x) represents the area of the cross-section over the point x (with unit m 2 ); • n is the number of distinct types of ion species (with unit 1); • for the jth ion species; c j is the number density (with unit 1/m 3 ); -z j is the valence (the number of charges per particle with unit 1); µ j is the electrochemical potential (with unit J = CV); -J j is the number flux density (with unit 1/s) -the number of particles across each cross-section per unit time; -D j (x) is the diffusion coefficient (with unit m 2 /s).

Permanent Charges
It is known that the spatial distribution of side chains in a specific channel defines the permanent charge of the channel. Although some information could be obtained without considering the permanent charge and by focusing on the effects of other system parameters, such as boundary conditions, ion valences, ion sizes, etc., we believe that different channel types differ mainly in the distribution of permanent charge ( [2]). To better understand the dependence of ionic flows on permanent charges, we demonstrate that the role of permanent charges in membrane channels is similar to the role of doping profiles in semiconductor devices. Semiconductor devices are similar to membrane channels in the way that they both use atomic-scale structures to control macroscopic flows from one reservoir to another. Ions move a lot like quasi-particles move in semiconductors. Roughly, holes and electrons are the cations and anions of semiconductors. Semiconductor technology depends on the control of migration and diffusion of quasi-particles of charge in transistors and integrated circuits. Doping is the process of adding impurities into intrinsic semiconductors to modulate its electrical, optical, and structural properties ( [42,43]). Roughly speaking, one may understand in the following sense, doping provides the charges that acid and basic side chains provide in a protein channel. There is no doubt that, for both ion channels and semiconductors, permanent charges add an additional component-probably the most important one-to their rich behavior. In general, the permanent charge Q(x) is modeled by a piecewise constant function, that is, we assume, for a partition x m ] are viewed as the reservoirs where there is no permanent charge).
In [8], under the framework of geometric singular perturbation theory, the existence and uniqueness (local) was established for the boundary value problem (1) and (2) with one cation and one anion and the permanent charge function modeled by where Q 0 is some nonzero constant. Due to the challenge in obtaining explicit expressions of the I-V relation with nonzero permanent charges, in [11], the author studied the case with Q 0 in (3) being small and employed regular perturbation analysis (viewing Q 0 as a small perturbation to the solutions of the system (1) and (2)) to further study the effects on ionic flows from the permanent charges. The analysis in [11] (Proposition 4.11 and its following discussion) indicates that to optimize the effect of the permanent charge, a short and narrow neck within which the permanent charge is confined, is expected. This indicates the critical role that the permanent charge plays in the study of ionic flow properties of interest.

Relaxed Electroneutrality Boundary Conditions
To describe the actual behavior of channels or useful transistors, macroscopic reservoirs linked by ion channels must be included ( [30,31,44,45]). Macroscopic boundary conditions that describe such reservoirs introduce boundary layers of concentration and charge. If those boundary layers reach into the part of the device that performs atomic control, they prominently influence its behavior. Particularly, boundary layers of charge are probably to produce artifacts over long distances because the electric field spreads a long way.
The boundary layer should be handled more carefully in the study of such problems, particularly for ion channel problems. In [8,13,26], the boundary layer is characterized partially, mainly in establishing the existence and local uniqueness result of the PNP system. However, the effects from the boundary layers on ionic flows, which in general carry more rich information, are not analyzed. This is because, very often, when examine the qualitative properties of ionic flows through ion channels, electroneutrality boundary conditions are naturally enforced at both ends of the channel (see, e.g., [5,6,11,12,20,21,[46][47][48][49][50][51]). They are defined as However, under the condition (4), the two boundary layers disappear (see [52] for more detailed discussion).
To better understand the mechanism of ionic flows through membrane channels, the boundary layer effects should be carefully considered during the study. Meanwhile, due to the sensitivity of electric potentials on boundary layers, a first but natural step is to study the state that is not neutral but close to. More precisely, one may assume (taking n = 2 in (4), for example) where σ and ρ are some positive constants close to but not equal to 1 simultaneously (σ = 1 = ρ in (5) implies neutral state). Following this idea, some recent works (see [22,52,53] for examples) have shown that more rich qualitative properties of ionic flows were observed while boundary layers are involved. Particularly, the authors in [52] analyzed the boundary layer effects on individual fluxes via PNP system with nonzero but small permanent charges. All the works indicate the importance of the role played by the boundary layer in the study of ionic flow properties of interest.

Problem Set-Up
In this work, we take the same setting as that in [11] but without assuming electroneutrality boundary conditions (4) for n = 2, which includes: (A1).We consider two charged particles (n = 2) with z 1 > 0 and z 2 < 0; (A2).The PNP model only includes the ideal component µ id i (X) of the electrochemical potential defined by where c 0 is some characteristic number density. (A3).ε r (X) = ε r and D i (X) = D i .
We will assume (A1)-(A3) from now on. We first make the following dimensionless rescaling ( [11]). Let Correspondingly, the boundary value problem (1) and (2) becomes with the boundary conditions

Methods
We take the advantage of the work done in [11], and further employ regular perturbation analysis on the parameters σ and ρ introduced in (5).

Previous Results
To get started, we first recall some results from [11], which are fundamental for our following discussion. Treating the nonzero permanent charge |Q 0 | small compared to the boundary concentrations L k 's and R k 's, the authors in [11] where J k = D k J k . It follows that J k0 = D k J k0 and J k1 = D k J k1 , where with Here, where, with H(x) = x 0 1 h(s) ds, We define the following function, which will be used often in our analysis. For t > 0, set One establishes easily that

Main Interest and Regular Perturbation Analysis
The most basic function of membrane channels is to regulate the permeability of membranes for a given species of ions and to select the types of ions and to facilitate and modulate the diffusion of ions across cell membranes. Currently, the permeation and selectivity properties of ion channels are usually characterized by the I-V relations measured experimentally [27,54]. Our main interest is to analyze the qualitative properties of the I-V relations and characterize the nonlinear interplays between physical parameters. More precisely, we consider where For convenience, we define five functions P 00 = P 00 (σ, ρ), P 01 = P 01 (σ, ρ), P 10 = P 10 (σ, ρ), P 11 = P 11 (σ, ρ) and P 12 = P 12 (σ, ρ) by , Then, I 0 and I 1 can be rewritten as Recall that our main interest in this work is to examine the qualitative properties of the I-V relations close to the state of electroneutrality, more precisely, based on our set-ups, it is the case as (σ, ρ) → (1, 1). We now employ the regular perturbation analysis and expand I k (V; σ, ρ) for k = 0, 1 at (σ * , ρ * ) = (1, 1) up to the first order and neglect higher orders. Careful calculations give where Here where ω(α; 1, 1) = (1 − α)L 1 + αR 1 and ω(β; 1, 1) = (1 − β)L 1 + βR 1 .

Results
In this section, based on the regular perturbation analysis provided in Section 2.2, we further examine the qualitative properties of ionic flows with boundary layers. To provide more intuitive illustrations of the analytical results obtained in Sections 3.1 and 3.2, numerical simulations under different set-ups are performed in Section 3.3.

Analysis on I 0d
Note that I 0d is a linear function in the potential V. The sign of the slope is of particular interest for us, and the following result can be established.
Proof. The proof is straightforward, and we omit it here.
This completes the proof.
We now turn to the sign of I 1d (V). We first consider the sign of g 4 (α), which will be used in the proof of Theorem 3. Here, g 4 (α) is defined by For the function g 4 (α), the following result can be established.
Proof. The discussion is then straightforward, and we omit it here.

Remark 1.
The function g 4 (α) defined above is very complicated to analyze. The number of zeros depends on the parameter t sensitively. In this work, we just consider the simplest case.
Proof. It is clear that ∆ has the same sign as that of ∆ 1 . For 1 < t < 2, one has From Lemma 4, we can easily obtain the sign of lim β→α ∆ 1 , and hence know the sign of ∆ 1 .

Remark 2.
In Theorem 3, the result established further depends on the order of some critical values related to channel geometry, α * , α 1 , α 3 and α 4 . We demonstrate that the order is not unique, which further depends on the nonlinear interaction between other system parameters, particularly, a, b, the jumping points of the permanent charge, and L 1 , R 1 , the boundary concentrations. However, similar result should be obtained. Moreover, the result stated in Theorem 3 indicates that under different conditions on the channel geometry, the small positive permanent charge Q 0 can either enhance the current I 0d (V; σ, ρ) or reduce it. Taking the statement (i) as an example, for α ∈ (0, α * ) with V d 11 < V d 12 , one has the small permanent charge Q 0 enhances I 0d . Furthermore, the result stated in Theorem 2 provides an efficient way to adjust the effects from the leading term I 1d that contains small permanent charge.

Numerical Simulations
In this part, numerical simulations are performed to provide more intuitive illustrations of some analytical results. To be specific, we numerically identify the critical potentials V d 01 , V d 11 , V d 12 and V c 1 , which characterized the effects caused by the appearance of boundary layers. To further illustrate the boundary layer effects on ionic flows, we also numerically obtain the zeroth-order (respectively, the first-order) I-V relations in small positive Q 0 without and with boundary layers, respectively. Corresponding critical potentials for each setup are identified, from which one is able to observe the effects on ionic flows from boundary layers clearly.
To get started, we rewrite the system (7) and (8) as a system of first-order ordinary differential equations. Upon introducing u = ε d dx φ, one has with boundary conditions In our simulations to system (23) and (24), we take z 1 = −z 2 = 1, L 1 = 12, R 1 = 8, ε = 0.01, Q 0 = 0.008, a = 0.4, b = 0.48, D 1 = 2.032, D 2 = 1.334, We comment that the choice of h(x) is based on the fact that the ion channel is cylinderlike, and the variable cross-section area is chosen to reflect the fact that the channel is not uniform and much narrower in the neck (0.4 < x < 0.48) than other regions ( [11]). We further take r 0 = 0.5 and the function h(x) is then continuous at the jumping points a = 0.4 and b = 0.48. Different models for h(x) may be chosen, and similar numerical results should be obtained.
Our numerical simulations show that: (A) The term I 0d (V) is increasing in the potential V, and there exists a unique zero V d . This is consistent with the first statement in Theorem 1 (see the left figure in Figure 1); (B) The term I 1d (V), as a quadratic function, under our setup, it is concave up with two zerosV d 11 and V d 12 . At V c 1 , the critical point, I 1d (V) attains its global minimum.
The numerical results are consistent with the first statement in Theorem 2 and Theorem 3, respectively (see the right figure in Figure 1).
To further illustrate the boundary layer effects on the I-V relations, we also performed numerical simulations to the system with electroneutrality boundary concentra-tion conditions (no boundary layers, the left figures in both Figures 2 and 3) and relaxed electroneutrality boundary conditions (appearance of boundary layers, the right figures in both Figures 2 and 3), respectively. In each setup, the corresponding critical potentials are identified. To be specific, V EN 01 , V EN 11 and V EN 12 are the critical potentials detected under electroneutrality boundary conditions, while V bd 01 , V bd 11 and V bd 12 are the ones with boundary layers. From the numerical simulations, one observes V EN 01 < V bd 01 , V EN 11 > V bd 11 and V EN 12 > V bd 12 . The observation indicates the important role played by the boundary layers. We take I 0 (V) as an example. For many studies on PNP type models, electroneutrality boundary conditions are applied, based on our numerical result, this means I 0 (V; 1, 1) > 0 for V > V EN 01 . However, if the boundary layers are considered, I 0 (V; σ, ρ) < 0 for V EN 01 < V < V bd 01 , even for (σ, ρ) → (1, 1). The dynamics of ionic flows are totally different.

Concluding Remarks
We study a one-dimensional steady-state Poisson-Nernst-Planck system with two oppositely charged ion species and small but nonzero permanent charges. The main purpose is to understand the problem from the mathematical point of view, which should provide some insights into related studies of ion channel problems. Particularly, we focus on the qualitative properties of the I-V relations with boundary layers due to the violation of the electroneutrality boundary conditions, which are studied from two directions: (i) Boundary layer effects on the zeroth-order (in Q 0 ) I-V relations in terms of the function I 0d (V; σ, ρ); (ii) Boundary layer effects on the first-order (in Q 0 ) I-V relations in terms of the function I 1d (V; σ, ρ).
Detailed analysis along each direction is provided, which includes the signs of I 0d (V; σ, ρ) and I 1d (V; σ, ρ) and their monotonicity. From the study, one can better understand the mechanism of ionic flows through membrane channels, particularly the internal dynamics of ionic flows, which are non-intuitive and cannot be detected by current technology. Critical potentials that balance the small permanent charge effects on the I-V relations are identified, and their critical roles played in the study of ionic flow properties are characterized. Numerical simulations are performed to provide more intuitive illustrations of the analytical results, and they are consistent. Among others, we find (I) The monotonicity of I 0d depends sensitively on the parameter t defined by L 1 /R 1 and the boundary layers through the parameters σ and ρ; (II) The sign of I 1d depends sensitively on the interplays among system parameters, particularly, the parameter t defined by the ratio L 1 /R 1 , and the parameter α = H(a)/H(1) representing the channel geometry, while the monotonicity of I 1d is only sensitive on α.
Finally, we comment that the setup in this work is relatively simple, and the study in this work is the first step in the analysis of more realistic models. The simple model considered allows us to obtain a more explicit expression of the I-V relations in terms of physical parameters of the problem so that we are able to extract concrete information of the effects from boundary layers and small but nonzero permanent charges. Moreover, the analysis in this simpler setting provides further understanding of the qualitative properties of the I-V relations through membrane channels, and detailed characterization of the nonlinear interplay among different system parameters. The critical potentials identified in this work are critical for one to observe different ionic flow properties through membrane channels. More importantly, some of them can be approximated experimentally.