ION SIZE EFFECTS ON INDIVIDUAL FLUXES VIA POISSON-NERNST-PLANCK SYSTEMS WITH BIKERMAN’S LOCAL HARD-SPHERE POTENTIAL: ANALYSIS WITHOUT ELECTRONEUTRALITY BOUNDARY CONDITIONS

A quasi-one-dimensional steady-state Poisson-Nernst-Planck model with Bikerman’s local hard-sphere potential for ionic flows of two oppositely charged ion species through a membrane channel is analyzed. Of particular interest is the qualitative properties of ionic flows in terms of individual fluxes without the assumption of electroneutrality conditions, which is more realistic to study ionic flow properties of interest. This is the novelty of this work. Our result shows that i) boundary concentrations and relative size of ion species play critical roles in characterizing ion size effects on individual fluxes; ii) the first order approximation Jk1 = DkJk1 in ion volume of individual fluxes Jk = DkJk is linear in boundary potential, furthermore, the signs of ∂V Jk1 and ∂2 V λJk1, which play key roles in characterizing ion size effects on ionic flows can be both negative depending further on boundary concentrations while they are always positive and independent of boundary concentrations under electroneutrality conditions (see Corollaries 3.2-3.3, Theorems 3.4-3.5 and Proposition 3.7). Numerical simulations are performed to identify some critical potentials defined in (2). We believe our results will provide useful insights for numerical and even experimental studies of ionic flows through membrane channels.

1. 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 electro-diffusion of ions (mainly Na + , K + , Ca ++ and Cl − ) through biological membranes, establishing communications among cells and the external environment.In this way, ion channels λ = ν 2 /ν and fixed boundary concentrations L k and R k as follows (will be further discussed in Section 2): This is crucial to further study the qualitative properties of ionic flows.It turns out that, for k = 1, 2, both J k1 (V ; λ, ε) and I 1 (V ; λ, ε) are linear functions in V (see (12) in Section 2.2).In particular, six critical potentials V kc , V c k , V c , V c are defined by J k1 (V kc ; λ, 0) =0, I 1 (V c ; λ, 0) = 0, d dλ J k1 (V c k ; λ, 0) = 0, ( The significance of the six critical potentials is apparent from their definitions (See Theorems 4.8, 4.9, 4.18 and 4.19 in [41] for details), and from (1), the sign of ∂ V J k1 (resp.∂ V λ J k1 , ∂ V I 1 and ∂ V λ I,) plays a key role in characterizing the effect from ion sizes.Furthermore, the signs of those terms depend sensitively on multiple physical parameters, such as the boundary potentials and concentrations, ion valences, diffusion coefficients and ion sizes.The characterization of the distinct effects of the nonlinear interplay between these physical parameters will provide detailed information for one to better understand the flow property of interest.
Two interesting questions arising immediately are i) in general, when ∂ V I 1 (V ; λ, 0) and ∂ V λ I 1 (V ; λ, 0) are negative?ii) Can the total flux be studied in terms of the individual fluxes J k = D k J k ?
To answer these questions, in this work, we study a quasi-one-dimensional PNP model with the same setting as that in [41] except for the assumption of electroneutrality conditions.Of particular interest is the sign study of ∂ V J 11 (V ; λ, 0) and ∂ V λ J 11 (V ; λ, 0) in terms of the individual flux.We take particular advantage of the work in [41] to provide a detailed explanation of how these physical parameters interact to produce a wide spectrum of behaviors for ionic flows.Our rigorous analysis (see Theorems 3.2-3.3,Corollaries 3.4-3.5 and Proposition 3.7) shows that (i) boundary concentrations and relative size of ion species play critical roles in characterizing ion size effects on individual fluxes; (ii) the first order approximation J k1 = D k J k1 in ion volume of individual fluxes J k = D k J k is linear in boundary potential, furthermore, the signs of ∂ V J k1 and ∂ 2 V λ J k1 which play key roles in characterizing ion size effects on ionic flows can be both negative depending further on boundary concentrations while they are always positive and independent of boundary concentrations under electroneutrality conditions.This is our main contribution in this work.To the best of the authors' knowledge, this work is the first analysis on roles that ion size plays in individual fluxes without electroneutrality conditions.
We emphasize 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.
It should be pointed out that the quasi-one-dimensional PNP model and the local hard-sphere model (see (7) below) adopted in [41] and in this paper are rather simple.Aside the trivial fact that they will miss the three-dimensional features of the problem, a major weakness is the missing of the excess electrostatic component in the excess potentials.Important phenomena such as charge inversion and layering may not be detected by this simple model.
The rest of the paper is organized as follows.In Section 2, we describe the quasi-one-dimensional PNP model of ion flows, a local model for hard-sphere (HS) potentials, the formulation of the boundary value problem of the singularly perturbed PNP-HS system, and the basic assumptions.Some results from [48] are recalled, and these will be the starting point of our study.In Section 3, without the assumption of electroneutrality conditions, we study the sign of ∂ V J k1 (resp.∂ V λ J k1 , ∂ V I 1 and ∂ V λ I), which plays a key role in characterizing the effect from ion sizes.It also turns out that the signs of those terms depend sensitively on multiple physical parameters such as boundary concentrations, ion valence and ion sizes.Partial orders of critical potentials identified in (2) are provided.Numerical simulations are performed to system ( 9)- (10) to identify some critical potentials defined in (2).The paper ends with a concluding remark provided in Section 4.
2. Models and some previous results.In this section, we briefly recall the model of Poisson-Nernst-Planck (PNP) systems, and some main results obtained in [41] related to ion size effects on ionic flows.

