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

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 electrodiffusion of ions (mainly Na$^+$, K$^+$, Ca$^{++}$ and Cl$^{-}$) through biological membranes, establishing 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,28,29,30,31,40,41,42,43,59,60,68,93]}:

where $X\in\Omega$ with $\Omega$ 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, $\varepsilon_r(X)$ is the relative dielectric coefficient, $\varepsilon_0$ is the vacuum permittivity; $e$ is the elementary charge, $k_B$ is the Boltzmann constant, $T$ is the absolute temperature; $\Phi$ is the electric potential, for the $i$th ion species, $C_i$ is the concentration, $z_i$ is the valence (the number of charges per particle), ${\cal J}_i$ is the flux density, ${\cal D}_i(X)$ is the diffusion coefficient, and $\mu_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,24,25,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,42,43,47,56,57,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, \cdots, 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{\cal J}_k$ is the individual flux of charge of the $k$th ion species. For fixed boundary concentrations $L_k$'s and $R_k$'s, ${\cal J}_k$'s depend on $V$ only and formula (1.4) provides a dependence of the current ${\cal I}(V; \varepsilon)$ on the voltage $V$. In terms of applications, it is important to study properties of individual fluxes ${\cal J}_k$ due to the fact that most experiments (with some exceptions) can only measure the total current ${\cal 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 quasi-one-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.

2.   Problem setup

2.1.   Excess potential and a local hard sphere model

The electrochemical potential $\mu_i(X)$ for the $i$th ion species consists of the ideal component $\mu_i^{id}(X)$ and the excess component $\mu_i^{ex}(X)$: $\mu_i(X) = \mu_i^{id}(X)+\mu_i^{ex}(X),$ where

