Dynamics of ionic flows via Poisson-Nernst-Planck systems with local hard-sphere potentials: Competition between cations

We study a quasi-one-dimensional steady-state Poisson-Nernst-Planck type model for ionic flows through a membrane channel with three ion species, two positively charged with the same valence and one negatively charged. Bikerman’s local hard-sphere potential is included in the model to account for ion sizes. The problem is treated as a boundary value problem of a singularly perturbed differential system. Under the framework of a geometric singular perturbation theory, together with specific structures of this concrete model, the existence and uniqueness of solutions to the boundary value problem for small ion sizes is established. Furthermore, treating the ion sizes as small parameters, we derive an approximation of individual fluxes, from which one can further study the qualitative properties of ionic flows and extract concrete information directly related to biological measurements. Of particular interest is the competition between two cations due to the nonlinear interplay between finite ion sizes, diffusion coefficients and boundary conditions, which is closely related to selectivity phenomena of open ion channels with given protein structures. Furthermore, we are able to characterize the distinct effects of the nonlinear interplays between these physical parameters. Numerical simulations are performed to identify some critical potentials which play critical roles in examining properties of ionic flows in our analysis.


Introduction
One of the fundamental concerns of physiology is the function of ion channels. Ion channels are approximately cylindrical, hollow proteins with a hole down their middle that provides a controllable path for electrodiffusion of ions (mainly Na + , K + , Ca ++ and Cl − ) through biological membranes, es-tablishing communications among cells and the external environment. This way, ion channels control a wide range of biological functions. The study of ion channels consists of two related major topics: Structures of ion channels and ionic flow properties. With a given structure of an open channel, the main interest is to understand its electrodiffusion property.
Beyond general electrodiffusion phenomema for electrolytic solutions in bulks or near charged walls, ionic flows through ion channels have more specifics, that is, the study of ionic flows has to take into considerations of global constraints, including boundary conditions in addition to the structure of an ion channel. As already demonstrated by the celebrated works [51][52][53][54][55] of Hodgkin and Huxley for neurons consisting of a population of ion channels and by works in the volume Single-Channel Recording [94] (edited by B. Sakmann and E. Neher) and many works afterwards, the properties of ion channels depend in an extremely rich way on different regions of boundary concentrations and boundary electric potentials. It is the global constraints and the internal structures of membrane channels that make the relevant electrodiffusion properties specific for ion channel problems.
Taking the structural characteristics into consideration, the basic continuum model for ionic flows is the Poisson-Nernst-Planck (PNP) system, which treats the aqueous medium (within which ions are migrating) as a dielectric continuum [9, 10, 14, 16, 19, 27-31, 40-43, 59, 60, 68, 93]: (1.1) where X ∈ Ω with Ω being a three-dimensional cylindrical-like domain representing the channel, Q(X) is the distribution of the permanent charge along the interior wall of the channel, ε r (X) 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, for the ith ion species, C i is the concentration, z i is the valence (the number of charges per particle), J i is the flux density, D i (X) is the diffusion coefficient, and µ i (X) is the electrochemical potential. Under some reasonable conditions, the PNP system can be derived as a reduced model from molecular dynamics [95], from Boltzmann equation [3], and from variation principles [56][57][58]. The simplest PNP model is the classic PNP (cPNP) system that contains only the ideal component of electrochemical potential. The cPNP system treats ions essentially as point-charges, and neglects ion size effects. It has been simulated [13,14,16,19,40,48] and analyzed [1, [4][5][6]32,39,62,69,72,73,85,97,98,102,104,105] to a great extent. However, a major weak point of the cPNP model is that ions are treated as point of charges, which is only reasonable in the extremely dilute setting. Many extremely important properties of ion channels, such as selectivity, rely on ion sizes critically, in particular, for ions that have the same valence, such as sodium Na + and potassium K + , the main difference is their ionic sizes. Selectivity of ion channels has been a central issue in the study of cells and tissues at least since Hodgkin, Huxley, and Katz discovered the role of Na + and K + currents in the action potential [18,49,50]. Selectivity of ion channels has been theorized to occur is in a variety of ways, ranging from structural discussions [23-26, 79, 80, 108, 109] to scaling models [33], kinetic models [45,46] and electrostatic models [20][21][22]. However, all these different theories of selectivity fail to reproduce experimental measurements over a wide range of conditions.
On the contrary, reduced models, such as the PNP model, using only physical variables, which include the radius of the pore (the channel geometry) and the dielectric constant of the protein, are successful in interpreting the mechanisms that govern selectivity of the L-type Ca ++ , Na + channels and also their mutations [9][10][11][12]. It has come as a surprise that a model of selectivity that includes only a few features of the atomic structure has been able to describe the selectivity properties of calcium and sodium channels very well, in all solutions over a wide range of conditions, with only two adjustable parameters using crystal radii of ions. Examining the finite ion size effects on ionic flows, especially for multiple cations with the same valence should provide deep insights into the selectivity phenomenon of ion channels, and this is the motivation of our work. Ionic flows through ion channels exhibit extremely rich phenomena, which is why ion channels are nano-scale valves for essentially all activities of living organisms. This is the very reason that it is a great challenge to understand the mechanisms of ion channel functions. To study effects on ionic flows from finite ion sizes, one has to consider excess (beyond the ideal) component in the electrochemical potential. One way is to include uncharged hard-sphere potentials to partially account for ion size effects. Physically, this means that each ion is approximated as a hard-sphere with its charges at the center of the sphere. Both local and nonlocal models for hard-sphere potentials were introduced for this purpose [8,91,92]. Nonlocal models give the hard-sphere potentials as functionals of ion concentrations while local models depend pointwise on ion concentrations. The PNP models with ion sizes have been investigated analytically [2,7,61,63,71,76,77,99,106] for two ion species, one positively charged and one negatively charged and computationally [41-43, 47, 56-58, 68, 84, 88, 96, 107] for ion channels and have shown great success. Existence and uniqueness of minimizers and saddle points of the free-energy equilibrium formulation with ionic interaction have also been mathematically analyzed [34,70].
For mathematical analysis, the challenge lies in the fact that specific dynamics depend on complicated nonlinear interplays of multiple physical parameters such as boundary concentrations and potentials, diffusion coefficients, ion sizes and valences, permanent charge distributions, etc. In general, there is no hope to have explicit solution formulae for such a complicated problem even with simple boundary values. The recent development in analyzing classic PNP models [32,72,73] sheds some lights on the voltage-current relationship in simplified settings. This development is based heavily on modern invariant manifold theory of nonlinear dynamical systems, particularly, the geometric theory of singular perturbations. Together with specific structures to PNP models, we are able to, far beyond the existence results, obtain a more or less explicit approximation formula for solutions from which one can further study the qualitative properties of ionic flows and extract concrete information directly related to biological measurements.
In this work, we will study the PNP model with three ion species, two positively charged with the same valence and one negatively charged. Bikerman's local hard-sphere model is included to account for finite ion size effects. Of particular interest is the competition between two cations due to finite ion sizes, which provides useful and deep insights into the selectivity phenomena for open ion channels with given protein structures. Our analysis is based on a further reduction of PNP models. On the basis that ion channels have narrow cross sections relative to their lengths, PNP systems defined on three-dimensional ion channels are further reduced to quasi-one-dimensional models first proposed in [82] and, for a special case, the reduction is rigorously justified in [75]. A quasi-one-dimensional (time-evolution) PNP type model for ion flows of n ion species is where A(X) is the area of the cross-section of the channel at location X. The boundary conditions are, for i = 1, 2, · · · , n, For ion channels, an important characteristic is the I-V relation. Given a solution of the boundary value problem (1.2)-(1.3), the current is where z k J k is the individual flux of charge of the kth ion species. For fixed boundary concentrations L k 's and R k 's, J k 's depend on V only and formula (1.4) provides a dependence of the current I(V; ε) on the voltage V. In terms of applications, it is important to study properties of individual fluxes J k due to the fact that most experiments (with some exceptions) can only measure the total current I while individual fluxes contain much more information on channel functions. This is another reason that we mainly focus on the study of individual fluxes in this work. We would like to point out that in general, the I-V relation is not unique for the PNP system with nonzero permanent charge Q(x), even for the classic PNP system with two ion species [32]. However, with the permanent charge Q(x) being zero, the I-V relation is unique, that is, the cPNP system has a unique solution [72]. In our work, we assume the permanent charge to be zero, and view the volume of the ion species as a regular perturbation to the system, and since the linearized system at zero ion-size is nondegenerate, the solution to the PNP system is unique; but the dynamics of the ionic flows is still rich due to the nonlinear interplay between distinct physical parameters. We comment that our results, for the relatively simple setting and assumptions of our model, are rigorous. We believe these results will provide useful insights for numerical and even experimental studies of ionic flows through membrane channels. At the same time, we point out that the quasione-dimensional PNP model and the local hard-sphere model adopted in this paper are rather simple. Aside the fact that they will miss the three-dimensional features of the problem, a major weakness is the omission of the excess electrostatic component in the excess potentials. As a result, important phenomena such as charge inversion and layering may not be detected by this model. The rest of this paper is organized as follows. In section 2, we set up our problem with further assumptions. In section 3, the existence and (local) uniqueness result for the boundary value problem is established under the framework of the geometric singular perturbation theory. Based on the analysis, treating the ion sizes as small parameters, approximations of individual fluxes are derived, from which the effects on ionic flows from ion sizes are analyzed in details. This leads to our main focus, the competition between two cations due to the nonlinear interplay between finite ion sizes, diffusion coefficients and boundary conditions, studied in section 4. Some concluding remarks are stated in section 5. A detailed proof of Lemma 3.5 is provided in Appendix A.

