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

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

1. Introduction

Ion channels are large proteins embedded in cell membranes with a hole down their middle that provides a controllable path for electrodiffusion of ions (mainly Na, K, Ca and Cl) through biological membranes, establishing communications among cells and the external environment [1,2,3]. In general, the study of ion channels consists of two related major topics: structures of ion channels and ionic flow properties.

The physical structure of ion channels is defined by the channel shape and the spacial distribution of permanent and polarization charge. Very often, the shape of a typical ion channel is approximated by a cylindrical-like domain with a non-uniform cross-section area. Within a large class of ion channels, amino acid side chains are distributed mainly over a relatively “short” and “narrow” portion of the channel, where acidic side chains contribute permanent negative charges and basic side chains contributes permanent positive charges, and this is analogous to the doping of semiconductor devices, e.g., bipolar PNP and NPN transistors [1,4].

With a given structure of an open channel, the main interest is to understand its electrodiffusion property. Mathematical analysis plays important and unique roles for generalizing and understanding the principles that allow control of electrodiffusion, explaining mechanics of observed biological phenomena and for discovering new ones, assuming a more or less explicit solution of the associated mathematical model can be obtained. However, in general, the latter is too much to expect. Recently, there have been some successes in mathematical analysis of Poisson–Nernst–Planck (PNP) models for ionic flows through membrane channels [5,6,7,8,9,10,11,12,13,14].

1.1. Poisson–Nernst–Planck Models for Ionic Flows

Considering the structural characteristics, the basic continuum model for ionic flows is the Poisson–Nernst–Planck system, which treats the aqueous medium as a dielectric continuum ([15,16,17,18] etc.). The PNP system can be derived as a reduced model from molecular dynamics [19], from Boltzmann equations [20], and from variational principles [21,22,23].

It is known that more sophisticated models [24,25,26,27,28,29,30,31] have also been developed which can model the physical problem more accurately, however, it is very challenging to examine their dynamics analytically and even computationally. Considering the key feature of the biological system, the PNP system represents an appropriate model for both analysis and numerical simulations of ionic flows.

The simplest PNP system is the classical Poisson–Nernst–Planck (cPNP) system that includes the ideal component in (4) only. The ideal component contains contributions by considering ion particles as point charges and ignoring the ion-to-ion interaction. It has been shown by some numerical studies that classical PNP models provide good qualitative agreement with experimental data for I-V relations [20,32]. The classical PNP models have been simulated and analyzed extensively (see, e.g., [7,12,13,14,20,32,33,34,35,36,37,38,39]).

For ionic solutions with n ion species, the PNP system reads where with being a three-dimensional cylindrical-like domain representing the channel, is the permanent charge density, is the relative dielectric coefficient, is the vacuum permittivity, e is the elementary charge, is the Boltzmann constant, T is the absolute temperature; is the electric potential. Furthermore, for the kth ion species, is the concentration, is the valence, is the electrochemical potential depending on and , is the flux density, and is the diffusion coefficient.

Based on the fact that ion channels have narrow cross-sections relative to their lengths, reduction of the three-dimensional steady-state PNP systems (1) to quasi-one-dimensional models was first proposed in [40] and was rigorously justified in [36] for special cases. A quasi-one-dimensional steady-state PNP model takes the form where is the coordinate along the axis of the channel, is the area of cross-section of the channel over the location X.

Equipped with system (2), we impose the following boundary conditions (see, Ref. [7] for a reasoning), for ,

1.2. Electrochemical Potential

The electrochemical potential for the ith ion species consists of the ideal component and the excess component : , where with some characteristic number density defined by

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

1.3. Permanent Charge