with some characteristic number density $C_0$. The classical PNP system takes into consideration of the ideal component $\mu_i^{id}(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 $\mu_i^{ex}(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 $\mu_i^{ex}(X)$ consists of two components: the (uncharged) hard-sphere component $\mu_i^{HS}(X)$ and the excess electrostatic component $\mu_i^{ES}$. While there are very successful models for $\mu_i^{HS}(X)$, modeling of $\mu_i^{ES}(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 $\mu_i^{ex}(X)$

where $\nu_j$ is the volume of the $j$th ion species. In particular, in the following analysis, we will take $\nu_n = \nu$ so that $\lambda_n = 1.$

2.2.   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 $\mu_i$, in addition to the ideal component $\mu_i^{id}$, we also include the local hard-sphere potential $\mu_i^{bik}$ in (2.2).

(A4) The relative dielectric coefficient and the diffusion coefficient are constants, that is, $\varepsilon_r(X) = \varepsilon_r$ and ${\cal D}_i(X) = {\cal 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

The boundary conditions become, for $i = 1, 2, 3$,

We first make a dimensionless rescaling following ^[39]. Set $C_0 = \max\{{\mathcal L}_i, {\mathcal R}_i:i = 1, 2\}$ and let

We would like to point out that the dimensionless parameter $\varepsilon$ defined in (2.5) as $\varepsilon = \frac{1}{l}\sqrt{\frac{\varepsilon_r\varepsilon_{0}k_BT}{e^2C_0}}$ is directly related to the ratio $\kappa_D/l,$ where $\kappa_D = \sqrt{\frac{\varepsilon_r\varepsilon_{0}k_BT}{\sum_j(z_je)^2C_j}}$ is the Debye length; in particular, $\varepsilon = \kappa_D/l$ when $z_j^2 = 1$ and $C_j = C_0$. Typically, the parameter $\varepsilon$ 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 $\varepsilon$ 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 $\frac{1}{\int_0^1h^{-1}(x)dx}$ (for example, the individual flux will be $\frac{D_k J_k}{\int_0^1h^{-1}(x)dx}$, 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 = \varepsilon\dot{\phi}$ and $\tau = x$. System (2.8) becomes

System (3.1) will be treated as a dynamical system of phase space $\mathbb{R}^{9}$ with state variables $(\phi, u, C, J, \tau)$.

For $\varepsilon > 0$, the rescaling $x = \varepsilon\xi$ of the independent variable $x$ gives rise to the fast system,

where prime denotes the derivative with respect to $\xi$.

Let $B_L$ and $B_R$ be the subsets of the phase space $\mathbb{R}^{9}$ defined by

Then the original boundary value problem is equivalent to a connecting problem, namely, finding a solution of (3.1) or (3.2) from $B_L$ to $B_R$.

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 $\varepsilon > 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.1.   Limit fast dynamics and boundary layers

By setting $\varepsilon = 0$ in (3.1), we obtain the so-called slow manifold

By setting $\varepsilon = 0$ in (3.2), we get the limiting fast system

Note that the slow manifold $\mathcal{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 $\mathcal{Z}$ is normally hyperbolic.

We would like to point out that the normal hyperbolicity of the slow manifold ${\cal 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 $\nu > 0$ small, treating (3.6) as a regular perturbation of that with $\nu = 0$, we look for solutions $\Gamma(\xi; \nu) = (\phi(\xi; \nu), u(\xi; \nu), c_1(\xi; \nu), c_2(\xi; \nu), c_3(\xi; \nu), J_1(\nu), J_2(\nu), J_3(\nu), \tau)$ of system (3.6) of the form

Lemma 3.2. For $\nu > 0$ small, one has, up to the first order in $\nu,$

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 $\nu$:

and the first order limiting fast system in $\nu$:

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({\cal Z})$ be the stable manifold of ${\cal Z}$ that consists of points approaching ${\cal Z}$ in forward time, and let $W^u({\cal Z})$ be the unstable manifold of ${\cal Z}$ that consists of points approaching ${\cal 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\cap W^s(\mathcal{Z})$ and the boundary layer at $x = 1$ must lie in $N_R = M_R\cap W^u(\mathcal{Z})$. Let $\omega(N^L)$ be the $\omega$-limit set of $N^L$, and $\alpha(N^R)$ be the $\alpha$-limit set of $N^R$, which are defined by

Proposition 3.3. Assume that $\nu\ge 0$ is small. Under electronutrality conditions (3.4), one has

(ⅰ) The stable manifold $W^s(\mathcal{Z})$ intersects $B_L$ transversally at points

and the $\omega$-limit set of $N^{L} = M^{L} \bigcap W^s(\mathcal{Z})$ is

where $J_{i}(\nu) = J_{i0}+J_{i1}\nu+o(\nu), i = 1, 2, 3$, can be arbitrary, and

(ⅱ) The unstable manifold $W^u(\mathcal{Z})$ intersects $B_R$ transversally at points

and the $\alpha$-limit set of $N^{R} = M^{R} \bigcap W^u(\mathcal{Z})$ is

where $J_{i}(\nu) = J_{i0}+J_{i1}\nu+o(\nu), \ i = 1, 2, 3$, can be arbitrary, and

Proof. The result for system (3.9), the zeroth order system in $\nu$, has been obtained in ^[6], we will not repeat it here. For system (3.10), the first order system in $\nu$, it has the following four nontrivial first integrals:

We now establish the results for $\phi_1^L, c_{11}^L, c_{21}^L, c_{31}^L$ and $u_1^l$ for system (3.10). Those for $\phi_1^R, c_{11}^R, c_{21}^R$ and $c_{31}^R$ can be established in a similar way.

Note that $\phi(0) = c_{11}(0) = c_{21}(0) = c_{31}(0) = 0.$ Using the integrals $H_1, H_2$ and $H_3,$ one has

Let $\xi\to\infty,$ careful calculations give

Under electroneutrality conditions (3.4), one has $c_{10}^L = L_1, \ c_{20}^L = L_2$ and $c_{30}^L = L_3.$ It follows that $c_{11}^L = c_{21}^L = c_{31}^L = 0.$ In view of $H_4(0) = H_4(\infty),$ together with the above analysis, one has $u_1^l = 0.$ This completes the proof.

3.2.   Limit slow dynamics and regular layer over $[0, 1]$

We now construct the regular layer $\Lambda$ on $\mathcal{Z}$ that connects $\omega(N^{L})$ and $\alpha (N^{R})$. Following the ideas in ^{[32,71,72,73]}, we make a rescaling $u = \varepsilon p$ and $c_1+c_2-c_3 = -\varepsilon q$ in system (3.1). In term of the new variables, system (3.1) becomes

where ${c_3 = c_1+c_{2}+\varepsilon q}.$

It is again a singular perturbation problem and its limiting slow system is

where ${c_3 = c_1+c_2}.$

For system (3.12), the slow manifold is

Therefore, the limiting slow system on $\mathcal{S}$ is

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 $\omega(N_L)$ and $\alpha(N_L)$ given in Proposition 3.3; in particular, for $j = 0, 1,$

For convenience, we define quantities $T_0^c, \ T_1^c, \ T_0^m, \ T_1^m, \ F_1, \ F_2$ and a function $A(x)$ as, for $k = 0, 1,$

Substituting (3.14) into system (3.13), together with (3.8), by careful calculations, we obtain the zeroth order limiting slow system in $\nu$

and the first order limiting slow system in $\nu$

For the zeroth order system (3.16) (see ^[6] for details), one has

Lemma 3.4. Under electroneutrality conditions (3.4), there is a unique solution $(\phi_0(x), c_{10}(x), c_{20}(x), J_{10}, J_{20}, J_{30}, \tau(x))$ of (3.16) such that $(\phi_0(0), c_{10}(0), c_{20}(0), \tau(0)) \; = (\phi_0^L, c_{10}^L, c_{20}^L, 0, )$ and $(\phi_0(1), c_{10}(1), c_{20}(1), \tau(1)) = (\phi_0^R, c_{10}^R, c_{20}^R, 1),$ where $\phi_0^L, c_{10}^L, c_{20}^L, \phi_0^R, c_{10}^R,$ and $c_{20}^R$ are given in Proposition 3.3. It is given by

where

Note that, for $x\in(0, 1),$ under (3.4), one has $A(x) = 1-x+\frac{R_3}{L_3}x > 0.$ For convenience, we define functions $M_k(x) = M_1(L_1, L_2, R_1, R_2;\lambda_1, \lambda_2; x)$ as, $k = 1, 2, 3, 4,$

For the first order system (3.17), we have (see Appendix A for a detailed proof)

Lemma 3.5. Under electroneutrality conditions (3.4), there is a unique solution $(\phi_1(x), c_{11}(x), c_{21}(x), J_{11}, J_{21}, J_{31}, \tau(x))$ of (3.17) such that $(\phi_1(0), c_{11}(0), c_{21}(0), \tau(0)) = (\phi_1^L, c_{11}^L, c_{21}^L, 0) \; = (0, 0, 0, 0),$ and $(\phi_1(1), c_{11}(1), c_{21}(1), \tau(1)) = (\phi_1^R, c_{11}^R, c_{21}^R, 1) = (0, 0, 0, 1),$ where $\phi_1^L, \ c_{11}^L, \ c_{21}^L, \ \phi_1^R, \ c_{11}^R$ and $c_{21}^R$ are given in Proposition 3.3. It is given by

where $J_{10}, J_{20}, J_{30}, {\cal T}^m_0, {\cal T}^c_0$ and $M_k(x), \ k = 1, 2, 3, 4$ are given in (3.15), (3.18), (3.20), and

The next result can be justified directly.

Proposition 3.6. Assume electroneutrality conditions (3.4). Fixing $V$ and treating both the zeroth order approximations ${\cal J}_{k0}$'s (in $\nu$) and first order approximations ${\cal J}_{k1}$'s (in $\nu$) as functions of boundary concentrations $L_i$ and $R_i$, one has ${\cal J}_{k0}$ is homogeneous of degree one and ${\cal J}_{k1}$ is homogeneous of degree two, that is, for any $\alpha > 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 ${\cal P}(V) = D_1J_{11}-D_2J_{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.

The slow orbit, up to $O(\nu)$,

given in Lemmas 3.4 and 3.5 connects $\omega (N_L)$ and $\alpha (N_R)$. Let $\bar{M}_L$ (resp., $\bar{M}_R$) be the forward (resp., backward) image of $\omega (N_L)$ (resp., $\alpha (N_R)$) under the slow flow (3.13) on the seven-dimensional slow manifold ${\mathcal S}$. Following the idea in ^[71,72], we have

Proposition 3.8. There exists $\nu_0 > 0$ small depending on boundary conditions so that, if $\, 0\le \nu\le \nu_0$, then, on the seven-dimensional slow manifold ${\mathcal S}$, $\bar{M}_L$ and $\bar{M}_R$ intersects transversally along the unique orbit $\Lambda(x; \nu)$ given in (3.21).

3.3.   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 $\Gamma^0$ from the point

to the point $(\phi^L, 0, c_1^L, c_2^L, c_3^L, J_1, J_2, J_3, 0) \in \omega(N_L) \subset \mathcal{Z}$ and $\Gamma^1$ from the point $(\phi^R, 0, c_1^R, c_2^R, c_3^R, J_1, J_2, J_3, 1) \in z_1(N_R) \subset \mathcal{Z}$ to the point

and a regular layer $\Lambda$ on $\mathcal{Z}$ that connects the two landing points

of the two boundary layers.

We now establish the existence of a solution of (2.8)–(2.9) near the singular orbit constructed above which is a union of two boundary layers and one regular layer $\Gamma^0 \cup \Lambda \cup \Gamma^1$. The proof follows the same line as that in ^{[6,32,71,72,73]} and the main tool used is the Exchange Lemma ^{[64,65,66,100]} of the geometric singular perturbation theory.

Theorem 3.9. Let $\Gamma^0 \cup \Lambda \cup \Gamma^1$ be the singular orbit of the connecting problem system (3.1) associated to $B_L$ and $B_R$ in system (3.3). Let $\nu_0 > 0$ be as in Proposition 3.8. Then, there exists $\varepsilon_0 > 0$ small (depending on the boundary conditions and $\nu_0$) so that, if $0\le \nu\le \nu_0$ and $0 < \varepsilon\le \varepsilon_0$, then the boundary value problem (2.8)–(2.9) has a unique smooth solution near the singular orbit $\Gamma^0 \cup \Lambda \cup \Gamma^1$.

Proof. Let $\nu_0$ be as in Proposition 3.8. For $0\leq \nu\leq \nu_0$, we define $u^l = u_0^l+u_1^l\nu, \ J_1(\nu) = J_{10}+J_{11}\nu, \ J_2(\nu) = J_{20}+J_{21}\nu$ and $J_3(\nu) = J_{30}+J_{31}\nu.$ Fix $\delta > 0$ small to be determined. Let

For $\varepsilon > 0$, let $M_L(\varepsilon, \delta)$ be the forward trace of $B_L(\delta)$ under the flow of system (3.1) or equivalently of system (3.2) and let $M_R(\varepsilon)$ be the backward trace of $B_R$. To prove the existence and uniqueness statement, it suffices to show that $M_L(\varepsilon, \delta)$ intersects $M_R(\varepsilon)$ transversally in a neighborhood of the singular orbit $\Gamma^0 \cup \Lambda \cup \Gamma^1$. The latter will be established by an application of the Exchange Lemma.

Notice that $\dim B_L(\delta)$ = 4. It is clear that the vector field of the fast system (3.2) is not tangent to $B_L(\delta)$ for $\varepsilon\ge 0$, and hence, $\dim M_L(\varepsilon, \delta)$ = 5. We next apply Exchange Lemma to track $M_L(\varepsilon, \delta)$ in the vicinity of $\Gamma^0 \bigcup \Lambda \bigcup \Gamma^1$. First of all, the transversality of the intersection $B_L(\delta) \bigcap W^s(\mathcal{Z})$ along $\Gamma^0$ in Proposition 3.3 implies the transversality of intersection $M_L(0, \delta) \bigcap W^s(\mathcal{Z})$. Secondly, we have also established that $\dim \omega(N_L) = \dim N_L-1 = 3$ in Proposition 3.3 and that the limiting slow flow is not tangent to $\omega(N_L)$ in Section 3.2. Under these conditions, the Exchange Lemma ^{[64,65,66,100]} states that there exist $\rho > 0$ and $\varepsilon_1 > 0$ so that, if $0 < \varepsilon\le \varepsilon_1$, then $M_L(\varepsilon, \delta)$ will first follow $\Gamma^0$ toward $\omega(N_L) \subset \mathcal{Z}$, then follow the trace of $\omega(N_L)$ in the vicinity of $\Lambda$ toward $\{\tau = 1\}$, leave the vicinity of $\mathcal{Z}$, and, upon exit, a portion of $M_L(\varepsilon, \delta)$ is $C^1 \; O(\varepsilon)$-close to $W^u(\omega(N_L) \times (1-\rho, 1+\rho))$ in the vicinity of $\Gamma^1$. Note that $\dim W^u(\omega(N_L) \times (1-\rho, 1+\rho)) = \dim M_L(\varepsilon, \delta) = 5$.

It remains to show that $W^u(\omega(N_L) \times (1-\rho, 1+\rho))$ intersects $M_R(\varepsilon)$ transversally since $M_L(\varepsilon, \delta)$ is $C^1 \; O(\varepsilon)$-close to $W^u(\omega(N_L) \times (1-\rho, 1+\rho))$. Recall that, for $\varepsilon = 0$, $M_R$ intersects $W^u({\cal Z})$ transversally along $N_R$ (Proposition 3.3); in particular, at $\gamma_1: = \alpha(\Gamma^1)\in \alpha(N_R)\subset {\cal Z}$, we have

where, $T_{\gamma_1}W^u(\gamma_1)$ is the tangent space of the one-dimensional unstable fiber $W^u(\gamma_1)$ at $\gamma_1$ and the vector $V_s\not\in T_{\gamma_1}W^u({\cal Z})$ (the latter follows from the transversality of the intersection of $M_R$ and $W^u({\cal Z})$). Also,

where the vector $V_{\tau}$ is the tangent vector to the $\tau$-axis as the result of the interval factor $(1-\rho, 1+\rho)$. From Proposition 3.8, $\omega(N_L)\cdot 1$ and $\alpha(N_R)$ are transversal on ${\cal Z}\cap \{\tau = 1\}$. Therefore, at $\gamma_1$, the tangent spaces $T_{\gamma_1}M_R$ and $T_{\gamma_1}W^u(\omega(N_L) \times (1-\rho, 1+\rho))$ contain seven linearly independent vectors: $V_s$, $V_{\tau}$, $T_{\gamma_1}W^u(\gamma_1)$ and the other four from $T_{\gamma_1}(\omega(N_L)\cdot 1)$ and $T_{\gamma_1}\alpha(N_R)$; that is, $M_R$ and $W^u(\omega(N_L) \times (1-\rho, 1+\rho))$ intersect transversally. We thus conclude that, there exists $0 < \varepsilon_0\le \varepsilon_1$ such that, if $0 < \varepsilon\le \varepsilon_0$, then $M_L(\varepsilon, \delta)$ intersects $M_R(\varepsilon)$ transversally.

For uniqueness, note that the transversality of the intersection $M_L(\varepsilon, \delta) \cap M_R(\varepsilon)$ implies $\dim(M_L(\varepsilon, \delta) \cap M_R(\varepsilon)) = \dim M_L(\varepsilon, \delta)+\dim M_R(\varepsilon)-9 = 1$. Thus, there exists $\delta_0 > 0$ such that, if $0 < \delta\le \delta_0$, the intersection $M_L(\varepsilon, \delta) \cap M_R(\varepsilon)$ consists of precisely one solution near the singular orbit $\Gamma^0 \cup \Lambda \cup \Gamma^1$.

4.   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 ${\cal P} = D_1J_{11}-D_2J_{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 ${\cal S}_I$ to denote the cation associated with $c_1$ and $J_1$, and ${\cal S}_{II}$ to denote the one associated with $c_2$ and $J_2$.

Recall from (3.7) that $J_{k1}$'s ($k = 1, 2$) are the main terms that contain ion size effects. From (3.15), Lemmas 3.4 and 3.5, one has

For simplicity in our following discussion, we further assume that $\lambda_2 > \lambda_1$. Define $\mathcal{P}(V) = \mathcal{P}(L_1, L_2, R_1, R_2; V; \lambda_1, \lambda_2; D_1, D_2)$ as

From Lemma 3.5, $\mathcal{P}(V)$ can be rewritten as

with $\mathcal{P}_0(V) = \mathcal{P}_0(L_1, L_2, R_1, R_2; V; \lambda_1, \lambda_2; D_1, D_2)$ being defined by

where $\mathcal{N}_1(V) = \mathcal{N}_1(L_1, L_2, R_1, R_2, V; \lambda_1, \lambda_2)$ and $\mathcal{N}_2(V) = \mathcal{N}_2(L_1, L_2, R_1, R_2, V; \lambda_1, \lambda_2)$ are defined as follows:

4.1.   Competition at singularities

For the function ${\cal 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.

Lemma 4.1. Assume that $V\neq 0$. Under electroneutrality conditions, $L_3 = R_3$, $L_3 = R_3e^{-V}$ and $L_3 = R_3e^{V}$ are three removable singularities of $\mathcal{P}(V)$.

Remark 4.2. From the definition of the function $\mathcal{P}(V)$, $L_3 = R_3, \ L_3 = R_3e^{-V}$ and $L_3 = R_3e^{V}$ are removable singularities of the first order approximations $J_{k1}, \ k = 1, 2$ in $\nu$ of the individual fluxes.

We next consider the competition of two positive charged ion species as $L_3\to R_3$ and $L_3\to R_3e^{-V}$, respectively.

Theorem 4.3. Suppose $L_1 > R_1.$ For small $\varepsilon > 0$ and $\nu > 0$, there exists a critical potential $V^* > \ln 2$ such that, for $V > V^*$, one has $\mathcal{P}(V) > 0$ as $L_3\to R_3,$ that is, the ion channel prefers cation ${\cal S}_I$ over cation ${\cal S}_{II}$.

Proof. Direct calculation gives

Under our assumption, one has $\frac{(\lambda_2-\lambda_1)(L_1-R_1)V}{e^{-V}-1} < 0, \ \mbox{for}\ V > 0.$ Upon introducing $K_1(V)$ and $K_2(V)$ defined as

Then, as $L_3\to R_3,$ $\mathcal{P}(V) > 0$ is equivalent to

Rewrite $K_1(V)$ as $K_1(V) = \frac{e^{-V}-1}{2V}\tilde{K}_1(V)$ with $\tilde{K}_1(V) = V+2e^{-V}-2.$ A direct computation gives $\tilde{K}_1(0) = 0$ and $\tilde{K}_1(V) > 0$ for $V > \ln 2$; and $\tilde{K}_1(V) < 0$ for $V < \ln 2.$ Together with $\frac{e^{-V}-1}{2V} < 0$, one has $K_1(V) < 0$ for $V > \ln 2$ and $K_1(V) > 0$ for $V < \ln 2.$ Note that $K_2(V) = K_1(V)+\frac{V}{e^{V}-1}$ with $\frac{V}{e^{V}-1} > 0.$ Therefore, there must exist a $V^* > \ln 2$ such that $K_2(V) < 0$ for $V > V^*$. With $L_1 > R_1$ and $V > 0$, we have

(ⅰ) $R_1-L_1e^{V} = R_1-L_1+L_1(1-e^{V}) < 0,$ which implies $\frac{(D_1+D_2)(R_1-L_1e^{V})}{1-e^{V}} > 0.$

(ⅱ) $R_2-R_1 < L_2-L_1 < (L_2-L_1)e^{V}$ gives $\frac{D_2V}{(1-e^{V})^2}\left(R_2-R_1-(L_2-L_1)e^{V}\right) < 0.$

(ⅲ) $\frac{(D_1+D_2)(L_1-R_1)}{1-e^{V}} < 0.$

Therefore, with $L_1 > R_1,$ one has, for $V > V^*,$

This completes the proof.

Theorem 4.4. Suppose $\frac{D_1}{D_2} > \frac{L_2}{L_1}$ and $V > 0$ (resp. $\frac{D_1}{D_2} < \frac{L_2}{L_1}$ and $V < 0$). For small $\varepsilon > 0$ and $\nu > 0$, one has $\mathcal{P}(V) > 0$ as $L_3\to R_3e^{-V},$ that is, the ion channel prefers cation ${\cal S}_I$ over cation ${\cal S}_{II}$ as $L_3\to R_3e^{-V}.$

Proof. It follows from ${\cal P}(V)\to\frac{(D_1L_1-D_2L_2)(e^{V}-1)^2(\lambda_1L_1+\lambda_2L_2+R_3e^{-V})}{V}$ as $L_3\to R_3e^{-V}.$

4.2.   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_3e^{V}$ and $L_3 = R_3e^{-V}$ split the $L_3$-region into four subregions: $\left(0, R_3e^{-V}\right), \ \left(R_3e^{-V}, R_3\right), \ \left(R_3, R_3e^{V}\right)\ $ and $\ \left(R_3e^{V}, +\infty\right),$ over which the sign of $\mathcal{P}$ can be analyzed. Our study will mainly focus on the sign of ${\cal P}$ over the first subregion $\left(0, R_3e^{-V}\right).$

To get started, we establish some results which are crucial to study the sign of ${\cal P}(V)$, and will be frequently used in our following analysis. Note that, under electroneutrality conditions (3.4), from (3.18), one has

and from (3.15), (5.3) and (5.4), together with $A(1) = R_3/L_3,$ we have

where

It then follows that

Lemma 4.5. Under electroneutrality conditions, one has

(ⅰ) For $L_3 < R_3e^{-V}$, one has $F_1 < 0$ (resp. $F_1 > 0$) if $\frac{L_1}{L_2} < \frac{R_1}{R_2}$ (resp. $\frac{L_1}{L_2} > \frac{R_1}{R_2}$).

(ⅱ) For $L_3 > R_3e^{-V},$ one has $F_1 < 0$ (resp. $F_1 > 0$) if $\frac{L_1}{L_2} > \frac{R_1}{R_2}$ (resp. $\frac{L_1}{L_2} < \frac{R_1}{R_2}$).

Lemma 4.6. Fixing $V > 0.$ Then,

(Ⅰ) For $R_3 < 1,$ one has

(ⅰ) $q_{1}(V; L_3, R_3) > 0$ if $L_3 < R_3;$

(ⅱ) There exists an $L_3^{*} > R_3$ such that $q_{1}(V; L_3, R_3) < 0$ if $R_3 < L_3 < L_3^{*}$ and $q_1(V; L_3, R_3) > 0$ if $L_3 > L_3^{*}.$

(Ⅱ) For $R_3 > 1,$ one has $q_{1}(V; L_3, R_3) > 0\ \mathit{\mbox{if}}\ L_3 < R_3;\ q_{1}(V; L_3, R_3) < 0\ \mathit{\mbox{if}}\ L_3 > R_3.$

Lemma 4.7. $q_2(L_3, R_3) > 0$ if $L_3 > R_3;$ and $q_2(L_3, R_3) < 0$ if $L_3 < R_3.$

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_3e^{-V}: \; (a)\ L_1 < R_1e^{-V}$ and $L_2 < R_2e^{-V};$ $(b)\ L_1 < R_1e^{-V}$ and $L_2 > R_2e^{-V};$ and $(c)\ L_1 > R_1e^{-V}$ and $L_2 < R_2e^{-V}.$ In this section, we will only consider the case (a), which, with $V > 0$, is equivalent to

and study the sign of ${\cal 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, \ \ {\cal T}_0^m < 0$ and ${\cal T}_0^c > 0.$ Furthermore, if $\frac{L_1}{L_2} > \frac{R_1}{R_2}$, then ${\cal T}_1^m < 0$ and ${\cal T}_1^c < 0.$

Our main result now follows.

Theorem 4.9. Assume (4.6), $\frac{L_1}{L_2} > \frac{R_1}{R_2}$ and $\frac{R_1-L_1}{R_2-L_2} < \frac{D_2}{D_1}$. Then, for small $\varepsilon > 0$ and $\nu > 0$, there exists a critical potential $V_* > 0$ such that for $0 < V < V_*,$ one has $\mathcal{P}(V) > 0,$ that is, the ion channel prefers cation ${\cal S}_I$ over cation ${\cal S}_{II}$.

Proof. Note that, under our assumption, $\mathcal{P}\mathcal{P}_0 > 0.$ Instead of showing that $\mathcal{P} > 0,$ we will prove that $\mathcal{P}_0 > 0.$ Recall from (4.2),

where ${\cal N}_1$ and ${\cal N}_2$ as defined in (4.3). From Lemma 4.8, one can easily check that for the second term in ${\cal P}_0$, one has $-\left(D_1J_{10}+D_2J_{20}\right)\frac{(\lambda_2-\lambda_1)F_1{\cal T}^c_0L_3}{{\cal T}^m_0}(A(1)-1) > 0.$ Next we show that, under our assumption, $(D_2J_{20}-D_1J_{10})\mathcal{N}_1+(D_1+D_2)(J_{10}+J_{20})\mathcal{N}_2 > 0$ for $0 < V < V_*,$ where $V_*$ is identified at the end of our proof. For convenience, we rewrite ${\cal N}_1$ and ${\cal N}_2$ as

where

Notice that, with $L_1 < R_1e^{-V}$ and $L_2 < R_2e^{-V},$

Together with Lemma 4.8, one has ${\cal N}_{10} > 0\ \ \mbox{and}\ \ {\cal N}_{20} > 0.$

Finally we prove that

From (3.15), (3.18) and (5.3),

which is true under our assumption. It suffices to show that

which is, from Lemma 3.4 and Lemma 4.8, equivalent to $\omega(V) < 0$, where $\omega(V): = \omega(V; L_1, L_2, L_3, R_1, R_2, R_3; D_1, D_2)$ is defined as

With $\frac{R_1-L_1}{R_2-L_2} < \frac{D_2}{D_1},$ a careful calculation leads to $\omega(0) = 0;\ \omega'(0) < 0$ and $\omega''(V) > 0$ for $V > 0.$ It follows that the curve $\omega(V)$ is concave up for $V > 0$ with its minimum being negative. Therefore, there exists a unique $V_g$ such that $\omega(V) < 0$ for $0 < V < V_g$. Define $V_* = \min\left\{V_g, \ln\frac{R_1}{L_1}, \ln\frac{R_2}{L_2}\right\}$, one has $\mathcal{P}(V) > 0$ for $0 < V < V_*$.

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 ${\cal 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 $(\lambda_1, \lambda_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), $\frac{L_1}{L_2} < \frac{R_1}{R_2}$, and $\frac{R_1-L_1}{R_2-L_2} < \frac{D_2}{D_1}$, one is able to identify some critical potential $V_*$, such that for $0 < V < V_*,$ the ion channel prefers cation ${\cal S}_I$ over cation ${\cal 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 $\alpha$ in Proposition 3.6.

To end this section, we perform numerical simulations for the system (3.1) directly with small $\varepsilon$ and $\nu$ to detect two critical potentials identified in Theorem 4.3 and Theorem 4.9, respectively (Figures 2 and 3).

5.   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 $\nu$) 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 $(\lambda_1, \lambda_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.

Appendix A. We provide detailed proof for Lemma 3.5 in section 3.2.

Proof of Lemma 3.5 From $c_{11}$ and $c_{21}$ equations in (3.17), one has

Together with Proposition 3.3, one has

From (3.18), together with $x = \frac{2L_3}{{\cal T}^m_0}(1-A(x)),$ we have

Therefore,

Evaluating Eq (5.2) at $x = 1$, we have

It follows that

Solving the $\phi_{1}$ equations in (3.17), we have

Note that, from (3.18) and (5.2), together with $1-A(x) = \frac{{\cal T}^m_0}{2L_3}x,$ one has

Therefore,

Evaluating $\phi_1(x)$ at $x = 1$ gives

By the variation of constant formula, one has, with $c_{11}^L = c_{21}^L = 0,$

where, from (3.18), (5.1) and (5.2),

Evaluating $\phi_1(x), \ c_{11}(x)$ and $c_{21}(x)$ at $x = 1,$ the expressions for $J_{11}, J_{21}$ and $J_{31}$ follows directly.

Acknowledgments

The authors are grateful to the anonymous referees whose suggestions have in our opinion, significantly improved the paper. J. Chen and M. Zhang were partially supported by MPS Simons Foundation (No. 628308).

Conflict of interest

The author declares no conflicts of interest in this paper.