Excess potential and a local hard sphere model
The electrochemical potential µ i (X) for the ith ion species consists of the ideal component µ id i (X) and the excess component with some characteristic number density C 0 . The classical PNP system takes into consideration of the ideal component µ id i (X) only. This component reflects the collision between ion particles and the water molecules.
It has been accepted that the classic PNP 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 i (X) accounts for finite sizes of ions, which is critically important for selectivity−which ion a channel prefers−of ion channel. While channel protein structures play an important role for channel selectivity, the electrodiffusion of ions is clearly another critical component. For different ions with the same valence, their sizes are a crucial factor for the selectivity. The excess potential µ ex i (X) consists of two components: the (uncharged) hard-sphere component µ HS i (X) and the excess electrostatic component µ ES i . While there are very successful models for µ HS i (X), modeling of µ ES i (X) is itself extremely challenging and is still a very active research area in liquid-state theory of chemistry and physics [34,41,90,91].
In this paper, we will take the following Bikerman's local hard-sphere model to approximate µ ex i (X) 1 where ν j is the volume of the jth ion species. In particular, in the following analysis, we will take ν n = ν so that λ n = 1.

The steady-state boundary value problem and assumptions
The main goal of this paper is to examine the qualitative effect of ion sizes via the steady-state boundary value problem of (1.2)-(1.3) with the local hard-sphere model (2.2) for the excess potential. For definiteness, we will take the following settings: (A1) We consider three ion species (n = 3) with z 1 = z 2 = −z 3 = 1. (A2) The permanent charge is set to be zero: Q(X) = 0. (A3) For the electrochemical potential µ i , in addition to the ideal component µ id i , we also include the local hard-sphere potential µ bik i in (2.2). (A4) The relative dielectric coefficient and the diffusion coefficient are constants, that is, ε r (X) = ε r and D i (X) = D i .
Remark 2.1. In the study of ion channel problems, the selectivity of cations Na + and K + is extremely important. This is exactly the reason that we choose z 1 = z 2 = 1 together with z 3 = −1 for the anion Cl − .
In the sequel, we will assume (A1)-(A4). Under the assumptions (A1)-(A4), the steady-state system of (1.2) reads (2. 3) The boundary conditions become, for i = 1, 2, 3, We first make a dimensionless rescaling following [39]. Set C 0 = max{L i , R i : i = 1, 2} and let (2.5) We would like to point out that the dimensionless parameter ε defined in (2.5) 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 )).
The boundary value problem (2.3)-(2.4) becomes with the boundary condition Remark 2.2. We will take h(x) = 1 over the whole interval [0, 1] in our analysis. This is because for ion channels with zero permanent charge, it turns out that the variable h(x) contributes through an average, explicitly through the factor , see [71]), which does not affect our analysis of qualitative properties of ionic flows.