The spatial distribution of side chains in a specific channel defines the permanent charge of the channel. While some information may be obtained by ignoring the permanent charge and focusing on the effects of boundary conditions, the valences and sizes of ions, etc., we believe that different channel types differ mainly in the distribution of permanent charge [3]. To better understand the importance of the relation of ionic flows and permanent charges, we remark that the role of permanent charges in membrane channels is similar to the role of doping profiles in semiconductor devices. Semiconductor devices are similar to membrane channels in the way that they both use atomic-scale structures to control macroscopic flows from one reservoir to another. Ions move a lot like quasi-particles move in semiconductors. Roughly, holes and electrons are the cations and anions of semiconductors. Semiconductor technology depends on the control of migration and diffusion of quasi-particles of charge in transistors and integrated circuits. Doping is the process of adding impurities into intrinsic semiconductors to modulate its electrical, optical, and structural properties [45,46]. One may roughly understand in the following sense, doping provides the charges that acid and basic side chains provide in a protein channel. For both ion channels and semiconductors, permanent charges add an additional component—probably the most important one—to their rich behavior.

In general, the permanent charge is modeled by a piecewise constant function, that is, we assume, for a partition of into m subintervals, for where ’s are constants with (the intervals and are viewed as the reservoirs where there is no permanent charge). In this work, we take the following model for where and is some nonzero constant.

1.4. Comparison with Some Existing Works

The current work follows a similar dynamical system framework as that employed in [6,7,9,10,11] to establish the existence and uniqueness result of the problem. However, compared to these works, our set-ups are much more challenging and more realistic, more importantly, the specific structure of our model allows us to obtain detailed description of the nonlinear interplay among different system parameters. This is far beyond the existence and uniqueness result. To be specific, our model includes three ion species, two positively charged with the same valences, and one negatively charged (in [7,9,10], only two oppositely charged particles are included, selectivity of cations, one of the most relevant biological properties of ion channels cannot be described); and a profile of nonzero but small permanent charges (in [6], it includes three ion species but with zero permanent charges, the effects on ionic flows from the two key structures of ion channels: channel geometry and distribution of permanent charges, cannot be examined, while this could provide crucial insights for the selectivity of cations through membrane channels). In [11], the author extended the work in [7] and established the existence and local uniqueness of the classical PNP system with n ion species.

The work, in some sense, is motivated by [9], and there are some similarities in the treatment. More precisely, both of the works employ regular perturbation analysis to derive the explicit expressions of the individual fluxes up to the first order in the small permanent charge, which is reflected in Section 2.2 in the current work. However, the derivation and the following analysis is much more challenging due to the nonlinearity of the individual fluxes in the potential V (in [9], the individual fluxes are linear in the potential V). The nonlinearity of the individual fluxes in the potential provides much more rich dynamics of ionic flows, and demonstates more complicated nonlinear interaction among the system parameters, which is addressed in Section 3.1. Meanwhile, this indicates that our work provides a better understanding of the mechanism of ionic flows through single ion channels, which is necessary and important for future studies of ion channel problems.

1.5. Main Results

For convenience, we briefly summarize our main results as follows with .

Constructing a singular orbit of the limiting PNP system () over the whole interval , which is a union of singular orbits over the subintervals and . Over each subinterval, the singular orbit consists of two boundary/internal layers and a regular layer; see Proposition 2 and Lemma 2 in Section 2.1.1 for the singular over , in Proposition 4 and Lemma 3 for the one over in Section 2.1.2, and Proposition 5 and Lemma 4 for the one over in Section 2.1.3.

Establishing the existence and local uniqueness result of the underlying PNP system ( but small); see Theorem 1 in Section 2.1.5.

Obtaining the zeroth order and first order (in ) solutions of system (40) and (41), crucial to derive explicit expressions of the individual fluxes up to the first order in ; see Propositions 6 and 7 in Section 2.2.

The sign of A and , critical for our analysis in Section 3.1; see Lemmas 7 and 8 in Section 2.2.

Study on competition between cations in terms of from two directions: the sign of and the monotonicity of in the electric potential V, based on distinct interplays among and consisting of three cases

Case study with ; see Theorems 2–4 in Section 3.1.1.

Case study with ; see Theorems 5–7 in Section 3.1.2.

