Competition between Cations via Classical Poisson–Nernst–Planck Models with Nonzero but Small Permanent Charges

We study a one-dimensional Poisson–Nernst–Planck system for ionic flow through a membrane channel. Nonzero but small permanent charge, the major structural quantity of an ion channel, is included in the model. Two cations with the same valences and one anion are included in the model, which provides more rich and complicated correlations/interactions between ions. The cross-section area of the channel is included in the system, and it provides certain information of the geometry of the three-dimensional channel, which is critical for our analysis. Geometric singular perturbation analysis is employed to establish the existence and local uniqueness of solutions to the system for small permanent charges. Treating the permanent charge as a small parameter, through regular perturbation analysis, we are able to derive approximations of the individual fluxes explicitly, and this allows us to study the competition between two cations, which is related to the selectivity phenomena of ion channels. Numerical simulations are performed to provide a more intuitive illustration of our analytical results, and they are consistent.


Introduction
Ion channels are large proteins embedded in cell membranes with a hole down their middle that provides a controllable path for electrodiffusion of ions (mainly Na + , K + , Ca ++ and Cl − ) through biological membranes, establishing communications among cells and the external environment [1][2][3]. In general, the study of ion channels 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 and polarization charge. Very often, the shape of a typical ion channel is approximated by a cylindrical-like domain with a non-uniform cross-section area. Within a large class of ion channels, amino acid side chains are distributed mainly over a relatively "short" and "narrow" portion of the channel, where acidic side chains contribute permanent negative charges and basic side chains contributes permanent positive charges, and this is analogous to the doping of semiconductor devices, e.g., bipolar PNP and NPN transistors [1,4].
With a given structure of an open channel, the main interest is to understand its electrodiffusion property. Mathematical analysis plays important and unique roles for generalizing and understanding the principles that allow control of electrodiffusion, explaining mechanics of observed biological phenomena and for discovering new ones, assuming a more or less explicit solution of the associated mathematical model can be obtained. However, in general, the latter is too much to expect. Recently, there have been some successes in mathematical analysis of Poisson-Nernst-Planck (PNP) models for ionic flows through membrane channels [5][6][7][8][9][10][11][12][13][14].
It is known that more sophisticated models [24][25][26][27][28][29][30][31] have also been developed which can model the physical problem more accurately, however, it is very challenging to examine their dynamics analytically and even computationally. Considering the key feature of the biological system, the PNP system represents an appropriate model for both analysis and numerical simulations of ionic flows. The simplest PNP system is the classical Poisson-Nernst-Planck (cPNP) system that includes the ideal component µ id k (X) in (4) only. The ideal component µ id k contains contributions by considering ion particles as point charges and ignoring the ion-to-ion interaction. It has been shown by some numerical studies that classical PNP models provide good qualitative agreement with experimental data for I-V relations [20,32]. The classical PNP models have been simulated and analyzed extensively (see, e.g., [7,[12][13][14]20,[32][33][34][35][36][37][38][39]).
For ionic solutions with n ion species, the PNP system reads ∇ · ε r (r)ε 0 ∇Φ = −e n ∑ s=1 z s C s + Q(r) , where r ∈ Ω with Ω being a three-dimensional cylindrical-like domain representing the channel, Q(r) is the permanent charge density, ε r (r) is the relative dielectric coefficient, ε 0 is the vacuum permittivity, e is the elementary charge, k B is the Boltzmann constant, T is the absolute temperature; Φ is the electric potential. Furthermore, for the kth ion species, C k is the concentration, z k is the valence, µ k is the electrochemical potential depending on Φ and {C j }, J k is the flux density, and D k (r) is the diffusion coefficient. Based on the fact that ion channels have narrow cross-sections relative to their lengths, reduction of the three-dimensional steady-state PNP systems (1) to quasi-one-dimensional models was first proposed in [40] and was rigorously justified in [36] for special cases. A quasi-one-dimensional steady-state PNP model takes the form where X ∈ [0, l] is the coordinate along the axis of the channel, A(X) is the area of crosssection of the channel over the location X.

Electrochemical Potential
The electrochemical potential µ k (X) for the ith ion species consists of the ideal component µ id k (X) and the excess component µ ex k (X): µ k (X) = µ id k (X) + µ ex k (X), where µ id k (X) = z k eΦ(X) + k B T ln with some characteristic number density C 0 defined by The cPNP system takes into consideration of the ideal component µ id k (x) only. This component reflects the collision between ion particles and the water molecules. It has been accepted that the cPNP system is a reasonable model in, for example, the dilute case under which the ion particles can be treated as point particles and the ion-to-ion interaction can be more or less ignored. The excess chemical potential µ ex k (x) accounts for the finite size effect of ions. We point out that, among many limitations, such as the "gating" phenomena, may not be captured by the simple cPNP model. However, the basic findings on dynamics of ionic flows and their dependence on the system parameters, particularly, the permanent charges, the channel geometry, the ratios of boundary concentrations, and the ratios of diffusion constants, provide important insights into the mechanism of ion channels and better understandings of ionic flow properties. More importantly, some of them are non-intuitive, and deserve further studies. More structural detail and more correlations between ions should be taken into considerations in PNP models such as those including various potentials for ion-to-ion interaction accounting for ion size effects ( [5,21,[41][42][43][44] etc.).

Permanent Charge
The spatial distribution of side chains in a specific channel defines the permanent charge of the channel. While some information may be obtained by ignoring the permanent charge and focusing on the effects of boundary conditions, the valences and sizes of ions, etc., we believe that different channel types differ mainly in the distribution of permanent charge [3]. To better understand the importance of the relation of ionic flows and permanent charges, we remark 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 [45,46]. One may roughly understand in the following sense, doping provides the charges that acid and basic side chains provide in a protein channel. 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 are viewed as the reservoirs where there is no permanent charge). In this work, we take the following model for Q(x) where X 0 = 0, X 3 = l and Q 0 is some nonzero constant.

Comparison with Some Existing Works
The current work follows a similar dynamical system framework as that employed in [6,7,[9][10][11] to establish the existence and uniqueness result of the problem. However, compared to these works, our set-ups are much more challenging and more realistic, more importantly, the specific structure of our model allows us to obtain detailed description of the nonlinear interplay among different system parameters. This is far beyond the existence and uniqueness result. To be specific, our model includes three ion species, two positively charged with the same valences, and one negatively charged (in [7,9,10], only two oppositely charged particles are included, selectivity of cations, one of the most relevant biological properties of ion channels cannot be described); and a profile of nonzero but small permanent charges (in [6], it includes three ion species but with zero permanent charges, the effects on ionic flows from the two key structures of ion channels: channel geometry and distribution of permanent charges, cannot be examined, while this could provide crucial insights for the selectivity of cations through membrane channels). In [11], the author extended the work in [7] and established the existence and local uniqueness of the classical PNP system with n ion species.
The work, in some sense, is motivated by [9], and there are some similarities in the treatment. More precisely, both of the works employ regular perturbation analysis to derive the explicit expressions of the individual fluxes up to the first order in the small permanent charge, which is reflected in Section 2.2 in the current work. However, the derivation and the following analysis is much more challenging due to the nonlinearity of the individual fluxes in the potential V (in [9], the individual fluxes are linear in the potential V). The nonlinearity of the individual fluxes in the potential provides much more rich dynamics of ionic flows, and demonstates more complicated nonlinear interaction among the system parameters, which is addressed in Section 3.1. Meanwhile, this indicates that our work provides a better understanding of the mechanism of ionic flows through single ion channels, which is necessary and important for future studies of ion channel problems.

Main Results
For convenience, we briefly summarize our main results as follows with j, k = 1, 2, 3. (i) Constructing a singular orbit of the limiting PNP system (ε → 0) over the whole ; see Theorems 8 and 9 in Section 3.1.3.

Remark 1.
In (v), there are actually another three cases: . The results and arguments are similar to those corresponding to the case stated in (vi-1)-(vi-3), and are not included in this work. Interested readers can study them following our discussions detailed in Section 3.1.

Problem Set-Up
For definiteness, we will take the following setting in this work: (A1).We consider three ion species (n = 3) with z 1 = z 2 = z > 0 and z 3 < 0. (A2).The permanent charge is defined as in (6). (A3).For the electrochemical potential µ i , we only include the ideal component µ id i as in (4). (A4).We assume ε r (X) = ε r and D i (X) = D i .
In the sequel, we will assume (A1)-(A4). We first make a dimensionless rescaling following [3]. With C 0 given in (5), let The BVP (2) and (3) then becomes (noting that z 1 = z 2 = z) , with the boundary conditions, for i = 1, 2, 3, We comment that the dimensionless parameter ε defined in (1.6) as ε = 1 l ε r ε 0 k B T e 2 C 0 is directly related to the ratio κ D /l, where κ D = ε r ε 0 k B T ∑ j (z j e) 2 C j is the Debye length; in particular, ε = κ D /l when z 2 j = 1 and C j = C 0 . Typically, the parameter ε is small due to the fact that the two variables l, the length of the channel, and C 0 , some characteristic number density could be very large. For many cases, the value of ε is of order O(10 −3 ) (see [47] for a more detailed description). (8) and (9) We first rewrite system (8) into a standard form for singularly perturbed systems and convert the boundary value problem (8) and (9) to a connection problem.

Geometric Singular Perturbation Theory for
Upon introducing u = εφ and τ = x. System (8) becomes where overdot denotes the derivative with respect to the variable x. System (10) will be treated as a singularly perturbed system with ε as the singular parameter. Its phase space is R 9 with state variables (φ, u, c 1 , c 2 , c 3 , J 1 , J 2 , J 3 , τ).
For ε > 0, the rescaling x = εξ of the independent variable x gives rise to where prime denotes the derivative with respect to the variable ξ.
We comment that for ε > 0, systems (10) and (11) have exactly the same phase portrait. However, their limiting systems at ε = 0 are different. The limiting system of (10) is called the limiting slow system, whose orbits are called slow orbits or regular layers. The limiting system of (11) is the limiting fast system, whose orbits are called fast orbits or singular (boundary and/or internal) layers. By a singular orbit of system (10) or (11), we mean a continuous and piecewise smooth curve in R 9 that is a union of finitely many slow and fast orbits. Very often, limiting slow and fast systems provide complementary information on state variables. Correspondingly, the main task of singularly perturbed problems is to patch the limiting information together to form a solution for the entire ε > 0 system. Let B L and B R be the subsets of the phase space R 9 defined by Then the original boundary value problem is equivalent to a connection problem, namely, finding a solution of (10) or (11) from B L to B R (see, for example, [48]).
Due to the jumps of the permanent charge function (6) These eight unknown values will be determined along our construction of a singular orbit on the whole interval [0, 1]. (i) The singular orbit on [0, a] consists of two boundary layers Γ 0 l and Γ a l and one regular layer Λ l with (φ, c 1 , c 2 , c 3 , τ) being 3 , a) at x = a.
In particular, given (φ [a] , c 3 ), the flux densities J l k and the value u l (a) are uniquely determined (see Section 2.1.1). (ii) The singular orbit on [a, b] consists of two boundary layers Γ a m and Γ b m and one regular layer Λ m with (φ, c 1 , c 2 , c 3 , τ) being In particular, given (φ [a] , c 3 ), the flux densities J m k and the value u m (a) and u m (b) are uniquely determined (see Section 2.1.2). (iii) The singular orbit on [b, 1] consists of two boundary layers Γ b r and Γ 1 r and one regular layer Λ r with (φ, c 1 , c 2 , c 3 , τ) being In particular, given (φ [b] , c 3 ), the flux densities J r k and the value u r (b) are uniquely determined (see Section 2.1.3).
To obtain a singular orbit on the whole interval [0, 1], one need This consists of eight conditions. The number of conditions is exactly the same as the number of unknowns in (13).
The singular orbit constructed for problem (10) associated to B L and B R consists of nine pieces Γ 0 Once a singular orbit is constructed, one then can apply the geometric singular perturbation theory, such as Exchange Lemma, to show that, for ε > 0 but small, there is a unique solution that is close to the singular orbit.
In this work, we will examine the competition between cations due to the nonlinear interaction among channel geometry, small permanent charges, diffusion coefficients and boundary conditions, which can be extracted from the matching conditions (14) (for simplicity, we still use (φ, u, c 1 , c 2 , c 3 , J 1 , J 2 , J 3 ) for the zeroth order system in ε).
We remark that u l (a), u m (a), u m (b), u r (b), J l k , J m k , J r k are actually the functions of the unknowns φ [a] , c k with the parameter Q 0 . Furthermore, for simplicity, in the following analysis, we will use J k to denote J l k , J m k and J r k , respectively. Once a solution for (φ [a] , c 3 ) is determined, one then can derive the zeroth order (in ε) individual fluxes J k (Q 0 ) = D k J k (Q 0 ). Through our following discussions, we always assume the so-called electroneutrality boundary concentration conditions For simplicity, we also introduce Following (13), we introduce 3 , J 1 , J 2 , J 3 , a) ∈ R 9 : u, J 1 , J 2 , J 3 , arbitrary .
We now construct a singular orbit on [0, a] that connects B L to B a , which generally consist of two boundary layers and a regular layer (see [7,10,11,43]). Over the subinterval [0, a], the permanent charge is zero because we review [0, a] as one of the reservoirs. However, the nonzero Q 0 over the subinterval [a, b] will affect the solution on [0, a] and [b, 1] (another reservoir with zero permanent charge) through the matching conditions 3 to construct the singular orbit over the whole interval [0, 1].
Limiting fast dynamics and boundary layers on [0, a] Setting ε = 0 in (10), we get the so-called slow manifold, Setting ε = 0 in (11), we get the limiting fast system, Note that the slow manifold Z l is the set of equilibria of (17). The following can be established directly [11]. Lemma 1. For system (17), the slow manifold Z l is normally hyperbolic.
Proof. The slow manifold Z l is precisely the set of equilibria of (17). The linearization of (17) at each point of (φ, 0, c 1 , c 2 , c 3 , J 1 , J 2 , J 3 , τ) ∈ Z l has seven zero eigenvalues whose generalized eigenspace is the tangent space of the seven-dimensional slow manifold Z l of equilibria, and the other two eigenvalues are ± z 2 1 c 1 + z 2 2 c 2 + z 2 3 c 3 , whose eigenvectors are not tangent to Z l (Recall that c i 's are concentrations and we are only interested in positive ones). Thus, Z l is normally hyperbolic.
The theory of normally hyperbolic invariant manifolds [49] states that there exists eight-dimensional stable manifold W s (Z l ) of Z l that consists of points approaching Z l in forward time; and there exists eight-dimensional unstable manifold W u (Z l ) of Z l that consists of points approaching Z l in backward time. Let M 0 l be the collection of orbits from B L in forward time under the flow of system (17) and M a l be the collection of orbits from B a in backward time under the flow of system (17). Then, for a singular orbit connecting B L to B a , the boundary layer Γ 0 In this subsection, we will determine the boundary layers N 0 l and N a l , and their landing points ω(N 0 l ) and α(N a l ) on the slow manifold Z l . The regular layer, determined by the limiting slow system, will lie in Z l and connect the landing points ω(N 0 l ) at x = 0 and α(N a l ) at x = a. First, one has Proposition 1. The following functions are the first integrals of system (17), Proof. This can be verified directly.
For the landing points ω(N 0 l ) and α(N a l ), following the similar outline as those in [7,11], one has Proposition 2. Assume the condition (15), one has (i) The stable manifold W s (Z l ) intersects B L transversally at points V, u 0 r , L 1 , L 2 , L 3 , J 1 , J 2 , J 3 , 0 , and the ω-limit set of N 0 where J i for i = 1, 2, 3 are arbitrary, and where J i for i = 1, 2, 3 are arbitrary, and for k = 1, 2, (iii) The boundary layer Γ 0 l at x = 0 is determined up to (J 1 , J 2 , J 3 ) as follows: the φ-component satisfies the Hamiltonian system Similarly, the boundary layer Γ a l at x = a is determined in the following way: the φ-component satisfies the Hamiltonian system Limiting slow dynamics and regular layers on [0, a] For convenience, we introduce T m , T c and H(x) defined as Next we construct the regular layer on Z l that connects ω(N 0 r ) and α(N a l ). Note that, for ε = 0, system (10) loses most information. To remedy this degeneracy [7,10,11], we make a rescaling u = εp and −zc 1 − zc 2 − z 3 c 3 = εq in system (10). In term of the new variables, system (10) becomeṡ It is again a singular perturbation problem and its limiting slow system iṡ For system (20), the slow manifold is . Therefore, the limiting slow system on S l is given bẏ For system (21), one has where J 1 , J 2 and J 3 are uniquely determined as Proof. Adding the second equation to the third one in (21), one haṡ which gives Substituting (24) into the second equation in (21) to geṫ .
By the variation of constants formula, we obtain Similarly, c 2 (x) can be obtained. Substituting (24) into the first equation in (21) to geṫ . Evaluating c 1 (x), c 2 (x) and φ(x) at x = a yield the formulas for J 1 , J 2 and J 3 .
The slow orbit given in Lemma 2 connects ω(N 0 r ) and α(N a l ). LetM 0 r (resp.,M a l ) be the forward (resp., backward) image of ω(N 0 r ) (resp., α(N a l )) under the slow flow (21). One has the following result (the proof follows exactly the same line as Proposition 3.7 in Section 3.1.2 of [43]).

Proposition 3.
On the nine-dimensional slow manifold S l ,M 0 r andM a l intersect transversally along the unique orbit Λ l (x) given in (25).
Following (13), we let and the corresponding limiting fast system, For system (26), similar argument shows that the slow manifold Z m is normally hyperbolic. We denote the stable (resp. unstable) manifold of Z m by W s (Z m ) (resp. W u (Z m )). Let M a m be the collection of orbits from B a in forward time under the flow of system (26) and M b m be the collection of orbits from B b in backward time under the flow of system (26). Then, for a singular orbit connecting B a to B b , the boundary layer Γ a m at In this section, we will determine the boundary layers N a m and N b m , and their landing points ω(N a m ) and α(N b m ) on the slow manifold Z m . The regular layer, determined by the limiting slow system, will lie in Z m and connect the landing points ω(N a m ) at x = a and α(N b m ) at x = b. Similarly, we have the following result.
where J i for i = 1, 2, 3 are arbitrary, and φ [a,m] is the unique solution of where c where J i for i = 1, 2, 3 are arbitrary, and φ [b,m] is the unique solution of where c Similarly, the boundary layer Γ b m at x = b is determined in the following way: the φ-component satisfies the Hamiltonian system

Limiting slow dynamics and regular layers on [a, b]
We now turn to the study of the flow in the vicinity of the slow manifold Z m . Following a similar argument as that in Section 2.1.1, we make a scaling u = εp and −zc 1 It is again a singular perturbation problem and its limiting slow system iṡ For system (31), the slow manifold is Therefore, the limiting slow system on S m is given bẏ Note that, on Z m , one has z(c 1 + c 2 ) + z 3 c 3 + Q 0 = 0. It follows that (32) has the same phase portrait as that of the following system obtained by multiplying z(z − z 1 )(c 2 + c 3 ) − z 1 Q 0 h(τ) on the right-hand side of system (32): where A 3 (y) = Q 0 zT m 1 − e zz 3 T m y − C [a,m] J 1 +J 2 e zz 3 T m y , and J 1 , J 2 and J 3 are uniquely determined as, for some y 0 > 0, Remark 2. The proof is similar as that of Lemma 2. As for the system (35), note that we are looking for solutions to reach α(N a m ); that is, whenever τ(y) = b, we require φ(y) = φ [b,m] , c 2 (y) = c . Evaluating system (34) at y = y 0 , by a careful calculation, one has system (35).
It follows that the regular layer solution Λ m on [a, b] is given by (34) with J 1 , J 2 and J 3 determined by (35). Together with the boundary layers Γ a m and Γ b m described in statement (iv) of Proposition 4, this gives the singular orbit on the interval [a, b].

Singular Orbit on
Th construction of singular orbits over [b, 1] from B b to B R is virtually the same as the construction of singular orbits on [0, a] studied in Section 2.1.1. Instead of repeating the same process, we will state only the results for later use.
Limiting fast dynamics and boundary layers on [b, 1] The limiting fast system reads The slow manifold is Z r = {u = 0, z(c 1 + c 2 ) + z 3 c 3 = 0}, which consists of the equilibria of system (36) and is normally hyperbolic with a eight-dimensional center-stable manifold W s (Z r ) and a eight-dimensional center-unstable manifold W u (Z r ). For the boundary layers, one has the following proposition.

Proposition 5.
Under the condition (15), one has (i) System (36) has the following integrals: 3 , J 1 , J 2 , J 3 , b , and the ω-limit set of N b where J i for i = 1, 2, 3 are arbitrary, and where J i for i = 1, 2, 3 are arbitrary, and (iv) The boundary layer Γ b r at x = b is determined up to (J 1 , J 2 , J 3 ) as follows: the φ-component satisfies the Hamiltonian system Similarly, the boundary layer Γ 1 r at x = 1 is determined in the following way: the φ-component satisfies the Hamiltonian system Limiting slow dynamics and regular layers on [b, 1] We now examine the existence of the regular layer on Z r that connects ω(N b r ) and α(N 1 r ). It follows from exactly the same analysis as that in Section 2.1.1, the limiting slow dynamics readṡ For system (38), one has are given in Proposition 5. It is given by , and J 1 , J 2 and J 3 are uniquely determined as , The regular solution Λ r that connects ω(N b r ) and α(N 1 l ), together with Γ b r and Γ 1 l in statement (iv) of Proposition 5 yields the singular orbit on [b, 1].

Matching and Singular Orbits on the Whole Interval [0, 1]
A singular orbit over the whole interval [0, 1] will be the union of the singular orbits constructed on each of the subintervals (see Figure 1). The matching conditions are u l (a) = u m (a), u m (b) = u r (b), and J 1 , J 2 and J 3 have to be the same on all subintervals; that is, from (18), (23), (27)-(30), (35), (37) and (39), , , where, for k = 1, 2, Recall that (φ [a] , c 3 ) and (φ [b] , c 3 ) are the unknown values preassigned at x = a and x = b, J 1 , J 2 and J 3 are the unknown values for the flux densities of the three ion species. There are also three auxiliary unknowns φ [a,m] , φ [b,m] and y 0 in (40). The total number of unknowns in (40) is fourteen, which matches the total number of equations.
A qualitative important question is whether the set of nonlinear Equation (40) (a governing system) has a unique solution. This can be studied through bifurcation analysis and numerical simulations. However, this is beyond the aim of this work.

Existence of Solutions near the Singular Orbit
Note that any solution of the set of algebraic equations determines a singular orbit for the connection problem. Once a singular orbit is constructed, one can apply geometric singular perturbation theory to show that, for ε > 0 small, there is a unique solution that is close to the singular orbit.
For our case, the singular orbit consists of nine pieces: two boundary layers Γ 0 l and Γ 1 r ; four internal layers Γ a l , Γ a m , Γ b m and Γ b r ; and three regular layers Λ l , Λ m and Λ r (see Figure 1). More precisely, with J = (J 1 , J 2 , J 3 ), 1.
The boundary layer Γ 0 l connects the point V, u 0 r , L 1 , The regular layer Λ 1 connects the point The internal layer Γ a l connects the point φ [ The internal layer Γ b r connects the point φ [b] , u m (b), c The regular layer Λ r connects the point φ [b,r] , 0, c The boundary layer Γ 1 r connects the point φ R , 0, c R 1 , c R 2 , c R 3 , J, 1 ∈ α(N 1 r ) ⊂ Z r to the point 0, u 1 l , R 1 , R 2 , R 3 , J, 1 ∈ B R . Figure 1. Schematic picture of a singular orbit projected to the space of (u, ∑ z k c k , τ) with Q(x) defined in (6).
r be the singular orbit of the connecting problem system (10) associated with B L and B R in system (12). There exists ε 0 > 0 small, so that if 0 < ε < ε 0 , then the boundary value problem (8) and (9) has a unique solution near the singular orbit Γ 0 Proof. Fix δ > 0 small to be determined. Let To track this evolution of M 0 l (ε), we apply the exchange lemma successively (three times) along the stages in order described above. During the first stage, we track the evolution of M 0 l (ε) along the singular orbit Γ 0 l ∪ Λ l ∪ Γ a l . The Exchange Lemma ( [48,50,51], etc.) indicates that , at the end of the first stage and near Γ a ) for some ρ > 0 independent of ε, provided that the following conditions are satisfied: (i) M 0 l (0) ∩ W s (Z l ) is transversal along Γ 0 l , which is established in Proposition 2; (ii) the vector field on Z l is not tangent to ω(N 0 l ) at (φ L , 0, c L 1 , c L 2 , c L 3 , J l 1 , J l 2 , J l 3 , 0) ∈ Z l , which follows fromτ = 1 in (21). Let Σ l = W u (α(N a l ) · (−ρ, ρ) ∩ {τ = a}. then, near Γ a l , M 0 l (ε) is close to the forward trace of Σ l under the flow of system (10) with Q = Q 0 . We can then apply the Exchange Lemma again to M l (ε) along Γ a m ∪ Λ m ∪ Γ b m 0ver [a, b]. At the end of this stage, one has that M l (ε) is is close to the forward trace of Σ m under the flow of system (10) with Q = 0. We can then apply the Exchange Lemma again to 1]. At the end of this stage, one has that M l (ε) is is a singleton. This completes the proof.

Regular Perturbation Analysis: Expansions along Small Q 0
We expand all unknown quantities in the governing system (40) and (41) in Q 0 under the assumption that |Q 0 | is small, for example, for k = 1, 2, 3, we write We will determine the coefficients of the zeroth order and first order terms for dominating effects on ionic flows from the permanent charge. We introduce and, corresponding to (16) and (19), 2i , i = 0, 1.
Careful calculations (tedious but straightforward) give Proposition 6. The zeroth order solution in Q 0 of (40) and (41) is given by

Proposition 7.
The first order solution in Q 0 of (40) and (41) is given by, with k = 1, 2, 0 (ln L − ln R) For convenience, we introduce f 0 (L, R) and g 0 (L, R; V), which are defined by

Lemma 5. One has
From Propositions 6 and 7, one has Corollary 1. Assume the electroneutrality boundary conditions (15). For k = 1, 2, one has For convenience in our following analysis, we introduce a function γ(t) for t > 0 with For γ(t), one can easily established the following properties.

Lemma 7.
Set t = L/R. A has the same sign with that of R − L, that is, if t > 1, then A < 0, and if t < 1, then A > 0.
• For the case with α < γ(t), we first claim that there exists a unique β 1 ∈ (α, 1) such that 1 − B = 0 for β = β 1 . In fact, based on the facts that lim β→α (1 − B) = (α − γ(t)) ln t < 0 and 1 − B is strictly increasing on (α, +∞), it suffices to show that 1 − B > 0 for β = 1, which follows from g(1) > 0. For convenience, for t > 1, we set where J 0 1,2 = D 1 J 10 − D 2 J 20 and J 1 1,2 = D 1 J 11 − D 2 J 21 are given by where L d and R d is defined in (16). In particular, for B = 1, one has J 1 1,2 is the leading term that contains the effect from small permanent charges, and will be our main focus. The study on J 1 1,2 could provide important insights and better understanding of the selectivity of cations, and meanwhile demonstrates the critical role played by the small permanent charge in the selectivity of ion channels. The two expressions for J 1 1,2 (V) will be alternatively chosen for convenience in our following discussion.
In terms of the interplay between the diffusion coefficients and boundary concentrations, there are six cases to discuss, more precisely, one has In this work, we will just focus on case (i), (iii) and (v), since the discussions for case (ii), (iv) and (vi) are very similar. Interested readers can work on them following our arguments. In our following discussion, we examine J 1 1,2 (V) from two directions for each case: the sign of J 1 1,2 (V) and the monotonicity of J 1 1,2 (V) in the electric potential V.
In this section, we study J 1,2 (V) under the condition D 1 We first consider a special case.
), that is, the (small) positive Q 0 reduces (resp. enhances) J 1,2 . Furthermore, and B = 1, directly from (49), one has From Lemmas 5 and 7, together with z − z 3 > 0 and H(1) > 0, one has J 1 1,2 < 0 (resp. ). To see the monotonicity of J 1 1,2 (V) in the potential V, we calculate where h 0 (V) = L − Re −zV − L(ln L − ln R + zV). Direct calculation shows that h 0 (V) is concave downward for all V and attains its global maximum 0 at V 1 = 1 z ln R L . Note that lim . One then has ). This completes our proof.
We next consider the more general case with B = 1.
Proof. The statements follow directly from (49) and Lemma 5. Proof. We will just provide a detailed proof for the first statement, and the second one can be discussed similarly.
It follows that , and further It is easy to see that f d (V) has two zeros given by V 1 and This holds for the case with V 1 ≥ V d by a similar argument. Note that lim