Denote
Then, system (2.6) can be rewritten as Recall that the boundary condition is now 3. Geometric singular perturbation theory for (2.8)-(2.9) We will rewrite system (2.8) into a standard form for singularly perturbed systems and convert the boundary value problem (2.8) and (2.9) to a connecting problem.
Denote the derivative with respect to x by overdot and introduce u = εφ and τ = x. System (2.8) becomes System (3.1) will be treated as a dynamical system of phase space R 9 with state variables (φ, u, C, J, τ).
For ε > 0, the rescaling x = εξ of the independent variable x gives rise to the fast system, where prime denotes the derivative with respect to ξ. 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 connecting problem, namely, finding a solution of (3.1) or In what follows, instead of studying the boundary value problem, we will consider the equivalent connecting problem for system (3.1) or (3.2) and construct its solution from B L to B R . The construction process involves two main steps: the first step is to construct a singular orbit to the connecting problem, and the second step is to apply geometric singular perturbation theory to show that there is a unique solution near the singular orbit for small ε > 0.
Following the idea in [32,[71][72][73], we first construct a singular orbit on [0, 1] that connects B L to B R . In our coming analysis, we assume the so-called electroneutrality conditions (3.4)

Limit fast dynamics and boundary layers
By setting ε = 0 in (3.1), we obtain the so-called slow manifold By setting ε = 0 in (3.2), we get the limiting fast system Note that the slow manifold Z is precisely the set of equilibria of (3.6). It follows from [32,61,63,[71][72][73]76] that Lemma 3.1. For system (3.6), the slow manifold Z is normally hyperbolic.
We would like to point out that the normal hyperbolicity of the slow manifold Z is crucial in our analysis, which guarantees the existence of the stable manifold and the unstable manifold that will be discussed below (see [32,35] for more details for the theory of normally hyperbolic invariant manifolds).
For ν > 0 small, treating (3.6) as a regular perturbation of that with ν = 0, we look for solutions Lemma 3.2. For ν > 0 small, one has, up to the first order in ν, Proof. A direct computation leads to the result. We omit it here.
Substituting (3.7) into system (3.6), together with (3.8), by careful calculations, we obtain the zeroth order limiting fast system in ν: (3.9) and the first order limiting fast system in ν: Let M L be the collection of orbits from B L in forward time under the flow of system (3.6) and M R be the collection of orbits from B R in backward time under the flow of system (3.6). Let W s (Z) be the stable manifold of Z that consists of points approaching Z in forward time, and let W u (Z) be the unstable manifold of Z that consists of points approaching Z in backward time. Then, for a singular orbit connecting B L to B R , the boundary layer at x = 0 must lie in N L = M L ∩ W s (Z) and the boundary layer at x = 1 must lie in N R = M R ∩ W u (Z). Let ω(N L ) be the ω-limit set of N L , and α(N R ) be the α-limit set of N R , which are defined by can be arbitrary, and can be arbitrary, and Proof. The result for system (3.9), the zeroth order system in ν, has been obtained in [6], we will not repeat it here. For system (3.10), the first order system in ν, it has the following four nontrivial first integrals: We now establish the results for φ L 1 , c L 11 , c L 21 , c L 31 and u l 1 for system (3.10). Those for φ R 1 , c R 11 , c R 21 and c R 31 can be established in a similar way. Note that φ(0) = c 11 (0) = c 21 (0) = c 31 (0) = 0. Using the integrals H 1 , H 2 and H 3 , one has Under electroneutrality conditions (3.4), one has c L 10 = L 1 , c L 20 = L 2 and c L 30 = L 3 . It follows that c L 11 = c L 21 = c L 31 = 0. In view of H 4 (0) = H 4 (∞), together with the above analysis, one has u l 1 = 0. This completes the proof.

Limit slow dynamics and regular layer over [0, 1]
We now construct the regular layer Λ on Z that connects ω(N L ) and α(N R ). Following the ideas in [32,[71][72][73], we make a rescaling u = εp and c 1 + c 2 − c 3 = −εq in system (3.1). In term of the new variables, system (3.1) becomeṡ (3.11) where It is again a singular perturbation problem and its limiting slow system iṡ where c 3 = c 1 + c 2 . where c 3 = c 1 + c 2 .
As for the regular layer problem, we look for solutions of (3.13) of the form to connect ω(N L ) and α(N L ) given in Proposition 3.3; in particular, for j = 0, 1, For convenience, we define quantities T c 0 , T c 1 , T m 0 , T m 1 , F 1 , F 2 and a function A(x) as, for k = 0, 1,  For the zeroth order system (3.16) (see [6] for details), one has , and c R 20 are given in Proposition 3.3. It is given by where (3.20) For the first order system (3.17), we have (see Appendix A for a detailed proof) and c R 21 are given in Proposition 3.3. It is given by The next result can be justified directly.
Proposition 3.6. Assume electroneutrality conditions (3.4). Fixing V and treating both the zeroth order approximations J k0 's (in ν) and first order approximations J k1 's (in ν) as functions of boundary concentrations L i and R i , one has J k0 is homogeneous of degree one and J k1 is homogeneous of degree two, that is, for any α > 0, Remark 3.7. This interesting observation actually provides a nice way to adjust the effects on ionic flows (in particular, for the individual flux) from finite ion size by suitably controlling its boundary concentrations. Furthermore, in our later study of the function P(V) = D 1 J 11 − D 2 J 21 , which provides detailed information that is related to the selectivity of cations due to electrodiffusion properties of ions with a given protein structures, this scaling law can either reduce or enhance the preference of the ion channel for cations.

Existence of solutions near the singular orbit
We have constructed a unique singular orbit on [0,1] that connects B L to B R . It consists of two boundary layer orbits Γ 0 from the point and a regular layer Λ on Z that connects the two landing points of the two boundary layers.
Proof. Let ν 0 be as in Proposition 3.8. For 0 ≤ ν ≤ ν 0 , we define u l = u l 0 + u l 1 ν, J 1 (ν) = J 10 + J 11 ν, J 2 (ν) = J 20 + J 21 ν and J 3 (ν) = J 30 + J 31 ν. Fix δ > 0 small to be determined. Let For ε > 0, let M L (ε, δ) be the forward trace of B L (δ) under the flow of system (3.1) or equivalently of system (3.2) and let M R (ε) be the backward trace of B R . To prove the existence and uniqueness statement, it suffices to show that M L (ε, δ) intersects M R (ε) transversally in a neighborhood of the singular orbit Γ 0 ∪ Λ ∪ Γ 1 . The latter will be established by an application of the Exchange Lemma.

Competition between cations
In this section, our main interest is the competition between two positively charged ion species due to the nonlinear interplays between finite ion sizes, diffusion coefficients and boundary conditions. Mathematically, this is characterized by the quantity P = D 1 J 11 − D 2 J 21 (J 11 and J 21 are the leading terms corresponding to two cations that include finite size effects), which provides important information of competitions between two cations. For convenience, in the following, we use S I to denote the cation associated with c 1 and J 1 , and S II to denote the one associated with c 2 and J 2 .

Competition at singularities
For the function P defined in (4.1), direct calculations lead to the following result, which is crucial for our further study on the qualitative properties of ionic flows.  We next consider the competition of two positive charged ion species as L 3 → R 3 and L 3 → R 3 e −V , respectively. Theorem 4.3. Suppose L 1 > R 1 . For small ε > 0 and ν > 0, there exists a critical potential V * > ln 2 such that, for V > V * , one has P(V) > 0 as L 3 → R 3 , that is, the ion channel prefers cation S I over cation S II .

Competitions away from singularities
For fixed boundary concentration R 3 and boundary electric potential V > 0, the removable singularities L 3 = R 3 , L 3 = R 3 e V and L 3 = R 3 e −V split the L 3 -region into four subregions: 0, R 3 e −V , R 3 e −V , R 3 , R 3 , R 3 e V and R 3 e V , +∞ , over which the sign of P can be analyzed. Our study will mainly focus on the sign of P over the first subregion 0, R 3 e −V .
To get started, we establish some results which are crucial to study the sign of P(V), and will be frequently used in our following analysis. Note that, under electroneutrality conditions (3.4), from (3.18), one has where It then follows that Lemma 4.6. Fixing V > 0. Then, Based on Lemma 4.6, we further assume R 3 > 1 in the following argument (the case with R 3 < 1 can be discussed similarly). From the electroneaurality conditions (3.4), one has the following possibilities for L 3 < R 3 e −V : (a) L 1 < R 1 e −V and L 2 < R 2 e −V ; (b) L 1 < R 1 e −V and L 2 > R 2 e −V ; and (c) L 1 > R 1 e −V and L 2 < R 2 e −V . In this section, we will only consider the case (a), which, with V > 0, is equivalent to and study the sign of P, which provides rich dynamical behavior of ionic flows. The other two cases can be analyzed in a similar way, and we leave them to interested readers. From Lemmas 4.5-4.7, we obtain Lemma 4.8. Under assumption (4.6), one has J 10 < 0, J 20 < 0, J 30 < 0, T m 0 < 0 and T c 0 > 0. Furthermore, if L 1 L 2 > R 1 R 2 , then T m 1 < 0 and T c 1 < 0. Our main result now follows. Theorem 4.9. Assume (4.6), L 1 L 2 > R 1 R 2 and R 1 −L 1 R 2 −L 2 < D 2 D 1 . Then, for small ε > 0 and ν > 0, there exists a critical potential V * > 0 such that for 0 < V < V * , one has P(V) > 0, that is, the ion channel prefers cation S I over cation S II .
Remark 4.10. In Theorem 4.9 (see also Theorems 4.3 and 4.4 in Subsection 4.1), under electroneutrality boundary concentration conditions, we study the competition between two cations in terms of the function P defined in (4.1), which provides useful insights for the study of ionic flows through membrane channels, in particular, the selectivity of cations due to the electrodiffusion property of ions with given protein structures of ion channels (another important factor in studying ion selectivity). Distinct effects of the nonlinear interplays between the physical parameters, such as relative ion size (λ 1 , λ 2 ), diffusion coefficients (D 1 , D 2 ), boundary concentrations (L 1 , L 2 , L 3 ), (R 1 , R 2 , R 3 ) and boundary potential V are characterized in a great detail. For example, under the assumption of (4.6), L 1 L 2 < R 1 R 2 , and R 1 −L 1 R 2 −L 2 < D 2 D 1 , one is able to identify some critical potential V * , such that for 0 < V < V * , the ion channel prefers cation S I over cation S II . Together with the scaling law stated in Proposition 3.6 and further illustrated in Remark 3.7, the preference of ion channels for distinct cations could be either reduced or enhanced by choosing some suitable positive parameter α in Proposition 3.6.
To end this section, we perform numerical simulations for the system (3.1) directly with small ε and ν to detect two critical potentials identified in Theorem 4.3 and Theorem 4.9, respectively (Figures 2  and 3).

Concluding remarks
We studied a quasi-one-dimensional PNP model with three ion species, two positively charged with the same valence and one negatively charged, and Bikerman's local hard-sphere potential accounted for ion size effects on ionic flows. Under the framework of a geometric singular perturbation theory, the existence and uniqueness of the boundary value problem was established. Furthermore, treating the ion sizes as small parameters, both zeroth order and first order approximations (in ν) to individual fluxes are derived. Of particular interest is to examine the competition between two positively charged ion species for open ion channels with given protein structures, which provides useful insights into the study of selectivity phenomena. Based on our rigorous analysis, we are able to characterize the distinct effects of the nonlinear interplay between physical parameters, such as relative ion size (λ 1 , λ 2 ), diffusion coefficients (D 1 , D 2 ), boundary concentrations (L 1 , L 2 , L 3 ), (R 1 , R 2 , R 3 ) and boundary potential V, which provides an efficient way to control ionic flows. These are the novelty and main contribution of our work. The results, although established for simple biology settings have demonstrated extremely rich dynamics of ionic flows and sensitive dependence on all those physical parameters. More complex phenomena for more realistic ion channel models are expected. We would like to point out that the topic studied in this work is closely related to artificial nanopores as well, where modeling mainly focuses on computing device functions: ionic currents as responses to the voltage as a stimulus (see [36,38,78,89,101] for more details). We believe this work will be useful for future numerical studies, stimulate further analytical studies of ionic flows concerning the selectivity of cations, and even provide useful insights into biological experiments.