Case study with ; see Theorems 8 and 9 in Section 3.1.3.

Analysis on the magnitude of , equivalent to examine the sign of , where ; see Theorems 10 and 11 in Section 3.1.4.

In (v), there are actually another three cases: (1)

1.6. Problem Set-Up

For definiteness, we will take the following setting in this work:

We consider three ion species () with and .

The permanent charge is defined as in (6).

For the electrochemical potential , we only include the ideal component as in (4).

We assume and .

In the sequel, we will assume (A1)–(A4). We first make a dimensionless rescaling following [3]. With given in (5), let

The BVP (2) and (3) then becomes (noting that ) with the boundary conditions, for

We comment that the dimensionless parameter defined in (1.6) as is directly related to the ratio where is the Debye length; in particular, when and . Typically, the parameter is small due to the fact that the two variables l, the length of the channel, and , some characteristic number density could be very large. For many cases, the value of is of order (see [47] for a more detailed description).

2. Methods

2.1. Geometric Singular Perturbation Theory for (8) and (9)

We first rewrite system (8) into a standard form for singularly perturbed systems and convert the boundary value problem (8) and (9) to a connection problem.

Upon introducing and . System (8) becomes where overdot denotes the derivative with respect to the variable x.

System (10) will be treated as a singularly perturbed system with as the singular parameter. Its phase space is with state variables .

For , the rescaling of the independent variable x gives rise to where prime denotes the derivative with respect to the variable .

We comment that for , systems (10) and (11) have exactly the same phase portrait. However, their limiting systems at are different. The limiting system of (10) is called the limiting slow system, whose orbits are called slow orbits or regular layers. The limiting system of (11) is the limiting fast system, whose orbits are called fast orbits or singular (boundary and/or internal) layers. By a singular orbit of system (10) or (11), we mean a continuous and piecewise smooth curve in that is a union of finitely many slow and fast orbits. Very often, limiting slow and fast systems provide complementary information on state variables. Correspondingly, the main task of singularly perturbed problems is to patch the limiting information together to form a solution for the entire system.

Let and be the subsets of the phase space defined by

Then the original boundary value problem is equivalent to a connection problem, namely, finding a solution of (10) or (11) from to (see, for example, [48]).

Due to the jumps of the permanent charge function (6) at and we split the interval into three subintervals and where the intervals and represent the reservoirs, and the interval represents the channel. To construct a singular orbit over the whole interval , we first construct a singular orbit on each of the subintervals. To get started, we preassign the values of and at and as follows:

These eight unknown values will be determined along our construction of a singular orbit on the whole interval

The singular orbit on consists of two boundary layers and and one regular layer with being

In particular, given the flux densities and the value are uniquely determined (see Section 2.1.1).

The singular orbit on consists of two boundary layers and and one regular layer with being

In particular, given and , the flux densities and the value and are uniquely determined (see Section 2.1.2).

The singular orbit on consists of two boundary layers and and one regular layer with being

In particular, given , the flux densities and the value are uniquely determined (see Section 2.1.3).

To obtain a singular orbit on the whole interval , one need

This consists of eight conditions. The number of conditions is exactly the same as the number of unknowns in (13).

The singular orbit constructed for problem (10) associated to and consists of nine pieces . Once a singular orbit is constructed, one then can apply the geometric singular perturbation theory, such as Exchange Lemma, to show that, for but small, there is a unique solution that is close to the singular orbit.

In this work, we will examine the competition between cations due to the nonlinear interaction among channel geometry, small permanent charges, diffusion coefficients and boundary conditions, which can be extracted from the matching conditions (14) (for simplicity, we still use for the zeroth order system in ).

We remark that , , , , , , are actually the functions of the unknowns with the parameter . Furthermore, for simplicity, in the following analysis, we will use to denote and , respectively.

Once a solution for is determined, one then can derive the zeroth order (in ) individual fluxes . Through our following discussions, we always assume the so-called electroneutrality boundary concentration conditions

For simplicity, we also introduce

