Geometric Mean of Concentrations and Reversal Permanent Charge in Zero-Current Ionic Flows via Poisson-Nernst-Planck Models
GGeometric Mean of Concentrations and Reversal PermanentCharge in Zero-Current Ionic Flows via Poisson-Nernst-PlanckModels
Hamid Mofidi ∗ September 22, 2020
Abstract
This work examines the geometric mean of concentrations and its behavior in varioussituations, as well as the reversal permanent charge problem, the charge sharing seen in x-ray diffraction. Observations are obtained from analytical results established using geometricsingular perturbation analysis of classical Poisson-Nernst-Planck models. For ionic mixtures ofmultiple ion species Mofidi and Liu [
SIAM J. Appl. Math. (2020), 1908-1935 ] centered twoion species with unequal diffusion constants to acquire a system for determining the reversalpotential and reversal permanent charge. They studied the reversal potential problem andits dependence on diffusion coefficients, membrane potential, membrane concentrations, etc.Here we use the same approach to study the dual problem of reversal permanent charges andits dependence on other conditions. We consider two ion species with positive and negativecharges, say Ca + and Cl − , to determine the specific conditions under which the permanentcharge is unique. Furthermore, we investigate the behavior of geometric mean of concentrationsfor various values of transmembrane potential and permanent charge. Key words.
Ionic flows, PNP, fast-slow systems, concentrations, permanent charge
The nervous system is a too complicated and intricate part of an animal that is especially criticalfor transmitting signals between different body parts. It recognizes environmental changes thataffect the body, then works together with other body systems to respond to detected changes.It contains a large number of neurons that are electrically excitable cells. Neurons send electricsignals to each other through thin fibers called axons, which generates chemicals known as neu-rotransmitters to be delivered at synapses. This electric signal, propagating along the axon, isa variation of polarization of transmembrane electrostatic potential called an action potential.An action potential is induced by the opening and closing many ion channels distributed on theaxon membrane. Ion channels, proteins embedded in membranes, provide a major pathway forcells to communicate with each other and with the outside to transform signals and to conductgroup tasks ([6, 14, 25, 26]). The essential structure of an ion channel is its shape and its perma-nent charge. The shape of a typical channel could be approximated as a cylindrical-like domain.Within an ion channel, amino acid side chains are distributed mainly over a “short” and “narrow”portion of the channel, with acidic side chains contributing negative charges and basic side chains ∗ Department of Mathematics, University of Iowa, Iowa City, IA 52242, ( [email protected] ). a r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p roviding positive charges. It is specific of side-chain distributions, which is referred to as thepermanent charge of the ion channel. The function of channel structures is to select the types ofions and to facilitate the diffusion of ions across cell membranes.The multi-scale feature of the problem with multiple physical parameters enables the system tohave high flexibility and to show rich phenomena/behaviors ([2, 13]). On the other hand, the samemulti-scale characteristic with multiple physical parameters presents a remarkably demanding taskfor anyone to derive meaningful information from experimental data, also given the fact that theinternal dynamics cannot be discerned with the present technique.To reveal mechanisms of perceived biological phenomena and explore new aspects, the role ofmathematical analysis is inevitable. There have been some successes recently in probing Poisson-Nernst-Planck (PNP) models for ionic flows through ion channels [16, 18, 32, 34, 39, 40, 42, 48].Centering specific critical characteristics of the biological systems, PNP models serve suitably foranalysis and numerical simulations of ionic flows. One can acquire PNP systems as diminishedmodels from molecular dynamic models, Boltzmann equations, and variational principles [3, 27,28, 51]. There are various types of PNP models:(i) The classical PNP treats dilute ionic mixtures, where no ion-to-ion interactions are involved. (ii)
The hard-sphere PNP reflects volume exclusive by employing ions as hard-spheres.
More sophisticated models have also been studied in [2, 10, 15, 50], etc. It is challenging, though,to achieve analytical or computational results from complicated models.In this work, we are engaged in geometric mean of concentrations and reversal permanentcharges that are determined by zero total currents. We study the connection of these quantitieswith the membrane potentials or diffusion constants. The total current I = I ( V, Q ) depends on thetransmembrane potential V and the permanent charge Q . For fixed transmembrane potential V ,a reversal permanent charge Q = Q rev ( V ) is a charge that generates zero current I ( V, Q rev ( V )) =0. We employ the classical PNP model and consider a cylinder-like channel to fulfill the basicunderstanding of plausible effects of general diffusion coefficients in ionic channels.To appreciate the significance of permanent charges in ionic channels, we emphasize that per-manent charges in ionic channels perform the role of doping profiles in semiconductor devices.Doping gives the charges what acid and base side chains provide in ionic channels. Both ionicchannels and semiconductor devices employ atomic-scale constructions to regulate macroscopicflows from one pool to another. Holes and electrons are the cations and anions of semiconduc-tors, respectively. Ions usually flow as quasi-particles flow in semiconductors that depends oncontrolling movement and diffusion of quasi-particles of charge in transistors and integrated cir-cuits. Doping is the process of adding impurities into primary semiconductors to strengthen itselectrical, optical, and structural features [8, 49, 54].The role of diffusion constants, D j ’s is also essential. The authors of [18] explored the problemof determining reversal permanent charges for the case when all diffusion constants are equal.However, the identical diffusion constants case is degenerate, known from the biological perspec-tive. The problem with unequal diffusion coefficients has been considered in some works. In [2],the authors discussed how mobilities and their spatial inhomogeneities are affected by other pa-rameters. In [1], the authors carried a perturbation inquiry from a time-independent and spatiallyhomogeneous equilibrium solution. Two time scales of the dynamics are classified from the O ( (cid:15) )terms. In particular, for the equal diffusion coefficient, the authors show that the diffusion processfor O ( ε ) terms does not occur – a vital effect of unequal ionic mobilities. In [23], the authorsreviewed the cases with unequal mobilities by calculations of a wholly nonlinear electrokineticmodel. They recognized the appearance of a steady long-range field due to unequal mobilities.In this work, we are mainly inspired by the effect of unequal diffusion coefficients and bound-2ry concentrations on the geometric mean of concentrations and reversal permanent charges forthe zero-current problem. In [44], the authors used the geometric singular perturbation (GSP)framework developed in [16, 39, 40] for analyzing PNP models for ionic flow to arrange a systemof algebraic equations for the problem. The difference between D and D makes the systembecomes a complex nonlinear algebraic system that is further reduced to two nonlinear equationsthat appeared to work satisfactorily and allow one to explore the zero-current problem. We utilizethe same structure in this paper.To underline the leading portions of this document (as well as some in [34, 44, 45]), wefirst desire to stress that employing the geometric analysis enables one to express and classifyquantities and properties that are crucial to biology, and additionally to present the quantitativeand qualitative perception and predictions. In this work, we show a derivation of a mathematicalsystem for the zero-current problem that we use to learn reversal permanent charge and thegeometric mean of concentrations and their dependence on the other parameters like membranepotential, boundary concentrations, and diffusion constants. Some numerical simulations havebeen provided throughout the text to support the theoretical conclusions and give the interestedreader a sharp comprehension of the claims.Throughout the paper, we obtain the numerical results from the algebraic systems (3.2) and(3.3) that have been obtained from the governing system (2.11). The nonlinear algebraic systemsare solved by Matlab R (cid:13) (Version 9.5) function fsolve that uses the trust-region dogleg algorithm,that is based on the interior-reflective Newton method defined in [11].This paper is constructed as follows. The classical PNP model for ionic flows is recalled inSection 1.1 to set the stage for analyses in the next sections. In Section 2, we apply the GSPtheory on the PNP system, with zero-current, to convert the BVP to a connecting system to endup with a nonlinear algebraic system of equations, called the matching system. In Section 3, wereduce the matching system to two nonlinear equations. In particular, we analyze the geometricmean of concentrations and reversal permanent charge. We review some concluding remarks inSection 4. The PNP equations have been simulated and computed to a considerable extent [9, 12, 4, 30, 31].From those simulations, one can see that mathematical boundary conditions, i.e., macroscopicreservoirs, need to be incorporated in the mathematical formulation to describe the actual behaviorof channels [22, 47]. For an ionic mixture of n ion species, the PNP model is, for k = 1 , , ..., n ,Poisson: ∇ · (cid:16) ε r ( −→ X ) ε ∇ Φ (cid:17) = − e (cid:16) n (cid:88) s =1 z s C s + Q ( −→ X ) (cid:17) , Nernst-Planck: ∂ t C k + ∇ · −→J k = 0 , −−→J k = 1 k B T D k ( −→ X ) C k ∇ µ k , (1.1)where −→ X ∈ Ω (which is a three-dimensional cylindrical-like domain) representing the channelof length ˆ L nm(= ˆ L × − m), ε r ( −→ X ) is the relative dielectric coefficient (with unit 1), ε ≈ . × − Fm − is the vacuum permittivity, e ≈ . × − C (coulomb) is the elementarycharge, Q ( −→ X ) represent the permanent charge density of the channel (in M = Molar = mol/L), k B ≈ . × − JK − is the Boltzmann constant, T is the absolute temperature ( T ≈ .
16 K =kelvin, for water), and Φ is the electric potential (with the unit V = Volt = JC − ).For the k -th ion species, C k is the concentration (with unit M), z k is the valence (the number of3harges per particle with unit 1), and µ k is the electrochemical potential (with unit J = CV). Theflux density −→J k ( −→ X ) (with unit mol m − s − ) is the number of particles across each cross-section inper unit time, D k ( −→ X ) is the diffusion coefficient (with unit m / s), and n is the number of distincttypes of ion species (with unit 1).Since ion channels have thin cross-sections comparative to their lengths, three-dimensionalPNP systems can be reduced to quasi-one-dimensional models ([41]). The quasi-one-dimensionalsteady-state PNP model is, for k = 1 , , ..., n, A ( X ) ddX (cid:18) ε r ( X ) ε A ( X ) d Φ dX (cid:19) = − e (cid:32) n (cid:88) s =1 z s C s + Q ( X ) (cid:33) ,d J k dX = 0 , −J k = 1 k B T D k ( X ) A ( X ) C k dµ k dX , (1.2)where X is the coordinate along the channel, A ( X ) is the area of cross-section of the channel overlocation X , and J k (with unit mol/s) is the total flux through the cross-section. We apply thefollowing boundary conditions to the system (1.2), for k = 1 , , · · · , n ,Φ(0) = V , C k (0) = L k >
0; Φ( ˆ L ) = 0 , C k ( ˆ L ) = R k > . (1.3)One often uses the electroneutrality conditions on the boundary concentrations because the solu-tions are made from electroneutral solid salts, n (cid:88) s =1 z s L s = n (cid:88) s =1 z s R s = 0 . (1.4)There is a sharp layer for electrical potential and ion concentrations near the encounter betweenthe zero and non-zero permanent charges. This occurs near x = a and x = b , in our calculations,where the permanent charge jumps (see how Q ( x ) is defined in (1.10)). Since our boundaryconditions are forced at x = 0 < a and at x = 1 > b , one can disregard the changes in theboundary concentrations. Nevertheless, our approach can regulate the boundary layers evenwhen the boundary conditions do not fit the electroneutrality condition [16].The electrochemical potential µ k ( X ) for the k -th ion species consists of the ideal component µ idk ( X ) and the excess component µ exk ( X ), i.e., µ k ( X ) = µ idk ( X ) + µ exk ( X ). The excess elec-trochemical potential µ exk ( X ) accounts for the finite size effect of ions. It is needed wheneverconcentrations exceed, say 50 mM, as they almost always do in technological and biological situ-ations and often reach concentrations 1M or more. The classical PNP model only deals with theideal component µ idk ( X ), which disregards ions-size and displays the dilute ions-entropy in water.Dilute solutions tend to approach ideality as they proceed toward infinite dilution; that is, µ k ( X ) = µ idk ( X ) = z k e Φ( X ) + k B T ln C k ( X ) C , (1.5)where one may take C = max ≤ k ≤ n (cid:8) L k , R k , sup X ∈ [0 , ˆ L ] |Q ( X ) | (cid:9) as the characteristic concentra-tion of the problems.For given V , Q ( X ), L k ’s and R k ’s, if (Φ( X ) , C k ( X ) , J k ) is a solution of the boundary valueproblem (BVP) of (1.2) and (1.3), then the electric current I is I = e (cid:80) ns =1 z s J s . For an analysisof the boundary value problem (BVP) (1.2) and (1.3), we work on a dimensionless form. Set D = max ≤ k ≤ n { sup X ∈ [0 , ˆ L ] D k ( X ) } and ¯ ε r = sup X ∈ [0 , ˆ L ] ε r ( X ). Then let4 = ¯ ε r ε k B Te ˆ L C , ˆ ε r ( x ) = ε r ( X )¯ ε r , x = X ˆ L , h ( x ) = A ( X )ˆ L , D k ( x ) = D k ( X ) D ,Q ( x ) = Q ( X ) C , φ ( x ) = e k B T Φ( X ) , c k ( x ) = C k ( X ) C , ˆ µ k = 1 k B T µ k , J k = J k ˆ LC D . (1.6)In terms of the new variables, the BVP (1.2) and (1.3) become, for k = 1 , , · · · , n , ε h ( x ) ddx (cid:18) ˆ ε r ( x ) h ( x ) ddx φ (cid:19) = − n (cid:88) s =1 z s c s − Q ( x ) ,dJ k dx = 0 , − J k = h ( x ) D k ( x ) c k ddx ˆ µ k , (1.7)with the boundary conditions φ (0) = V = e k B T V , c k (0) = l k = L k C ; φ (1) = 0 , c k (1) = r k = R k C . (1.8) Remark . It is reasonable to assume that ε > L =2 . C = 10 M, then we obtain ε ≈ − [17]. The smallness of ε will later let ustreat the system (2.7) of the dimensionless problem as a singularly perturbed problem that canbe analyzed by the GSP theory. As we will discuss in more detail later in Section 2, the GSPTheory employs the modern invariant manifold theory from nonlinear dynamical system theoryto examine the entire structure, i.e., the phase space portrait of the dynamical system, and shouldnot be confused with the classical singular perturbation theory that uses, for example, matchedasymptotic expansions.As seen in Fig 1, for any point x ∈ [0 , h ( x ) (the right panel) is the cross-section area of thechannel in a dimensionless form corresponding to X ∈ [0 , ˆ L ] in Ω, (the left panel) which is thethree-dimensional form of the channel.Figure 1: The left panel shows Ω( x ) , which is a three-dimensional cylindrical like domain repre-senting a channel of length ˆ L ; the right panel shows h ( x ) which is the dimensionless form of thearea of cross-section of the channel.Remark . At this step, we would like to emphasize some tips:(i) We first recall that dimensionless variables are suitable for explaining mathematical andgeneral physical relations, and dimensional quantities are required to reveal how evolution has5tilized those relations. We frequently switch from dimensional form to the dimensionless one andconversely throughout the script. The quantities
V, l, r, c k , Q, D k , and J k are dimensionless valuescorresponding to the dimensional quantities V , L, R, C k , Q , D k , and J k , respectively, acquired from(1.6).(ii) The genuine dimensional forms of quantities have been applied for all figures throughoutthe paper, that is, C k = C c k (M) , Q = C Q ( M ) , J k = ˆ LC D J k (mol / s) , and we take C = 10M, ˆ L = 2 . D = 2 . × − m / s . Also, for diffusion constants D k we consider combinations of the following [38, 53],1 . × − cm / s for Na + , . × − cm /s for Cl − , . × − cm /s for Ca , . × − cm / s for CO − , (1.9)where N a + , Cl − , Ca , and CO − respectively stand for Sodium, Chloride, Calcium and Car-bonate.We now designate the case we will study in this paper. We will investigate a simple setup, theclassical PNP model (1.7) with the boundary conditions (1.8), and ideal electrochemical potential(1.5). More clearly, we assume (A0) The ionic mixture consists of two ion species with valences z = − z = 1 ;(A1) D k ( x ) = D k for k = 1 , is a constant and ˆ ε ( x ) = 1 ;(A2) Electroneutrality boundary conditions (1.4) hold;(A3) The permanent charge Q is piecewise constant defined as, Q ( x ) = (cid:26) Q = Q = 0 , x ∈ (0 , a ) ∪ ( b, ,Q , x ∈ ( a, b ) , (1.10) where Q is a constant. A common struggle for nonlinear differential equations is that it is usually impractical to developspecific solutions since they are very diverse, and examining solutions depends on the problems.Time scale separation is a phenomenon that appears in many physical systems like chemicalreactions, particles in fluids, etc. In many applications, systems are naturally formulated infast-slow systems, that are typically high-dimensional systems of nonlinear ordinary differentialequations. The GSP Theory is a robust device to analyze multiple scale problems. It relieson advances of invariant manifold theory from nonlinear dynamical system theory and is analternative and complementary theory to the classical matched asymptotic theory for singularperturbations. 6 .1 Singularly Perturbed Problems: Fast-Slow Systems.
One may presume that singularly perturbed problems are more complex than regularly perturbedproblems. However, there are significant benefits of singularly perturbed problems over the reg-ularly perturbed problems. One is the approximate decoupling of the full singularly perturbedproblems. In general, IVP of singularly perturbed systems is exponentially (in ε ) sensitive in ini-tial conditions, but BVP is better behaved. There is no globally admitted classification to groupall perturbation problems into categories of regular and singular. Nonetheless, some formal andinformal tactics are listed below ([37]).A differential equation problem involving a small parameter 0 < ε (cid:28) ε , or if it is, the power series has a vanishingradius of convergence.2. The solution does not converge uniformly as ε → ε = 0.3. Substituting a power series expansion in ε yields problems of a “different type” from thatof the original differential equation.4. Substituting a (regular) power series expansion in ε fails, i.e., a problem is singular if it isnot regular.A fast-slow vector field (or ( m, n ) − fast-slow system) is a system of ordinary differential equationstaking the form dxdτ = . x = f ( x, y, ε ) ,ε dydτ = ε . y = g ( x, y, ε ) , (2.1)where f : R m × R n × R → R m , g : R m × R n × R → R n , and 0 < ε (cid:28)
1. Furthermore, the x variables are called slow variables , and the y variables are called fast variables . Setting t = τε gives the equivalent form dxdt = x (cid:48) = εf ( x, y, ε ) ,dydt = y (cid:48) = g ( x, y, ε ) . (2.2)We refer to t as the fast time scale and to τ as the slow time scale.The differential-algebraic equation obtained by setting ε = 0 in the formulation of the slowtime scale (2.1) is called the slow subsystem or slow vector field or limiting slow system: . x = f ( x, y, , g ( x, y, . (2.3)The flow generated by (2.3) is called the slow flow. Furthermore, the parameterized system ofODEs obtained by setting ε = 0 in the formulation of the fast time scale (2.2) is called the fastsubsystem or fast vector field or limiting fast system: x (cid:48) = 0 , y (cid:48) = g ( x, y, . (2.4)The flow generated by (2.4) is called the fast flow. The limiting slow system is also referred to asthe reduced problem and its flow as the reduced flow. The limiting fast system is also referred toas the layer equations or the layer problem. We call the set Z = (cid:8) ( x, y ) ∈ R m × R n : g ( x, y,
0) = 0 (cid:9) , (2.5)7he critical set, or slow manifold or critical manifold if Z is a submanifold of R m × R n .The relation between equilibrium points of the fast flow and the critical manifold Z is partic-ularly simple. Any orbit (or a portion of it) of the limiting slow system is called a singular sloworbit. Similarly, any orbit (or a portion of it) of the limiting fast system is called a singular fastorbit. By a singular orbit of ε > ε >
0, the slow system (2.1) and fast system (2.2) are equivalent; theirlimiting versions (2.3) and (2.4) are totally different but are complement to each other with lowerdimensions. In addition to being lower dimensions, as mentioned above, the limiting subsystemsoften correspond to ideal physical conditions that are easier to analyze in principle. Therefore,the main objective of the GSP theory is what information from the limiting systems can belifted to the full system, particularly how limiting fast dynamics and limiting slow dynamics areinterplayed.
In applying the nonlinear dynamical system theory, particularly the invariant manifold theory,to the study of singularly perturbed problems, one tries to understand essential structures of thephase portrait. A general systematic approach goes as follows. We recall that the slow manifold Z is the set of equilibria of the limiting fast system (2.4). The linearization at each point p onthe slow manifold Z is (cid:32) D x g ( p, D y g ( p, (cid:33) . (2.6)Therefore, λ = 0 is always an eigenvalue of multiplicity n . We call λ the trivial eigenvalue andall other m eigenvalues of (2.6) the nontrivial eigenvalues. Assume that the number of nontrivialeigenvalues in the right complex half-plane, on the imaginary axis, and in the left complex half-plane are m u , m c , and m s respectively. We denote the corresponding unstable, center, and stableeigenspaces by E u , E c , and E s withdim E u = m u , dim E c = m c + n, dim E s = m s . The set Z is called normally hyperbolic if the n × n matrix ( D x g )( p,
0) of first partial derivativeswith respect to the fast variables has no eigenvalues with zero real part for all p ∈ Z . Theorem 2.1. (Fenichel’s First and Second Theorems) ([24, 35])Suppose S is a compact normally hyperbolic submanifold (possibly with boundary) of the criticalmanifold Z of (2.1) and that f, g ∈ C r ( r < ∞ ) , that is they are smooth. Then for ε > sufficiently small:(1) Fenichel’s first theorem: There exists a manifold S ε , O ( ε ) close and diffeomorphic to S , thatis locally invariant under the flow of the full problem (2.1).(2) Fenichel’s second theorem: There exist manifolds W s ( S ε ) and W u ( S ε ) , that are O ( ε ) closeand diffeomorphic to W s ( S ) and W u ( S ) , respectively, and that are locally invariant under theflow of the full problem (2.1). Consider an ( m, n ) − fast-slow system (2.1) and suppose the critical manifold Z in (2.5) isnormally hyperbolic. Since Z is normally hyperbolic, it follows that D y g | Z is invertible. Hence,the implicit function theorem locally admits that Z in (2.5) becomes Z = (cid:8) ( x, h ( x )) ∈ R m + n (cid:9) ,8hich is a graph, where h : R n → R m is a map so that g ( x, h ( x ) ,
0) = 0. To state the ExchangeLemma, we need the following terminology. Two smooth (sub)manifolds M and M of R k intersect transversally at a point q if T q M and T q M together spans T q R k = T q M + T q M . Inthis case, the intersection N = M ∩ M is a smooth submanifold, and dim N = dim M +dim M − dim R k . Furthermore, if M is a submanifold of U , then U and M intersect transversally too.Also, C perturbations of M and M still intersect transversally.Jones and Kopell ([35, 36]), with deep insight, extracted an extremely useful consequence,collectedly called Exchange Lemma. They viewed the normally hyperbolic slow manifold togetherwith its invariantly foliated stable and unstable manifolds as a device and tested what the devicedoes to an invariant manifold as the invariant manifold enters a neighborhood of the slow manifold.Recall the setup of the singularly perturbed problem in the standard form, i.e. equations 2.1—2.4.Moreover, assume that the slow manifold Z is normally hyperbolic, and assume that k eigenvaluesof the n × n matrix g y ( x, h ( x ) ,
0) have positive real part and l eigenvalues of g y ( x, h ( x ) ,
0) havenegative real part with k + l = n . Let U be a neighborhood of Z within which the invariantfoliation structures holds. Let M ε be an invariant manifold of (2.1) and (2.2) for ε > M ε = k + σ with 1 ≤ σ ≤ m . The limiting manifold M is invariant under the limiting fastsystem (2.4). Furthermore, we assume(i). M intersects W s ( Z ) ∪ U transversally. Let N = M ∩ W s ( Z ), thendim N = dim M + dim W s ( Z ) − dim( R m + n ) = ( k + σ ) + ( m + l ) − ( m + n ) = σ. (ii). The ω − limit set ω ( N ) of the limiting fast system is a ( σ − − dimensional submanifold of Z . One should note that ω ( N ) is nonempty because σ − ≥ ω ( N ). We remark that ω ( N ) is not open in Z because σ − < m .The latter statement indicates that under the limiting slow flow denoted by dot, ω ( N ) · (0 , τ ) forany τ > Z of dimension σ = dim ω ( N ) + 1. Theorem 2.2. (Exchange Lemma for normally hyperbolic slow manifold)([35])Assume (A1), (A2) and (A3) hold. Let τ > τ > be given. Then, for ε > small, a portion of M ε is C O ( ε ) − close to W u (cid:0) ω ( N ) · ( τ , τ ) (cid:1) ∩ U . Note that, dim (cid:0) ω ( N ) · ( τ , τ ) (cid:1) = dim ω ( N ) + 1 = σ, and hence, W u (cid:0) ω ( N ) · ( τ , τ ) (cid:1) = dim ω ( N ) + k = σ + k = dim M ε . We remark that all conditions are required in limiting slow and fast systems, but the conclusionis for ε > small. The conclusion on C closeness stresses that not only the manifolds are O ( ε ) − close in the C sense, but also their tangent spaces are O ( ε ) − close to each other.Standard utilization of the Exchange Lemma is connecting orbits (BVP), heteroclinic andhomoclinic orbits. In connecting orbits, M ε could be the flow of the set B L defining the boundarycondition on the left. One then requires to discern whether M ε intersect the other set B R definingthe boundary condition on the right. That is what occurs to the PNP system, as we will addressin the next section. 9 .3 Application of GSP Theory to the BVP (1.7) and (1.8). We rewrite the classical PNP system (1.7) into a standard form of singularly perturbed systemsand turn the boundary value problem to a connecting problem. For more details one can read thepaper [16]. Denote the derivative with respect to x by overdot and introduce u = ε ˙ φ . System (1.7)becomes, for k = 1 , , ε ˙ φ = u, ε ˙ u = − (cid:88) s =1 z s c s − Q ( x ) − ε h x ( x ) h ( x ) u,ε ˙ c k = − z k c k u − ε J k D k h ( x ) , ˙ J k = 0 . (2.7)System (2.7) will be treated as a dynamical system with the phase space R and the indepen-dent variable x is viewed as time for the dynamical system.A GSP structure to examine the BVP of the classical PNP systems was formed first in [16, 39]for ionic compounds with two ion species. The model of ion channel properties involves couplednonlinear differential equations. The GSP theory enables one to make a conclusion about theBVP for ε > ε = 0 limit systems. Another unique structure is that astate-dependent scaling of the independent variable transforms the nonlinear limit slow systemto a linear system with constant coefficients. The coefficients depend on unknown fluxes to bedetermined as part of the whole problem, which is mathematical evidence for its strong dynamics.Consequently, the existence, multiplicity, and spatial profiles of the singular orbits– zeroth orderin ε approximations of the BVP– are reduced to a nonlinear system algebraic equations thatinvolve all relevant quantities together. This system of nonlinear algebraic equations accuratelydraws the physical framework of the problem. Furthermore, it confirms that all quantities interactwith each other, and we will show in this paper quantitatively how some of those transactionsoccur.With its extensions to include some of the effects of ion size, this geometric framework hasproduced several outcomes that are fundamental to ion channel properties [33, 34, 43, 44, 52, 55].The interested readers are referred to the papers mentioned above for more details on the GSPframework for PNP and concrete applications to ion channel problems.For simplicity, we use the letters l , r and Q where l = l = l , r = r = r , Q = 2 Q .Following the framework in [40], because of the jumps of the permanent charge Q ( x ) in (1.10) at x = a and x = b , we divide the formation of a singular orbit on the interval [0 ,
1] into that onthree subintervals [0 , a ] , [ a, b ] and [ b,
1] to convert the boundary value problem to a connectingproblem. We denote C = ( c , c ) T and J = ( J , J ) T , and preassign values of φ and C at x a = a and x b = b : φ ( x j ) = φ j and C ( x j ) = C j for j ∈ { a, b } . Now for j ∈ { l, a, b, r } , let B j be the subsets of the phase space R defined by B j = (cid:110) ( φ, u, C, J, w ) : φ = φ j , C = C j , w = x j (cid:111) . Note that the sets B l and B r are associated to the boundary condition in (1.8) at x = 0and x = 1 respectively. Thus, the BVP (1.7) and (1.8) is equivalent to the following connectingorbit problem: finding an orbit of (2.7) from B l to B r (See Figure 2). The construction wouldbe accomplished by finding first a singular connecting orbit – a union of limiting slow orbits andlimiting fast orbits, and then applying the exchange lemma to show the existence of a connecting10rbit for ε > B l to B a , one from B a to B b , and one from B b to B r with amatching of ( J , J ) and u at x = a and x = b (see [40] for details). Figure 2:
An illustration of a singular connecting orbit projected to the space of ( u ; z c + z c ; x ) .The solid line from the left boundary B l to the right boundary B r is the O ( ε ) estimate of theconnected problem obtained by Exchange Lemma (see [16]). By setting ε = 0 in system (2.7), we get the slow manifold, for k = 1 , , Z k = (cid:110) u = 0 , z c + z c + Q k = 0 (cid:111) . In terms of the independent variable ξ = x/ε , we obtain the fast system of (2.7), for k = 1 , φ (cid:48) = u, u (cid:48) = − z c − z c − Q j − ε h w ( w ) h ( w ) u,c (cid:48) k = − z k c k u − ε J k D k h ( w ) , J (cid:48) = 0 , w (cid:48) = ε, (2.8)where prime denotes the derivative with respect to ξ . The limiting fast system is, for k = 1 , φ (cid:48) = u, u (cid:48) = − z c − z c − Q j , c (cid:48) k = − z k c k u, J (cid:48) = 0 , w (cid:48) = 0 . (2.9)The slow manifold Z k is precisely the set of equilibria of (2.9) with dim Z k = 5. For thelinearization of (2.9) at each point on Z k j , there are 5 zero eigenvalues associated to the tangentspace of Z k and the other two eigenvalues are ± (cid:112) z c + z c . Thus, Z k is normally hyperbolic([19, 29]). We will denote the stable and unstable manifolds of Z k by W s ( Z k ) and W u ( Z k ),respectively. Let M [ k − , +] be the collection of all forward orbits from B j under the flow of (2.9)and let M [ k, − ] be the collection of all backward orbits from B j , where j , j are two of consecutiveletters in the set { l, a, b, r } (corresponding to B l , B a , B b , B r ) where the place of letters is fixed.Then the set of forward orbits from B j to the corresponding Z k is N [ j , +] = M [ j , +] ∩ W s ( Z k ),and the set of backward orbits from B j to Z k is N [ j , − ] = M [ j , − ] ∩ W u ( Z k ). Therefore, thesingular layer Γ [ j , +] at x j satisfies Γ [ j , +] ⊂ N [ j , +] and the singular layer Γ [ j , − ] at x j satisfiesΓ [ j , − ] ⊂ N [ j , − ] . All the essential geometric objects are explicitly delineated in [40].11he limiting fast (layer) dynamics preserve electrochemical potentials and do not depend ondiffusion constants ([40]). Thus, we directly apply the results on the fast dynamics from [18] andonly point out that we consider n = 2 ion species and need to keep φ a,l , φ a,m , φ b,m and φ b,r in thispaper, while in [18], the equality of diffusion constants D = D concludes φ a,l = V , φ a,m = φ b,m (denoted by V ∗ there) and φ b,r = 0. Lemma 2.3. ([The fast layer dynamics]) One has,a.) over x = a provides, for k = 1 , ,(i) relative to (0 , a ) where Q = 0 , z c a e z ( φ a − φ a,l ) + z c a e z ( φ a − φ a,l ) = 0 , c a,lk = c ak e z k ( φ a − φ a,l ) ; (ii) relative to ( a, b ) where Q (cid:54) = 0 , z c a e z ( φ a − φ a,m ) + z c a e z ( φ a − φ a,m ) + Q = 0 , c a,mk = c ak e z k ( φ a − φ a,m ) ; (iii) the matching u a − = u a + : c a,l + c a,l = c a,m + c a,m + Q ( φ a − φ a,m ) ;b.) over x = b provides, for k = 1 , ,(i) relative to ( a, b ) where Q (cid:54) = 0 , z c b e z ( φ b − φ b,m ) + z c b e z ( φ b − φ b,m ) + Q = 0 , c b,mk = c bk e z k ( φ b − φ b,m ) ; (ii) relative to ( b, where Q = 0 , z c b e z ( φ b − φ b,r ) + z c b e z ( φ b − φ b,r ) = 0 , c b,rk = c bk e z k ( φ b − φ b,r ) ; (iii) the matching u b − = u b + : c b,m + c b,m + Q ( φ b − φ b,m ) = c b,r + c b,r . As stated in the Introduction, the degeneracy of equal diffusion coefficients arises in theslow dynamics. Diffusion is a phenomenon in which the spatial distribution of solute particlesvaries due to their potential energy. It is a standard process that works to remove differences inconcentration and eventually turns a given mixture to a uniform structure state. The diffusion ofuncharged particles can be described by the Fick’s first law [20] through the equation ∂ t c = D ∂ xx c ,where c is the concentration, D is the diffusion constant, and t is time. The measurement ofdiffusion constants frequently involves measuring collections of simultaneous values of t , c , and x .These evaluated values are then implemented to a solution of Fick’s law to achieve the diffusionconstants. There are many approaches and techniques to maintain diffusion constants of ions inaqueous solutions [5, 7, 21, 38].Some kinds of selectivity depend on the non-equality of diffusion coefficients. Besides, manyelectrical phenomena wholly vanish. That is, the liquid coupling is zero when diffusion constantsare identical. Therefore, the equal diffusion constants case is degenerate. On the other hand, ev-erything becomes much more complicated (at least mathematically) when the diffusion constantsare unequal.Some ionic species’ diffusion constants may alter from one technique to another, even when allother parameters are held untouched. However, experimental measurements are directed underisothermal circumstances to keep away from a deviation of D values.12n this part, we briefly recall, with different notation, the analysis on the slow layer fromSection 2 of the paper [44]. For zero current I = z J + z J = 0 (so J = J ), and furtherassumption z = − z = 1, one has J D − J D = D − D D D J . Applying zero current condition in above,the limiting slow system becomes([44]),˙ φ = − ( D − D ) J D D h ( ω )(2 c + Q j ) , ˙ c = − ( D + D ) c + D Q j D D h ( ω )(2 c + Q j ) J , ˙ J = 0 , ˙ w = 1 . (2.10)The slow system (2.10) on (0 , a ) and ( b,
1) where Q = Q = 0 and on ( a, b ) where Q (cid:54) = 0 willresult in the following Lemma. Lemma 2.4.
The slow dynamics over each interval (slow manifold) is as follows:(a) Over (0 , a ) with c ( x ) − c ( x ) = 0 the slow dynamics system gives, c a,l = l − D + D D D J H ( a ) , φ a,l = V − D − D D + D ln c a,l l . (b) Over ( a, b ) with c ( x ) − c ( x ) + Q = 0 the slow dynamics system gives, φ b,m = φ a,m + D − D D D J y ∗ ,c b,m = e − D D D D J y ∗ c a,m + D Q D + D (cid:0) e − D D D D J y ∗ − (cid:1) ,J = − D D (cid:0) c b,m − c a,m (cid:1) − ( φ b,m − φ a,m ) Q ( D + D ) (cid:0) H ( b ) − H ( a ) (cid:1) . (c) Over ( b, with c ( x ) − c ( x ) = 0 the slow dynamics system gives, c b,r = r + D + D D D J (cid:0) H (1) − H ( b ) (cid:1) , φ b,r = D − D D + D ln r c b,r . Matching for Zero-current and Singular Orbits on [0 , . The last step to build a connecting orbit over the whole interval [0 ,
1] is to match the three singularorbits from the Lemmas 2.3 and 2.4 at the points x = a and x = b . The matching conditions are u a − = u a + , u b − = u b + , and that J (= J ) needs to be equal on all three subintervals. Thus, c a e ( φ a − φ a,l ) − c a e − ( φ a − φ a,l ) = c b e ( φ b − φ b,r ) − c b e − ( φ b − φ b,r ) = 0 ,c a e ( φ a − φ a,m ) − c a e − ( φ a − φ a,m ) + Q = 0 ,c b e ( φ b − φ b,m ) − c b e − ( φ b − φ b,m ) + Q = 0 , c a,l = c a e ( φ a − φ a,m ) + c a e − ( φ a − φ a,m ) + Q ( φ a − φ a,m ) , c b,r = c b e ( φ b − φ b,m ) + c b e − ( φ b − φ b,m ) + Q ( φ b − φ b,m ) ,J = J = − D D ( c a,l − l )( D + D ) H ( a ) = − D D ( r − c b,r )( D + D )( H (1) − H ( b )) , = − D D c b,m − c a,m ) − ( φ b,m − φ a,m ) Q ( D + D )( H ( b ) − H ( a )) ,φ b,m = φ a,m + D − D D D J y ∗ ,c b,m = e − D D D D J y ∗ c a,m + D Q ( D + D ) (cid:16) e − D D D D J y ∗ − (cid:17) . (2.11)13 emark . In (2.11), the unknowns are: φ a , φ b , c a , c a , c b , c b , J , φ a,m , φ b,m , y ∗ and Q thatis, there are eleven unknowns that matches the total number of equations on (2.11). This section investigates how some quantities, such as boundary concentrations and electric po-tential, work to make the current reverse. We analyze the outcomes for various diffusion constantsto reveal the impacts of diffusion constants on the geometric mean of concentrations and reversalpermanent charge and expose the opposites.We encountered a nonlinear governing system in (2.11) that is difficult to analyze as it is. In[18], for equal diffusion constants, D k ’s, the investigation on reversal permanent charges imple-mented for a general n . However, with the general D k ’s, the difficulty intensifies, as shown in[44], even for the case we will address in this work, where there are two ion species, i.e., n = 2with z = − z = 1.In [16] and [44], without and with zero current assumptions, respectively, the authors offeredtwo intermediate variables that support a notable reduction of the matching (2.11): A = (cid:112) c ( a ) c ( a ) , B = (cid:112) c ( b ) c ( b ) . (3.1)The variables A and B are the geometric mean of concentrations at x = a and x = b respectively.One may consider B as a function of A rather than an independent variable. In fact, it can beseen that B = B ( A ) = 1 − βα ( l − A ) + r . We can now reduce the Matching system (2.11) to anonlinear system with two equations and two unknowns. We skip the redued step here, for zerocurrent I = J − J = 0, as illustrated in detail in the Appendix of [44]. G ( Q , A, δ ) = V and G ( Q , A, δ ) = 0 , (3.2)where G ( Q , A, δ ) = δ (cid:16) ln S a + δQ S b + δQ + ln lr (cid:17) − (1 + δ ) ln AB + ln S a − Q S b − Q ,G ( Q , A, δ ) = δQ ln S a + δQ S b + δQ − N. (3.3)Furthermore, S a = (cid:113) Q + A , S b = (cid:113) Q + B , N = A − l + S a − S b , (3.4)and δ = D − D D + D , H ( x ) = (cid:90) x h ( s ) ds, α = H ( a ) H (1) , β = H ( b ) H (1) . (3.5)The function H ( x ) is the ratio of the length with the cross-section area of the portion of thechannel over [0 , x ]. The quantity H ( x ) initially has roots in Ohm law for a uniform resistor’sresistance. The quantities α and β , together with Q , are critical characteristics for the shapeand permanent charge of the channel structure [34].14artial derivatives of G and G with respect to Q and A are, ∂ A G ( A, Q , δ ) =(1 − δ ) Q (cid:16) A ( S a + δQ ) + 1 − βα B ( S b + δQ ) (cid:17) ,∂ Q G ( A, Q , δ ) = (1 − δ )( S a − S b )( S a + δQ )( S b + δQ ) ,∂ δ G ( A, Q , δ ) = ln S a + δQ S b + δQ + δQ S a + δQ − δQ S b + δQ + ln lr − ln AB ,∂ A G ( A, Q , δ ) = − − βα BS b + δQ − AS a + δQ − β − αα ,∂ Q G ( A, Q , δ ) = δ ln S a + δQ S b + δQ + (1 − δ ) Q ( S a − S b )( S a + δQ )( S b + δQ ) ,∂ δ G ( A, Q , δ ) = Q ln S a + δQ S b + δQ + δQ S a + δQ − δQ S b + δQ . (3.6)It follows directly from (3.6) that G ( Q , A, δ ) is increasing (decreasing) in A as Q is increasing(decreasing), and it is increasing (decreasing) in Q or δ when l − r increases (decreases). Moreover, G always decreases in A regardless of other values; however, its behavior with respect to Q and δ is complicated. The preceding observations have been confirmed in Lemma 3.2 of [44] in detail;however, the authors of that paper did not address how G acts (with respect to Q or δ ) when δQ <
0. We develop and complete the Lemma in the following Theorem to serve us in the nexttracks. In particular, it will help us later in Section 3.2, to establish the uniqueness of reversalpotential.
Theorem 3.1.
Assume that diffusion constants D and D are fixed:(a) if D < D , i.e. δ > , then there exists a Q − < , so that ∂ Q G ( Q , A, δ ) has the same signas that of l − r when Q > Q − , but it has the opposite sign when Q < Q − , and(b) if D > D , i.e. δ < , then there exists a Q +0 > , so that ∂ Q G ( Q , A, δ ) has the oppositesign as that of l − r when Q < Q +0 , but it has the same sign when Q > Q +0 .Proof. We prove (a). The other statement is similar. It follows from (3.3) that ∂ Q G for small(large) values of Q < l − r . And for any values of Q >
0, it always has the same sign as that of l − r . On the other hand, taking one morederivative, ∂ Q G can be written (after making common factors) in form of P ( Q ) P ( Q ) where P and P are two polynomials (in Q ) of degrees 2 and 4 respectively. It is easy to see that P ( Q ) (cid:54) = 0for any Q , and P ( Q ) has at most two roots. Thus, ∂ Q G may change sign at most twice. Thiscompletes the proof. Corollary 3.2.
A direct conclusion of Theorem 3.1 is that Q − and Q +0 take, respectively, theabsolute minimum and maximum of ∂ Q G for each case.Remark . We will see later in Theorem 3.4 that ∂ Q G and ∂ Q A have the same behavior.Hence, Figure 3 for A may also be trusted to help the reader understand Theorem 3.1. We established the geometric mean of concentration on (3.1) and stated why we need to specifythis parameter. The authors of [44] prove that, for fixed Q , D and D , one can solve for15 from G ( Q , A, δ ) = 0, where G is defined in (3.3). Moreover, A (0 , δ ) = (1 − α ) l + αr andlim Q →±∞ A ( Q , δ ) = l . Then, they showed that l, r, A and B satisfy one the following conditions: l < A ( Q , δ ) < B ( Q , δ ) < r or l > A ( Q , δ ) > B ( Q , δ ) > r . If δQ ≥
0, then ∂ Q A ( Q , δ ) hasthe same sign as that of ( l − r ) Q .We now present a few more features of the geometric mean of concentrations, A , through thefollowing Theorems and a Corollary. The same results can be achieved for the other parameter, B . The following Lemma is directly obtained from the equations (3.2)–(3.5). Lemma 3.3.
For any given Q ∈ R and for any A ∈ (cid:0) , A M (cid:1) where A M = l + α − β r is themaximum value of A , the functions G ( Q , A, δ ) and G ( Q , A, δ ) in (3.3) are continuous. Fur-thermore, ∂ A N > , ∂ A B < .Remark . It might be deduced, incorrectly from Corollary 3.2 and Lemma 3.3, that A M is equalto A ( Q ) evaluated at Q = Q − , or Q = Q + . However, the function A ( Q ) is not surjectivenecessarily, and in general, A ( Q ± ) ≤ A M . In next Theorem, we discuss about the bounds of thegeometric mean of concentration in more details. Theorem 3.4.
Assume that diffusion constants D and D are fixed. Let Q − < , Q +0 > be asdefined in Theorem 3.1. Subsequently,(a) if D < D , i.e. δ > , then ∂ Q A ( Q , δ ) has the same sign as that of l − r when Q > Q − ,but it has the opposite sign when Q < Q − , and(b) if D > D , i.e. δ < , then ∂ Q A ( Q , δ ) has the opposite sign as that of l − r when Q < Q +0 ,but it has the same sign when Q > Q +0 .Proof. This is the direct conclusion of Theorem 3.1 with the facts that ∂ Q A = − ∂ Q G /∂ A G and ∂ A G < Q . Remark . Note that Q − and Q +0 in Theorems 3.1 and 3.4 are the same. Similarly, we emphasizethat the same conclusions can be established for the case with D > D . Theorem 3.5.
For any given Q = Q there is a unique A ( Q ) such that G ( A ( Q ) , Q ) = 0 .Furthermore, lim Q →±∞ A ( Q ) = l , and A = A ( Q ) satisfies the following,(a) if l < r then l ≤ A ( Q ) ≤ A ( Q ± ) , (b) if l > r then A ( Q ± ) ≤ A ( Q ) ≤ l ,(c) if l = r then A ( Q ) = l = r = A ( Q ± ) ,where A ( Q ± ) = A ( Q + ) or A ( Q ± ) = A ( Q − ) depending on the sign of diffusion constants; also, Q +0 , and Q − were defined in Theorem 3.4.Proof. The uniqueness of A = A ( Q ) has been proven in [44]. To show the limit, since lim Q →±∞ S a ( Q ) = S b ( Q ) , it follows from G ( A ( Q ) , Q ) = 0 that, δ lim Q →±∞ Q ln S a + pQ S b + pQ = β − αα lim Q →±∞ ( A − l ) . On the other hand,lim Q →±∞ Q ln S a + pQ S b + pQ = − lim Q →±∞ ( Q S a + p )( S b + pQ ) − ( Q S b + p )( S a + pQ )( S a + pQ )( S b + pQ ) Q = 0 . Thus, lim Q →±∞ A ( Q ) = l . The parts (a)–(c) concludes from and Theorem 3.4.16 orollary 3.6. The geometric mean of concentrations A = A ( Q ) is always finite and boundedbetween l and A ( Q ± ) for any values of Q . Similarly, the other geometric mean of concentrations B ( Q ) is unique and bounded between r and B ( Q ± ) . Besides, lim Q →±∞ B ( Q ) = r. Furthermore,if A is surjective, then A ( Q ± = A M where A M was defined in Lemma 3.3; and B ( Q ± = B M correspondingly. Recall that the quantities B = B ( A ) and N = N ( A ) defined in (3.1) and (3.4). One cansimply see that B − A is decreasing over domain of A . Besides, B − A > A → + (consequently B > B − A < A → A − M (consequently B → + ). Therefore,there exists some A ∗ ∈ (0 , A M ) for which B ∗ := B ( A ∗ ) = A ∗ . Moreover, Lemma 3.3 along with N (0 + ) < , N ( A − M ) > A ( Q ) ∈ (0 , A M ) such that N ( ˆ A ) = 0. Thecorresponding quantities ˆ B and B ∗ are similarly defined. Lemma 3.7.
For any given Q one has,(i) If l < r then l < ˆ A < A ∗ = B ∗ < ˆ B < r , (ii) If l > r then l > ˆ A > A ∗ = B ∗ > ˆ B > r ,(iii) If l = r then ˆ A = ˆ B = A ∗ = B ∗ = l = r .Proof. Suppose r < l . Since N ( A = l ) > , N ( A = ˆ A ) = 0 and N is increasing with respect to A then we get ˆ A < l . Now, it follows from definition of B , N ( ˆ A ) = 0 and ˆ A < l that r < ˆ B < ˆ A .Now, set f ( A ) := B − A. Then, f ( ˆ A ) = ˆ B − ˆ A < f ( ˆ B ) = (cid:16) − βα ( l − ˆ B ) + r (cid:17) − (cid:16) − βα ( l − ˆ A ) + r (cid:17) = 1 − βα ( ˆ A − ˆ B ) > . Thus, it follows from ∂ A f < f ( A ∗ ) = 0 that ˆ B < A ∗ = B ∗ < ˆ A . The other cases aresimilar. Theorem 3.8.
For equal diffusion constant D = D , the solution A ( Q , of G ( A ( Q ) , Q ,
0) =0 is symmetric with respect to Q ; but for any ( Q , δ ) where δ (cid:54) = 0 , i.e. D (cid:54) = D , the function A ( Q , δ ) is non-symmetric with respect to Q .Proof. It directly follows from G = 0 in (3.2) and (3.3).Figures 3 and 4 confirms the discussions in Theorems 3.4, 3.5 and 3.8. In what follows,numerical simulations are conducted with the help of analysis on system (3.2). The combination ofnumerics and analysis gives a better understanding of the zero-current problems and complimentssome analytical results obtained in [44]. For our numerical simulations, we assume that h ( x ) = k for any x ∈ [ a, b ] with a = 1 / , b = 2 / < k <
1. Wefurther suppose that h (0) = h (1) = 1, and assume h ( x ) is approximated by two linear functionsover the non-constant intervals (0 , /
3) and (2 / , h ( x ) is defined as a piece-wise linearfunction over [0 , H ( x ) in (3.5) that H (0) = 0 , H ( a ) = a ln kk − , H ( b ) = a ln kk − b − ak , H (1) = (1 − b + a ) ln kk − b − ak . Thus, for 0 < k < α = α ( a, b, k ) = ak ln k (1 − b + a ) k ln k + ( k − b − a ) ,β = β ( a, b, k ) = ak ln k + ( k − b − a )(1 − b + a ) k ln k + ( k − b − a ) . (3.7)Since 0 < a < b <
1, it follows from above that 0 < α < β . The following Theorem is the directconclusion of (3.7). 17 heorem 3.9.
One has,(i) if a → + , then α → + , β → ( k − b (1 − b ) k ln k + ( k − b . For this case, k → + if and only if β → − , and k → − if and only if β → b ;(ii) if b → − , then β → − , α → ak ln kk − . For this case, k → + if and only if α → + , and k → − if and only if α → a .Furthermore, when a and b are close, one obtains that all values α, β, a and b are close too. -5 0 50.110.120.130.140.150.160.17 NaCl Na CO -5 0 50.230.240.250.260.270.280.29 NaCl Na CO Figure 3:
The function A ( Q , δ ) for two pairs of ( D , D ) corresponding to N aCl and
N a CO ,for various values of Q : left panel for L = 0 . M, R = 0 . M ; right panel for L = 0 . M, R = 0 . M . -30 -20 -10 0 10 20 3011.11.21.31.41.51.61.7 NaCl Na CO -30 -20 -10 0 10 20 302.32.42.52.62.72.82.93 NaCl Na CO Figure 4:
The function A ( Q , δ ) for two pairs of ( D , D ) corresponding to N aCl and
N a CO for various values of V : left panel for L = 0 . M, R = 0 . M ; right panel for L = 0 . M, R = 0 . M . .2 Results on Reversal Permanent Charge Q rev . The spatial distribution of side chains in a specific channel defines the permanent charge of thechannel, which forms most of the electrical structure of the channel protein. Thus, reversalpotentials should always exist within the ion channels. However, there is a simple necessarycondition for the existence of the reversal permanent charge Q rev , as one will see in the nextTheorem. The general result for reversal permanent charge with a given electric potential V is asfollows. Theorem 3.10.
For n = 2 , there exists a unique reversal permanent charge Q rev if and only if (cid:12)(cid:12) V (cid:12)(cid:12) < (cid:12)(cid:12)(cid:12) ln lr (cid:12)(cid:12)(cid:12) . (3.8) Proof.
Existence of permanent charge Q rev have been proved in [44]. The uniqueness is theconsequence of Theorems 3.4, 3.5 and ∂ Q G in (3.6). Theorem 3.11.
For any given ( V , l, r ) that satisfies the condition (3.8) one has,(a) if l < r , then ∂ V Q rev < and lim V →± ln lr Q rev ( V ) = ∓∞ ;(b) if l > r , then ∂ V Q rev > and lim V →± ln lr Q rev ( V ) = ±∞ Proof.
For any given V from G (cid:0) V , Q rev ( V ) , A ( Q rev ( V )) (cid:1) = 0 we have, − ∂ Q G ∂ V Q rev + ∂ A G ∂ Q A∂ V Q rev = 0 . It follows from above and ∂ Q A = − ∂ Q G ∂ A G that, ∂ V Q rev = ∂ A G ∂ Q G ∂ A G − ∂ A G ∂ Q G . (3.9)The parts (a) and (b)statements can be verified from (3.9), and Theorem 3.4.The numerical investigations in Figure 5 admits our results in Theorem 3.11. It shows thegraph of Q rev for values of V where − ln lr ≤ V ≤ ln lr . One can see that when V go to ± ln lr then Q rev becomes large unbounded. In this manuscript, we work on the classical PNP model allowing unequal diffusion constants andfor a single profile of permanent charges, to study the specific questions about reversal potentialsand reversal permanent charges that are among the central issues of biological concerns.A crucial assumption is that the dimensionless parameter ε of the ratio of the Debye lengthover the distance between the two applied electrodes is small. The assumption allows one to treatthe PNP system as a singularly perturbed system with ε as the singular parameter.Our study relies on a modern general geometric singular perturbation theory and some uniquestructures of the classical PNP models. Then we obtain a nonlinear matching system of algebraicequations (2.11) for the zero current condition that includes both the reversal potential andreversal permanent charge topics. We use an intermediate variable introduced in [16], to furtherreduce the matching system to an effective system of two algebraic equations with two unknowns.19
20 -10 0 10 20-200-150-100-50050100150200
NaCl Na CO -20 -10 0 10 20-200-150-100-50050100150200 NaCl Na CO Figure 5:
The function Q rev for two pairs of ( D , D ) corresponding to N aCl and
N a CO : leftpanel for L = 0 . M < R = 0 . M ; right panel for L = 0 . M > R = 0 . M . Several novel properties of biological significance have resulted from the analysis of these governingequations that some are not intuitive. In the future, we intend to maintain the study to extendthe analysis in this work and explore the problems numerically applying more advanced andcomplex models. More numerical observations can also be presented to study profiles of relevantphysical quantities, e.g., to numerically investigate the behavior of C k ( X ), Φ( X ) and µ k ( X )throughout the channel. There are many other inspiring projects one can initiate related tothis work. Another exciting novel project combines ion size to the problem using hard-sphereelectrochemical potentials to analyze ion size’s effects. Acknowledgement.
The author thanks Dr. Weishi Liu for beginning such an exciting topic,for his support and advice.
References [1] B. Balu and A. Khair, Role of Stefan-Maxwell fluxes in the dynamics of concentrated electrolytes.
Soft Matter. (2018), 8267-8275.[2] M.Z. Bazant, M.S. Kilic, B.D. Storey, and A. Ajdari, Towards an understanding of induced chargeelectrokinetics at large applied voltages in concentrated solutions. Adv. Coll. Interf. Sci. (2009),48-88.[3] V. Barcilon, Ion flow through narrow membrane channels: Part I.
SIAM J. Appl. Math. (1992),1391-1404.[4] V. Barcilon, D.-P. Chen, and R. S. Eisenberg, Ion flow through narrow membrane channels: Part II. SIAM J. Appl. Math. (1992), 1405-1425.[5] A. J. Bard and L.R. Faulkner, Electrochemical Methods, Fundamentals and Applications.
Wiley, NewYork, 1980.[6] D. Boda, W. Nonner, M. Valisko, D. Henderson, B. Eisenberg, and D. Gillespie, Steric selectivityin Na channels arising from protein polarization and mobile side chains.
Biophys. J. (2007),1960-1980.[7] R.E. Brooks, L.O. Heflinger, and R. F. Wuerker, Interferometry with a holographically reconstructedcomparison beam. J. Applied Physics letters (1965) 248-249.[8] M. Burger, H. W. Engl, P. A. Markowich, and P. Pietra, Identification of doping profiles in semicon-ductor devices. Inverse Problems (2001), 1765–1795.
9] D. P. Chen and R.S. Eisenberg, Charges, currents and potentials in ionic channels of one conforma-tion.
Biophys. J. (1993), 1405-1421.[10] D. Chen, R. Eisenberg, J. Jerome and C. Shu, Hydrodynamic model of temperature change in openionic channels. Biophysical J. (1995), 2304-2322.[11] T. F. Coleman and Y. Li, An interior, trust region approach for nonlinear minimization subject tobounds. SIAM J. Optimization (1996), 418-445.[12] S. Chung and S. Kuyucak, Predicting channel function from channel structure using Brownian dy-namics simulations. Clin. Exp. Pharmacol Physiol. (2001), 89-94.[13] B. Eisenberg, Ion Channels as Devices. J. Comp. Electro. (2003), 245-249.[14] B. Eisenberg, Crowded charges in ion channels. Advances in Chemical Physics (ed. S. A. Rice) (2011),77-223, John Wiley and Sons, Inc. New York.[15] B. Eisenberg, Y. Hyon, and C. Liu, Energy variational analysis of ions in water and channels: Fieldtheory for primitive models of complex ionic fluids.
J. Chem. Phys. (2010), 104104(1-23).[16] B. Eisenberg and W. Liu, Poisson-Nernst-Planck systems for ion channels with permanent charges.
SIAM J. Math. Anal. (2007), 1932-1966.[17] B. Eisenberg and W. Liu, Relative dielectric constants and selectivity ratios in open ionic channels. Mol. Based Math. Biol. (2017), 125-137.[18] B. Eisenberg, W. Liu, and H. Xu, Reversal permanent charge and reversal potential: case studies viaclassical Poisson-Nernst-Planck models. Nonlinearity (2015), 103-128.[19] N. Fenichel, Geometric singular perturbation theory for ordinary differential equations. J. DifferentialEquations (1979), 53-98.[20] A. Fick, On liquid diffusion. Philos. Mag. J. Sci. (1855), 3139.[21] G. Gerhardt, R.N. Adams, Determination of diffusion constants by flow injection analysis. Analyticalchemistry J. (1982), 2618-2620.[22] D. Gillespie, W. Nonner, and R. S. Eisenberg, Coupling Poisson-Nernst-Planck and density functionaltheory to calculate ion flux. J. Phys.: Condens. Matter (2002), 12129-12145.[23] S. M. H. Hashemi Amrei, S. C. Bukosky, S. P. Rader, W. D. Ristenpart, and G. H. Miller, OscillatingElectric Fields in Liquids Create a Long-Range Steady Field. Phys. Rev. Lett.
J. Math. Biol. (2010),347-386.[25] B. Hille, Ion Channels of Excitable Membranes (Third Edition).
Sinauer Associates, Inc., Sunderland,Massachusetts, USA 2001.[26] B. Hille, Transport Across Cell Membranes: Carrier Mechanisms, Chapter 2. Textbook of Physiology(ed. H. D. Patton, A. F. Fuchs, B. Hille, A. M. Scher and R. D. Steiner). Philadelphia, Saunders (1989), 24-47.[27] Y. Hyon, B. Eisenberg, and C. Liu, A mathematical model for the hard sphere repulsion in ionicsolutions. Commun. Math. Sci. (2010), 459-475.[28] Y. Hyon, J. Fonseca, B. Eisenberg, and C. Liu, Energy variational approach to study charge inversion(layering) near charged walls. Discrete Cont. Dyn. Syst. B (2012), 2725-2743.[29] M. Hirsch, C. Pugh, and M. Shub, Invariant Manifolds.
Lecture Notes in Math. , Springer-Verlag,New York, 1976.[30] U. Hollerbach, D.-P. Chen, and R.S. Eisenberg, Two- and Three-Dimensional Poisson-Nernst-PlanckSimulations of Current Flow through Gramicidin-A.
J. Comp. Science (2001), 373-409.
31] W. Im and B. Roux, Ion permeation and selectivity of OmpF porin: a theoretical study based onmolecular dynamics, Brownian dynamics, and continuum electrodiffusion theory.
J. Mol. Biol. (2002), 851-869.[32] S. Ji, B. Eisenberg, and W. Liu, Flux Ratios and Channel Structures.
J. Dynam. Differential Equa-tions (2017). https://doi.org/10.1007/s10884-017-9607-1.[33] S. Ji and W. Liu, Poisson-Nernst-Planck systems for ion flow with density functional theory forhard-sphere potential: I-V relations and critical potentials. Part I: Analysis.
J. Dynam. DifferentialEquations (2012), 955-983.[34] S. Ji, W. Liu, and M. Zhang, Effects of (small) permanent charge and channel geometry on ionicflows via classical Poisson-Nernst-Planck models, SIAM J. Appl. Math. (2015), 114-135.[35] C. Jones, Geometric singular perturbation theory, in: dynamical systems, montecatini terme, 1994,in: Lecture Notes in Math. Springer-Verlag, Berlin, (1995), 44-118.[36] C. Jones and N. Kopell, Tracking invariant manifolds with differential forms in singularly perturbedsystems.
J. Differential Equations (1994), 64-88.[37] C. Kuehn. Multiple Time Scale Dynamics, Springer, 2015.[38] J.-L. Liu and B. Eisenberg, Poisson–Nernst–Planck-Fermi theory for modeling biological ion channels,
J. Chem. Phys. (2014), 12B640.[39] W. Liu, Geometric singular perturbation approach to steady-state Poisson-Nernst-Planck systems.
SIAM J. Appl. Math. (2005), 754-766.[40] W. Liu, One-dimensional steady-state Poisson-Nernst-Planck systems for ion channels with multipleion species. J. Differential Equations (2009), 428-451.[41] W. Liu and B. Wang, Poisson-Nernst-Planck systems for narrow tubular-like membrane channels.
J.Dynam. Differential Equations (2010), 413-437.[42] W. Liu and H. Xu, A complete analysis of a classical Poisson-Nernst-Planck model for ionic flow. J.Differential Equations (2015), 1192-1228.[43] W. Liu, A flux ratio and a universal property of permanent charges effects on fluxes.
Comput. Math.Biophys. (2018), 28-40.[44] H. Mofidi and W. Liu, Reversal potential and reversal permanent charge with unequal diffusioncoefficients via classical Poisson–Nernst–Planck models. Submitted.
Preprint available at arXiv as arXiv:1909.01192. [45] H. Mofidi, B. Eisenberg and W. Liu, Effects of Diffusion Coefficients and Permanent Charge onReversal Potentials in Ionic Channels.
Entropy (2020), 325(1-23).[46] W. Nonner and R. S. Eisenberg, Ion permeation and glutamate residues linked by Poisson-Nernst-Planck theory in L-type Calcium channels. Biophys. J. (1998), 1287-1305.[47] W. Nonner, L. Catacuzzeno, and B. Eisenberg, Binding and selectivity in L-type Calcium channels:A mean spherical approximation. Biophys. J. (2000), 1976-1992.[48] J.-K. Park and J. W. Jerome, Qualitative properties of steady-state Poisson-Nernst-Planck systems:Mathematical study. SIAM J. Appl. Math. (1997), 609-630.[49] D. J. Rouston, Bipolar Semiconductor Devices.
McGraw-Hill, New York, 1990.[50] V. Sasidhar and E. Ruckenstein, Electrolyte osmosis through capillaries.
J. Colloid Interface Sci. (1981), 439-457.[51] Z. Schuss, B. Nadler, and R. S. Eisenberg, Derivation of Poisson and Nernst-Planck equations in abath and channel from a molecular model. Phys. Rev. E (2001), 1-14.
52] L. Sun and W. Liu, Non-localness of excess potentials and boundary value problems of Poisson-Nernst-Planck systems for ionic flow: A Case Study.
J. Dynam. Differential Equations (2018),779-797.[53] P. Vanysek, Ionic conductivity and diffusion at infinite dilution. Handbook of Chemistry and Physics,74th ed.; Linde, D. R., Ed.; CRC: Boca Raton, (1993), 90-92.[54] R. M. Warner, Jr., Microelectronics: Its unusual origin and personality. IEEE Trans. Electron.Devices (2001), 2457-2467.[55] L. Zhang, B. Eisenberg, and W. Liu, An effect of large permanent charge: Decreasing flux withincreasing transmembrane potential. Eur. Phys. J. Special Topics (2019), 2575-2601.(2019), 2575-2601.