2.1.
A quasi-one-dimensional PNP system.We assume the channel is narrow so that it can be effectively viewed as a one-dimensional channel that connects the interior and the exterior of the channel.A quasi-one-dimensional steady-state PNP model for ion flows of n ion species though a single channel is (see [50,52]) where X ∈ [0, l], e is the elementary charge, k B is the Boltzmann constant, T is the absolute temperature; Φ is the electric potential, Q(X) is the permanent charge of the channel, ε r (X) is the relative dielectric coefficient, ε 0 is the vacuum permittivity; A(X) is the area of the cross-section of the channel over the point X ∈ [0, l]; for the ith ion species, C i is the concentration (number of ith ions per volume), z i is the valence (number of charges per particle) that is positive for cations and negative for anions, µ i is the electrochemical potential, J i is the flux density, and D i (X) is the diffusion coefficient.
For system (3), we impose the following boundary conditions (see, [21] for justification), for k = 1, 2, • • • , n, For ion channels, an important characteristic is the so-called I-V relations (currentvoltage relations).For a solution of the steady-state boundary value problem (3)-( 4), the rate of flow of charge through a cross-section or current I is For fixed boundary concentrations L i 's and R i 's, J j 's depend on V only and formula (5) provides a relation of the current I on the voltage V .This relation is the I-V relation.
2.1.1.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 µ ex i (X): 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 classical 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 the finite size effect of charges (see, e.g., [60,61]).
In this paper, we will take Bikerman's local hard-sphere model for µ ex i (X) where ν j is the volume of a single jth ion species.We would like to point out that since C j is the number density of ith ion species, it follows that n j=1 ν j C j < 1.In this sense, Bikerman's LHS takes into consideration of nonzero ion sizes, however, it is not ion specific, more precisely, one has µ Bik 2.1.2.The steady-state boundary value problem and assumptions.The main goal of this paper is to examine the qualitative properties of the ion size effect on ionic flows via the steady-state PNP system (3)-(4).For definiteness, we will take essentially the same setting as that in [39], that is, (A1).We consider two ion species (n = 2) with z 1 > 0 and z 2 < 0. (A2).We assume the permanent charge Q(X) to be zero.(A3).For the electrochemical potential µ i , in addition to the ideal component µ id i (X) defined in (6), we also include a local hard-sphere potential (7) to approximate the excess component µ ex i (X).(A4).The relative dielectric coefficient and the diffusion coefficient are constants, that is, ε r (X) = ε r and D i (X) = D i .In the sequel, we will assume (A1)-(A4).We first make a dimensionless rescaling following ( [26]).Set C 0 = max{L i , R i : i = 1, 2} and let Using expressions (6) for the ideal component µ id i (X) and ( 7) for the excess chemical potential µ ex i (X), the boundary value problem ( 3)-( 4) becomes with boundary conditions, for i = 1, 2, To end this section, we would like to point out that (i) for typical ion channel problems, physical range for the parameter ε is 10 −2 − 10 −6 , which is smaller for crowded ionic mixtures (large C 0 ) and larger for less crowded ionic mixtures.It is further assumed that the dimensionless parameters ν i 's are small; typical physical range for with 10 −2 corresponding to crowded ionic mixtures, say, C 0 ∼ 10 M (molar) and with 10 −4 to less crowded ionic mixtures, say, C 0 ∼ 100 mM.. (ii) we 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 (for example, the individual flux will be , see [48]), which does not affect our analysis of the qualitative properties of the ionic flows.
2.2.Some previous results.We now recall some results obtained in [41], which are crucial for our study and which will be frequently used.In [41], the authors treat ε and ν = ν 1 as small parameters and derive approximations for the current I and the individual flux J k expanded in ν with λ = ν 2 /ν for k = 1, 2, where is the zeroth order expansion in ν, J k1 is the first order expansion in ν, and In particular, Here, where c L 10 and c L 20 are landing points given by Under our setups, the critical potentials V kc , V c k , V c , and V c identified in (2) are given explicitly by From ( 14) and (15), we observed that k , V c and V c are homogeneous of degree zero in (L 1 , R 1 ), that is, for any s > 0, Note that, from ( 13) and ( 14), one has To better answer the question imposed in the introduction, the first step is to study the the signs of ∂λ (V ; λ, 0) in terms of the individual flux without assuming the electroneutrality conditions.
3. Qualitative properties of ionic flows without electroneutrality conditions.In this section, we would like to study the ion size effects on ionic flows without electroneutrality conditions.More precisely, we would like to study the sign of , which plays a key role in characterizing the effects on ionic flows from ion sizes.
Due to the difficulty in analysis, we will consider a relatively simple but more general (compared to the study in [41]) case with for some positive constants σ, which is not equal to 1 (σ = 1 implies electroneutrality conditions on boundary concentrations).Under this assumption, one has c L 10 = σ z 1 It follows form ( 12) and ( 14) that where where To get started, we establish the following result, which will be used frequently in our analysis.
Proof.For convenience, for For n(x), one has n (x) = − (x−1) 2 Hence, m (x) > 0 if 0 < x < 1; and m (x) < 0 if x > 1.Note that m(x) → 0 as x → 0; and m(x) → 0 as x → ∞.Therefore, m(x) has its global maximum at x = 1 (that is, as , which is 1 2 .This completes the proof. 2 plays an important role in characterizing the size effects on ionic flows.It is worth analyzing it in details.Note that, from ( 12) and ( 14), one has Therefore, we will mainly focus on the study of the individual flux J 11 .We first consider the sign of ∂ V J 11 .
Proof.Note that the sign of ∂ V J 11 (V ; σ; λ, 0) is determined by the sign of f (σ) since e k B T > 0 and σ From equation (19), direct calculations give It is easy to see that f (σ) > 0 for σ > 0, which implies that f (σ) is concave up for σ > 0. It follows that f (σ) has a global minimum at some point σ c , which is the unique solution of f (σ) = 0.It is given by It is clear that f (σ) obtains its global minimum at σ = σ c , which is given by We now consider the sign of f (σ c ).For convenience, we redefine the quantity f (σ c ) as a function of m, denoted by g(m).Note that z 1 ≥ 1 and Note also that One can easily check that g(0+) > 0 if λ * 1 < λ < λ * 2 , and g(0+) < 0 if either 0 < λ < λ * 1 or λ > λ * 2 where λ * 1 and λ * 2 are two distinct real roots of We next study the sign of g(m) for two cases.
A similar argument leads to the result related to the sign of ∂ V λ J 11 .
is the unique real root of f (σ) = 0.
To end this section, we would like to comment that since α 11 = −β 11 in ( 14), one has ∂ Corresponding results related to J 21 can be easily obtained, and we leave them to the readers.

3.2.
Sign studies related to the total flux.Some qualitative properties of the total flux (the total flow rate of charge) in terms of [41].In particular, under the assumption of electroneutrality conditions, the author established that ∂ V I 1 (V ; λ, 0) > 0 (resp.∂ 2 V λ I 1 (V ; λ, 0) > 0 ), which implies that, for small ε > 0 and However, there are not too much discussion about the negativeness of these terms except a very special case stated in our introduction.
From (16), the study of the signs of ∂ 2 V λ I 1 (V ; λ, 0) and ∂ V I 1 (V ; λ, 0) follows directly from our analysis of the individual fluxes (Theorem 3.8 and Theorem 3.14).More precisely, Proposition 3.9.Under assumption (17), one has 3.3.Partial orders of critical potentials.Our main interest in this section is to provide a partial order for the critical potentials identified in (2) (see [41] for more details).We first consider the critical potentials V 1c , V 2c and V c that balance the ion size effects on both the individual flux and the total flux.For the six critical potentials, V kc , V c k related to the individual fluxes, and V c , V c related to the total flux, given by (15) (see also [41]), it is expected that V c and Vc depend on V 1c and V 2c , and V c and V c depend on V c 1 and V c 2 .More precisely, Corollary 3.10.We have Proof.It follows from the relation ( 14) and the expressions of the six critical potentials in (15) directly.We omit it here.
The following result can be established.
, there exist two distinct real roots σ1 and σ2 with σ1 < σ2 of p(σ) Proof.We only provide the proof for the case with α 11 > 0. Similar argument can be easily applied to the case with α 11 < 0.
From (15) and Corollary 3.10, one has where
Proposition 3.14.For fixed It follows that, for fixed (λ * , ν * ), one can numerically compute J 1 (V ; λ, ν * ), and hence H(V, λ) for any λ close to λ * , then apply Proposition 3.14 to estimate V c * 1 from the saddle point of H(V, λ).In our numerical simulations, we fix λ * = λ/3 and plot H(V, λ) as a function of V with three different λ values (see the second column in Figure 4 for J 1 , J 2 and I).Our analytical results for zeroth order in ε and first order in ν tell us that the graphs for these three different λ values should have a common intersection point with V = V c 1 .4. Concluding Remarks.In this work, we consider the PNP model with local excess chemical potentials to account for finite size effects on ionic flows for two ion species, one positively charged and one negatively charged.The qualitative properties of ionic flows, in terms of individual fluxes, are studied without the assumption of electroneutrality conditions, which is more realistic, and it turns out that these properties depend very sensitively on boundary concentrations (in terms of m(L 1 , R 1 )) and relative ion sizes (in terms of λ) in addition to other physical parameters such as boundary potentials, ion valences and ion sizes.Our result shows that the first order approximation J k1 = D k J k1 in ion volume of individual fluxes J k = D k J k is linear in boundary potential, and the signs of ∂ V J k1 and ∂ 2 V λ J k1 which play key roles in characterizing ion size effects on ionic flows can be both negative depending further on boundary concentrations while they are always positive and independent of boundary concentrations under electroneutrality conditions.This is the novelty of our work.For the relatively simple setting and assumptions of the model studied in this paper, we are able to characterize the distinct effects of the nonlinear interplay between these physical parameters.The dynamical behaviors of ionic flows studied in this work are much more rich compared to those in [41] under electroneutrality conditions.
Finally, we comment that the model studied in this paper is oversimplified and the specific setting (such as zero permanent charge assumption) of our problem may not reflect precisely any realistic biological settings.However, we believe the existence of these critical potentials are generally in valid (even for PNP model with nonzero permanent charge, see Figures 5 and 6 for V 1c and V c 1 with nonzero permanent charge Q(x) defined piecewise by Q(x) = 0 if x ∈ [0, 1/3] ∪ [2/3, 1] and Q(x) = Q 0 if x ∈ [1/3, 2/3], where Q 0 = 0.5) and the awareness of the potential existence of these critical voltage itself would be useful for further numerical studies and stimulate further analytical studies of ionic flows through ion channels.We believe our results will provide useful insights for numerical and even experimental studies of ionic flows through membrane channels.

Figure 1 .Figure 2 .
Figure 1.Numerical detection of critical values for λ (left graph) and m (right one) in Theorem 3.2.

Figure 5 .
Figure 5. Numerical identification of critical potentials V1c and V c 1 for individual flux J1 with z1 = −z2 = 1 and nonzero permanent charge.The x-axis for all figures actually represents e k B T V .

Figure 6 .
Figure 6.Numerical approximations of critical potentials V c 1 for individual flux J1 with z1 = −z2 = 1 and nonzero permanent charge as illustrated in Proposition 3.14.The x-axis for all figures actually represents e k B T V .