Case Study with
We study the term J 1 1,2 (V), which provides information of the preference of the ion channel over different cation under the condition D 1 We first characterize the sign of A(1 − B) under the further restriction D 1 , and the order of the critical potentials V 2 and V 3 .
We comment that the proofs of Corollaries 2-4 can be directly discussed based on Lemmas 7 and 8. We skip the details here.

Case Study with
In this section, we study the term J 1 1,2 (V) for the case with R 2 , which is equivalent to L d < 0 and R d > 0. By similar arguments as those in Sections 3.1.1 and 3.1.2, we obtain the following results.
3.1.4. Studies on the Magnitude of J 1,2 (V) For convenience in the argument, we let S 1 represent the first cation corresponding to the flux J 1 and S 2 be the cation corresponding to the flux J 2 .
Recall that, with Q 0 > 0 small, one has J 1,2 (V, , in Sections 3.1.1-3.1.3, we analyze the leading term J 1 1,2 (V) that contains the effects from small permanent charges, in particular, (i) the sign of J 1 1,2 (V), which characterizes the small positive permanent charge effects on the competition between two cations. To be specific, if J 1 1,2 (V) > 0 (resp. J 1 1,2 (V) < 0), then, the small positive permanent charge enhances (resp. reduces) J 1,2 (V; Q 0 ), and in either way, it affects the preference of the ion channel over different cation, which is closely related to the selectivity phenomena of the ion channel.
(ii) the monotonicity of J 1 1,2 (V) in the electric potential V, which further reduces or strengthens the preference by adjusting/controlling the boundary membrane potential.
Taking the case dJ 1 1,2 (V) dV > 0 for example, if J 1 1,2 (V) > 0, then, one can increase the boundary electric potential to further strengthen the individual flux J 1 (V), which indicates that more cation S 1 will go through the ion channel; while if J 1 1,2 (V) < 0, one then need to decrease the boundary electric potential for more cation S 1 to go through the ion channels.

Numerical Simulations
In this part, numerical simulations are performed to provide more intuitive illustrations of some analytical results. To get started, we rewrite the system (8) and (9) as a system of first order ordinary differential equations. Upon introducing u = εφ, one has with boundary conditions In our simulations to system (50) and (51), we take

Remark 3. The choice of h(x)
is based on the fact that the ion channel is cylindrical-like, and the variable cross-section area is chosen to reflect the fact that the channel is not uniform and much narrower in the neck (a < x < b) than other regions [9]. We further take r 0 = 0.5 and the function h(x) is then continuous at the jumping points x = a and x = b. Different models for h(x) may be chosen, and similar numerical results should be obtained.
Our numerical simulations consist of two parts focusing on some cases discussed in Sections 3.1.1-3.1.3, respectively, and further verify some analytical results stated in Theorems 3, 4, 6, 7 and 9 for some carefully selected system parameters. Other related results can also be numerically illustrated by choosing different parameter values, and we leave that to interested readers.
Recall that It turns out that our numerical simulations with nonzero but small ε are consistent with our analytical results. To be specific, (i) By choosing L 1 = 24, L 2 = 6, R 1 = 9, R 2 = 2, a = 1/3 and b = 0.7, one has α =0.265139, β = 0.751381, It follows that which is consistent with the conditions stated in (ii) of Lemma 8, and indicates that 1 − B < 0. Together with Lemma 7, we have A(1 − B) > 0. Therefore, A(1 − B)R d > 0, which is consistent with the condition (i) stated in Theorems 3 and 4, respectively.
Our numerical results show that (see Figure 4) J 1 1,2 (V; ε) has one zero V 2 and two critical points V 3 c and V 4 c with V 3 c < V 4 c such that J 1 1,2 (V; ε) > 0 if V > V 2 and J 1 1,2 (V; ε) < 0 if V < V 2 ; furthermore, J 1 1,2 (V; ε) increases in the potential V if V < V 3 c , decreases in the potential V if V 3 c < V < V 4 c , and increases in the potential V if V > V 4 c . We remark that the numerical result is consistent with our analytical result stated in (i) of Theorem 9.

Concluding Remarks
We study a one-dimensional steady-state Poisson-Nernst-Planck system with three ion species, two cations having the same valence and one anion. Nonzero but small permanent charge, the major structural quantity of an ion channel, is included in the model. The cross-section area of the channel is included in the system, and it provides certain information of the geometry of the three-dimensional channel, which plays crucial roles in our analysis. Two specific structures of the PNP model (i) the existence of a complete set of first integrals (Proposition 1, first statement (i) in Proposition 4 and Proposition 5, respectively); (ii) the observation made in Section 2.1.2 allows one to transform the limit slow system (32) to a linear system (33) with constant coefficients; allows one to reduce the singularly perturbed boundary value problem to an algebraic system-the governing system (40). The significance of the governing system is: (i) it includes almost all relevant physical parameters; (ii) once a solution of the governing system is obtained, the singular orbit (the zeroth order approximation (in ε) solution of the boundary value problem) can be readily determined.
Based on these specific structures of this concrete model, under the framework of the geometric singular perturbation theory, a singular orbit is obtained, from which explicit expressions of J k0 and J k1 are extracted. This makes it possible for one to further analyze the competition between cations, which further depend on the complicated nonlinear interplays among system parameters, such as the diffusion coefficients (D 1 , D 2 ), the channel geometry in terms of (α, β), the boundary conditions (L k , R k ; V), k = 1, 2, 3, particularly, the interaction among D 1 D 2 , R 2 R 1 and L 2 L 1 plays a critical role in characterizing the competition between cations (Sections 3.1.1-3.1.3). Among others, we find (i) As functions of the membrane potential V (fixing other system parameters), (i1) J 1 1,2 = D 1 J 11 − D 2 J 21 and ∂ V J 1 1,2 can be positive (resp. negative), which further depends on the nonlinear interaction among boundary concentrations and diffusion coefficients. The sign of J 1 1,2 provides critical information related to the selectivity phenomena of ion channels, while the sign of ∂ V J 1,2 provides efficient ways to further adjust/control the preference of the ion channel over distinct cation (Characterized in Theorems 2-9); (i2) |J 1,2 |, the magnitude of J 1,2 , which is equivalent to J 0 1,2 J 1 1,2 for small positive Q 0 , is analyzed, which provides further information of the ion channel over distinct cations (Theorems 10 and 11).
(ii) Critical potentials that balance the small permanent charge effects on J 1 1,2 (such as V 2 and V 3 ) are identified (Definition 1). Those critical potentials can be experimentally estimated. Taking V 2 for example, one can take an experimental curve as J 1,2 (V; Q 0 ) and numerically (or analytically) compute J 0 1,2 (V; 0) for ideal case that allows one to get an estimate of V 2 by considering the zeros of J 1,2 (V; Q 0 ) − J 0 1,2 (V; 0). The critical potentials play critical roles in characterizing permanent charge effects on ionic flows through membrane channels.
Finally, we comment that the setup in this work is relatively simple, and may raise the concern about the feasibility. Indeed, cPNP is known to be reliable when the ionic mixture is dilute, but with more ion species and nonzero permanent charges included, the ionic mixture would be crowded. On the other hand, the setup is reasonable for semiconductor problems and for synthetic channels. The detailed analysis in this work is not only important in understanding the behavior of such fluxes across channels but can also better help in understanding the fundamental mechanistic principles of metabolite fluxes in engineered synthetic biological channels [52]. Furthermore, the study in this work is the first step for analysis of more realistic models. The simple model considered in this work allows us to obtain more explicit expressions of the ionic fluxes in terms of physical parameters of the problem so that we are able to extract concrete information on small permanent charge effects. Moreover, the analysis in this simpler setting provides important insights for the analysis and numerical studies of more realistic models.