2.1.1. Singular Orbit on with

Following (13), we introduce

We now construct a singular orbit on that connects to which generally consist of two boundary layers and a regular layer (see [7,10,11,43]). Over the subinterval , the permanent charge is zero because we review as one of the reservoirs. However, the nonzero over the subinterval will affect the solution on and (another reservoir with zero permanent charge) through the matching conditions imposed on and to construct the singular orbit over the whole interval

Limiting fast dynamics and boundary layers on Setting in (10), we get the so-called slow manifold,

Setting in (11), we get the limiting fast system,

Note that the slow manifold is the set of equilibria of (17). The following can be established directly [11].

For system (

The slow manifold is precisely the set of equilibria of (17). The linearization of (17) at each point of has seven zero eigenvalues whose generalized eigenspace is the tangent space of the seven-dimensional slow manifold of equilibria, and the other two eigenvalues are , whose eigenvectors are not tangent to (Recall that ’s are concentrations and we are only interested in positive ones). Thus, is normally hyperbolic. □

The theory of normally hyperbolic invariant manifolds [49] states that there exists eight-dimensional stable manifold of that consists of points approaching in forward time; and there exists eight-dimensional unstable manifold of that consists of points approaching in backward time. Let be the collection of orbits from in forward time under the flow of system (17) and be the collection of orbits from in backward time under the flow of system (17). Then, for a singular orbit connecting to , the boundary layer at must lie in and the boundary layer at must lie in . In this subsection, we will determine the boundary layers and , and their landing points and on the slow manifold . The regular layer, determined by the limiting slow system, will lie in and connect the landing points at and at .

First, one has

The following functions are the first integrals of system (

This can be verified directly. □

For the landing points and , following the similar outline as those in [7,11], one has

Assume the condition (

The stable manifold where

The unstable manifold where and

The boundary layer together with

Similarly, the boundary layer together with

Limiting slow dynamics and regular layers on For convenience, we introduce and defined as

Next we construct the regular layer on that connects and . Note that, for , system (10) loses most information. To remedy this degeneracy [7,10,11], we make a rescaling and in system (10). In term of the new variables, system (10) becomes

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

For system (20), the slow manifold is Therefore, the limiting slow system on is given by

For system (21), one has

There is a unique solution where where

Adding the second equation to the third one in (21), one has which gives Substituting (24) into the second equation in (21) to get By the variation of constants formula, we obtain Similarly, can be obtained.

Substituting (24) into the first equation in (21) to get The solution is Evaluating , and at yield the formulas for and . □

The slow orbit given in Lemma 2 connects and . Let (resp., ) be the forward (resp., backward) image of (resp., ) under the slow flow (21). One has the following result (the proof follows exactly the same line as Proposition 3.7 in Section 3.1.2 of [43]).

On the nine-dimensional slow manifold

2.1.2. Singular Orbit on with

Following (13), we let and construct a singular orbit on connecting to .

Limiting fast dynamics and boundary layers on By a similar argument as in Section 2.1.1, one has the slow manifold and the corresponding limiting fast system,

For system (26), similar argument shows that the slow manifold is normally hyperbolic. We denote the stable (resp. unstable) manifold of by (resp. ). Let be the collection of orbits from in forward time under the flow of system (26) and be the collection of orbits from in backward time under the flow of system (26). Then, for a singular orbit connecting to , the boundary layer at must lie in and the boundary layer at must lie in . In this section, we will determine the boundary layers and , and their landing points and on the slow manifold . The regular layer, determined by the limiting slow system, will lie in and connect the landing points at and at .

Similarly, we have the following result.

System (

The stable manifold where where where

The unstable manifold where where where

The boundary layer together with

Similarly, the boundary layer together with

Limiting slow dynamics and regular layers on We now turn to the study of the flow in the vicinity of the slow manifold . Following a similar argument as that in Section 2.1.1, we make a scaling and System (10) becomes It is again a singular perturbation problem and its limiting slow system is For system (31), the slow manifold is Therefore, the limiting slow system on is given by

Note that, on one has It follows that

Since and system (32) has the same phase portrait as that of the following system obtained by multiplying on the right-hand side of system (32):

There is a unique solution where

The proof is similar as that of Lemma 2. As for the system (

It follows that the regular layer solution on is given by (34) with and determined by (35). Together with the boundary layers and described in statement (iv) of Proposition 4, this gives the singular orbit on the interval

2.1.3. Singular Orbit on with

Th construction of singular orbits over from to is virtually the same as the construction of singular orbits on studied in Section 2.1.1. Instead of repeating the same process, we will state only the results for later use.

Limiting fast dynamics and boundary layers on The limiting fast system reads The slow manifold is which consists of the equilibria of system (36) and is normally hyperbolic with a eight-dimensional center-stable manifold and a eight-dimensional center-unstable manifold . For the boundary layers, one has the following proposition.

Under the condition (

System (

The stable manifold where and

The unstable manifold where

The boundary layer together with

Similarly, the boundary layer together with

Limiting slow dynamics and regular layers on We now examine the existence of the regular layer on that connects and . It follows from exactly the same analysis as that in Section 2.1.1, the limiting slow dynamics reads For system (38), one has

There is a unique solution where where

The regular solution that connects and , together with and in statement (iv) of Proposition 5 yields the singular orbit on .

2.1.4. Matching and Singular Orbits on the Whole Interval

A singular orbit over the whole interval will be the union of the singular orbits constructed on each of the subintervals (see Figure 1). The matching conditions are , and and have to be the same on all subintervals; that is, from (18), (23), (27)–(30), (35), (37) and (39), where, for

Figure 1

Schematic picture of a singular orbit projected to the space of with defined in (6).

Recall that and are the unknown values preassigned at and , and are the unknown values for the flux densities of the three ion species. There are also three auxiliary unknowns and in (40). The total number of unknowns in (40) is fourteen, which matches the total number of equations.

A qualitative important question is whether the set of nonlinear Equation (40) (a governing system) has a unique solution. This can be studied through bifurcation analysis and numerical simulations. However, this is beyond the aim of this work.

2.1.5. Existence of Solutions near the Singular Orbit

Note that any solution of the set of algebraic equations determines a singular orbit for the connection problem. Once a singular orbit is constructed, one can apply geometric singular perturbation theory to show that, for small, there is a unique solution that is close to the singular orbit.

For our case, the singular orbit consists of nine pieces: two boundary layers and ; four internal layers and ; and three regular layers and (see Figure 1). More precisely, with ,

The boundary layer connects the point to the point

The regular layer connects the point to the point

The internal layer connects the point to the point

The internal layer connects the point to the point

The regular layer connects the point to the point

The internal layer connects the point to the point

The internal layer connects the point to the point

The regular layer connects the point to the point

The boundary layer connects the point to the point

The following result can be established by the exchange lemma (see, for example, [48,50,51]) of the geometric singular perturbation theory (see also [7,10,11,12]).

Let

Fix small to be determined. Let which is a neighborhood of in .

For , let be the forward trace of under the flow of system (10). We will show that intersects transversally near the point

The evolution of will undergo 3 stages one over each subinterval and . In the first stage over , it will start close to the point , follow the singular layer toward the slow manifold , move along the regular layer , and leave the vicinity of along the singular layer toward the point . It then continues the evolution over each subsequent subintervals in a similar fashion until it reaches the vicinity of the point .

To track this evolution of , we apply the exchange lemma successively (three times) along the stages in order described above. During the first stage, we track the evolution of along the singular orbit . The Exchange Lemma ([48,50,51], etc.) indicates that, at the end of the first stage and near , is -close to for some independent of , provided that the following conditions are satisfied:

is transversal along , which is established in Proposition 2;

the vector field on is not tangent to at , which follows from in (21).

Let then, near , is close to the forward trace of under the flow of system (10) with . We can then apply the Exchange Lemma again to along 0ver . At the end of this stage, one has that is -close to

Let then, near , is close to the forward trace of under the flow of system (10) with . We can then apply the Exchange Lemma again to along 0ver . At the end of this stage, one has that is -close to Since the latter is transversal to near the point , intersects transversally near . Note that and . Therefore, that is, the intersection near is a singleton. This completes the proof. □

2.2. Regular Perturbation Analysis: Expansions along Small

We expand all unknown quantities in the governing system (40) and (41) in under the assumption that is small, for example, for we write

We will determine the coefficients of the zeroth order and first order terms for dominating effects on ionic flows from the permanent charge. We introduce and, corresponding to (16) and (19),

Careful calculations (tedious but straightforward) give

The zeroth order solution in where

The first order solution in where

For convenience, we introduce and , which are defined by

One has

if

if

From Propositions 6 and 7, one has

Assume the electroneutrality boundary conditions ( where

For convenience in our following analysis, we introduce a function for with

For , one can easily established the following properties.

For

Set

We define the critical potentials where

We would like to point out that for the discussion in Section 3.1, the sign of the term is critical. While the sign of A can be determined by as stated in Lemma 7, we now establish the result, which characterizes the sign of .

Assume

if

if

if

The statement that follows from the fact that both A and have the same sign with that of .

Rewrite as where It follows directly that . Careful computations yield where and further

One has for all , that is, is concave upward for . Furthermore, one has and , which indicate that for . Note that for . Then, for , and is strictly increasing on . Note also that . One has, for ,  □

if , then, , which implies ; this completes the proof of statement (i).

For the case with , we first claim that there exists a unique such that for . In fact, based on the facts that and is strictly increasing on , it suffices to show that for , which follows from . For convenience, for , we set

Direct calculations give for all , and hence, is concave downward for . Further, is implied by , since . To prove , set

It is easy to check that and . Further one arrive at the conclusion that based on the facts that and for .

If , then one can easily obtain for all , and more specifically, , which implies , for , , for , , which indicates , for ; this completes the proof of statement (ii).

if , then, the straight line and have the unique intersection , which means that there exists a unique such that , which suggests , for , and further , for , , which yields , for , , for , , which hints , for ; this completes the proof of statement (iii).

3. Results

3.1. Competitions between Cations

We now consider the competition between two positively charged ion species due to the small permanent charges, which further depends on the nonlinear interplays with other system parameters, such as channel geometry , diffusion coefficients and boundary concentrations . The study is closely related to the selectivity phenomena of ion channels.

For convenience, we define as where and are given by where and is defined in (16). In particular, for , one has

is the leading term that contains the effect from small permanent charges, and will be our main focus. The study on could provide important insights and better understanding of the selectivity of cations, and meanwhile demonstrates the critical role played by the small permanent charge in the selectivity of ion channels. The two expressions for will be alternatively chosen for convenience in our following discussion.

In terms of the interplay between the diffusion coefficients and boundary concentrations, there are six cases to discuss, more precisely, one has In this work, we will just focus on case (i), (iii) and (v), since the discussions for case (ii), (iv) and (vi) are very similar. Interested readers can work on them following our arguments.

In our following discussion, we examine from two directions for each case: the sign of and the monotonicity of in the electric potential V.

3.1.1. Case Study with

In this section, we study under the condition which is equivalent to

We first consider a special case.

Suppose

With and , directly from (49), one has From Lemmas 5 and 7, together with and , one has (resp. ) if (resp. ). To see the monotonicity of in the potential V, we calculate where Direct calculation shows that is concave downward for all V and attains its global maximum 0 at . Note that One then has (resp. ) if , that is, (resp. , that is, ). This completes our proof. □

We next consider the more general case with .

Suppose that

if

if

The statements follow directly from (49) and Lemma 5. □

Suppose that

if

if

We will just provide a detailed proof for the first statement, and the second one can be discussed similarly. With , one has It follows that where and further

It is easy to see that has two zeros given by and . Furthermore, one can easily verify that if ; if ; and if .

If , then, decreases on , increases on , and decreases on . Note that , which is a local minimum of , and . The function has a second zero, denoted by , with . Furthermore, for , we have since , which can be obtained from the fact that decreases on and . Therefore, if , and if . This holds for the case with by a similar argument. Note that and . Statement (i) then follows. □

3.1.2. Case Study with

We study the term , which provides information of the preference of the ion channel over different cation under the condition .

We first characterize the sign of under the further restriction , and the order of the critical potentials and .

Let

Let

Let

We comment that the proofs of Corollaries 2–4 can be directly discussed based on Lemmas 7 and 8. We skip the details here.

Suppose that

Suppose that

For

For

For

Suppose that

If

If

3.1.3. Case Study with

In this section, we study the term for the case with , which is equivalent to and . By similar arguments as those in Section 3.1.1 and Section 3.1.2, we obtain the following results.

Suppose that

Suppose that

For

For

3.1.4. Studies on the Magnitude of

For convenience in the argument, we let represent the first cation corresponding to the flux and be the cation corresponding to the flux .

Recall that, with small, one has . Further depending on the interaction among and , in Section 3.1.1, Section 3.1.2 and Section 3.1.3, we analyze the leading term that contains the effects from small permanent charges, in particular,

the sign of , which characterizes the small positive permanent charge effects on the competition between two cations. To be specific, if (resp. ), then, the small positive permanent charge enhances (resp. reduces) , and in either way, it affects the preference of the ion channel over different cation, which is closely related to the selectivity phenomena of the ion channel.

the monotonicity of in the electric potential V, which further reduces or strengthens the preference by adjusting/controlling the boundary membrane potential. Taking the case for example, if , then, one can increase the boundary electric potential to further strengthen the individual flux , which indicates that more cation will go through the ion channel; while if , one then need to decrease the boundary electric potential for more cation to go through the ion channels.

On the other hand, (resp. ), indicates (resp. ), but it does not provide any information on the relation of and , which contains complementary information for the competition, and further depends on the sign of .

Recall form Lemma 7 that if , together with (48) and (49), one has

Suppose that

Suppose that

if either

for

for

if

for

for

From (48), one note that the sign of is determined by the sign of . Together with Lemmas 7 and 8, one can easily established the result. □

We point out that the studies provides further information of the preference of the ion channel over distinct cations. To be specific, we take for example.

if and , then, the ion channel prefers the cation over the cation , and the small positive permanent charge further enhances this preference;

if and , then, the ion channel prefers the cation over the cation , and the small positive permanent charge further enhances this preference.

Similar argument can be applied to the case .

3.2. Numerical Simulations

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

In our simulations to system (50) and (51), we take

The choice of

Our numerical simulations consist of two parts focusing on some cases discussed in Section 3.1.1, Section 3.1.2 and Section 3.1.3, respectively, and further verify some analytical results stated in Theorems 3, 4, 6, 7 and 9 for some carefully selected system parameters. Other related results can also be numerically illustrated by choosing different parameter values, and we leave that to interested readers.

Recall that It turns out that our numerical simulations with nonzero but small are consistent with our analytical results. To be specific,

By choosing and , one has

It follows that which is consistent with the conditions stated in (ii) of Lemma 8, and indicates that . Together with Lemma 7, we have . Therefore, , which is consistent with the condition (i) stated in Theorems 3 and 4, respectively.

Our numerical results show that (see Figure 2) , approximation of defined in (48), which is given by has a unique zero , and a critical point such that if and if ; furthermore, increases in the potential V if , and decreases in the potential V if .

Figure 2

Identification of critical potentials and and approximation of with for .

We remark that the numerical result is consistent with our analytical result, more precisely, statement (i) of both Theorems 3 and 4.

By choosing and , one has

It follows that which is consistent with the conditions stated in (ii) of Corollary 4.

Our numerical results show that (see Figure 3) , approximation of defined in (48), which is given by has two zeros and with , and a critical point between and such that if or , and if ; furthermore, decreases in the potential V if , and increases in the potential V if .

Figure 3

Identification of critical potentials and and approximation of with for .

We remark that the numerical result is consistent with our analytical result, more precisely, statement (iii) of Theorem 6 and statement (i) of Theorem 7.

By choosing and , one has

It follows that

From Lemma 7, one has and from Lemma 8 (statement (ii1)), one has , and hence, , which satisfies the condition required in the statement (i) of Theorem 9.

Our numerical results show that (see Figure 4) has one zero and two critical points and with such that if and if ; furthermore, increases in the potential V if , decreases in the potential V if , and increases in the potential V if .

Figure 4

Identification of critical potentials and and approximation of with for .

We remark that the numerical result is consistent with our analytical result stated in (i) of Theorem 9.

4. Concluding Remarks

We study a one-dimensional steady-state Poisson–Nernst–Planck system with three ion species, two cations having the same valence and one anion. Nonzero but small permanent charge, the major structural quantity of an ion channel, is included in the model. The cross-section area of the channel is included in the system, and it provides certain information of the geometry of the three-dimensional channel, which plays crucial roles in our analysis. Two specific structures of the PNP model allows one to reduce the singularly perturbed boundary value problem to an algebraic system-the governing system (40). The significance of the governing system is:

the existence of a complete set of first integrals (Proposition 1, first statement (i) in Proposition 4 and Proposition 5, respectively);

the observation made in Section 2.1.2 allows one to transform the limit slow system (32) to a linear system (33) with constant coefficients;

it includes almost all relevant physical parameters;

once a solution of the governing system is obtained, the singular orbit (the zeroth order approximation (in ) solution of the boundary value problem) can be readily determined.

Based on these specific structures of this concrete model, under the framework of the geometric singular perturbation theory, a singular orbit is obtained, from which explicit expressions of and are extracted. This makes it possible for one to further analyze the competition between cations, which further depend on the complicated nonlinear interplays among system parameters, such as the diffusion coefficients , the channel geometry in terms of , the boundary conditions , particularly, the interaction among and plays a critical role in characterizing the competition between cations (Section 3.1.1, Section 3.1.2 and Section 3.1.3). Among others, we find

As functions of the membrane potential V (fixing other system parameters),

and can be positive (resp. negative), which further depends on the nonlinear interaction among boundary concentrations and diffusion coefficients. The sign of provides critical information related to the selectivity phenomena of ion channels, while the sign of provides efficient ways to further adjust/control the preference of the ion channel over distinct cation (Characterized in Theorems 2–9);

, the magnitude of , which is equivalent to for small positive , is analyzed, which provides further information of the ion channel over distinct cations (Theorems 10 and 11).

Critical potentials that balance the small permanent charge effects on (such as and ) are identified (Definition 1). Those critical potentials can be experimentally estimated. Taking for example, one can take an experimental curve as and numerically (or analytically) compute for ideal case that allows one to get an estimate of by considering the zeros of . The critical potentials play critical roles in characterizing permanent charge effects on ionic flows through membrane channels.

Finally, we comment that the setup in this work is relatively simple, and may raise the concern about the feasibility. Indeed, cPNP is known to be reliable when the ionic mixture is dilute, but with more ion species and nonzero permanent charges included, the ionic mixture would be crowded. On the other hand, the setup is reasonable for semi-conductor problems and for synthetic channels. The detailed analysis in this work is not only important in understanding the behavior of such fluxes across channels but can also better help in understanding the fundamental mechanistic principles of metabolite fluxes in engineered synthetic biological channels [52]. Furthermore, the study in this work is the first step for analysis of more realistic models. The simple model considered in this work allows us to obtain more explicit expressions of the ionic fluxes in terms of physical parameters of the problem so that we are able to extract concrete information on small permanent charge effects. Moreover, the analysis in this simpler setting provides important insights for the analysis and numerical studies of more realistic models.