Reversal potential and reversal permanent charge with unequal diffusion coefficients via classical Poisson-Nernst-Planck models
RReversal potential and reversal permanent charge with unequaldiffusion coefficients via classical Poisson-Nernst-Planck models
Hamid Mofidi ∗ and Weishi Liu † September 4, 2019
Abstract
In this paper, based on geometric singular perturbation analyses of a quasi-one dimen-sional Poisson-Nernst-Planck model for ionic flows, we study the problem of zero currentcondition for ionic flows through membrane channels with a simple profile of permanentcharges. For ionic mixtures of multiple ion species, under equal diffusion constant condition,Eisenberg, et al [
Nonlinearity (2015), 103-128 ] derived a system of two equations fordetermining the reversal potential and an equation for the reversal permanent charge. Theequal diffusion constant condition is significantly degenerate from physical points of view.For unequal diffusion coefficients, the analysis becomes extremely challenging. This work willfocus only on two ion species, one positively charged (cation) and one negatively charged(anion), with two arbitrary diffusion coefficients. Mathematically, we identify two governingequations for the zero current which, for the first time, allow one to examine how the reversalpotential depends on the channel structure and diffusion coefficients; In particular, we areable to show, with a number of concrete results, that the possible different diffusion constantsindeed make significant differences. The inclusion of channel structures is also far beyondthe situation where the Goldman-Hodgkin-Katz (GHK) equation might be applicable. Acomparison of our result with the GHK equation is provided. The dual problem of reversalpermanent charges is briefly discussed too. Key words.
GSP for PNP, reversal potential, reversal permanent charge
Ion channels, proteins embedded in membranes, provide a major channel for cells to communicatewith each other and with the outside to transform signals and to conduct group tasks ([5, 9, 15,16]). The key structure of an ion channel is its shape and its permanent charge. The shape ofa 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 contributing positivecharges. It is the specific of side chain distributions that is referred to as the permanent chargeof the ion channel. The function of channel structures is to select the types of ions and tofacilitate the diffusion of ions across cell membranes. ∗ Department of Mathematics, University of Kansas, 1460 Jayhawk Blvd., Room 405, Lawrence, Kansas 66045,USA ( [email protected] ). † Department of Mathematics, University of Kansas, 1460 Jayhawk Blvd., Room 405, Lawrence, Kansas 66045,USA ( [email protected] ). a r X i v : . [ m a t h . D S ] S e p t present, these permeation and selectivity properties of an ion channel are mainly extractedfrom the I-V relation measured experimentally ([15, 16, 20, 21, 22, 23]). Individual fluxes carrymore information than the I-V relation but it is expensive and challenging to measure them([24, 27]). The I-V relation is a functional response of the channel structure on ionic fluxes butit depends on boundary conditions that are in fact driving forces of ionic transport. The multi-scale feature of the problem with multiple physical parameters allow the system to have a greatflexibility and to exhibit rich phenomena/behaviors – a great advantage of “natural devices” ([8]).On the other hand, the same multi-scale feature with multiple physical parameters presents anextremely challenging task for anyone to extract meaningful information from experimental data,given also the fact that the internal dynamics cannot be observed with present technique.Mathematical analysis plays important and unique roles for explaining mechanisms of ob-served biological phenomena and for discovering new ones, assuming a more or less explicitsolution of the associated mathematical model can be obtained. The latter is often too much tohope. Nonetheless, there have been some successes recently in analyzing Poisson-Nernst-Planck(PNP) models for ionic flows through ion channels ([10, 11, 27, 29, 31, 32, 35, 39], etc.).In this work, we are interested in reversal potentials (or Nernst potentials) as well as reversalpermanent charges . They are defined by zero total current: for fixed other physical quantities,the total current I = I ( V, Q ) depends on the transmembrane potential V and the permanentcharge Q . For fixed Q , a reversal potential V = V rev ( Q ) is a transmembrane potential thatproduces zero current I ( V rev ( Q ) , Q ) = 0. Likewise, for fixed transmembrane potential V , areversal permanent charge Q = Q rev ( V ) is a permanent charge that produces zero current I ( V, Q rev ( V )) = 0.Nernst was among the first who considered reversal potential and, for one ion species case,formulated an equation – now called the Nernst equation – for the reversal potential. Followinga treatment of Mott for electronic conduction in the copper-copper oxide rectifier ([36]), theNernst equation was generalized by Goldman ([14]), and Hodgkin and Katz ([19]) – called
Goldman-Hodgkin-Katz (GHK) equation - for reversal potentials involving multiple ion species.The derivations were based on the assumption that the electric potential φ ( x ) is linear in x –the coordinate along the longitude of the channel. The assumption was known wrong by thoseauthors and, unfortunately, there was no substitute yet for their equations.Recently in [11], the authors investigated the problem of determining reversal potentials andreversal permanent charges based on rigorous analysis on the Poisson-Nernst-Planck models.For the case when all diffusion constants are equal, the results are very complete. In particular,a system of two equations is derived that will lead to a determination of the reversal potential,and one equation is derived for the reversal permanent charge. On the other hand, the equaldiffusion constants case is quite degenerate, which is known from biological point of view evenfor ionic mixtures of two ion species. In this work, allowing different diffusion coefficients, westart our investigation on reversal potentials and reversal permanent charges for two ion species.We are particularly interested in the effect of ρ = D /D on the values of reversal potentialsand reversal permanent charges, where D k is diffusion constant for the k th ion species.The geometric singular perturbation framework developed in [10, 31, 32] particularly foranalyzing PNP models for ionic flow is again applied as in [11] to get a system of algebraicequations for the problem. The solution method of solving/analyzing the algebraic system issimply different from that in [11] due to the difference between D and D . The difficultyis overwhelmingly increased – more than technical. An important step in our analysis is areduction of the algebraic system to two nonlinear equations that turns out to work effectively.2s a consequence, this reduced system allows one to, for the first time, examine how the reversalpotential depends on the channel structure, boundary concentrations and diffusion coefficients.In particular, we are able to establish a number of precise differences that possible differentdiffusion constants make. Some of these results can be explained qualitatively in terms of physicalintuitions, for examples, the dependence of the sign of reversal potential on interplay betweendiffusion constants, boundary conditions and permanent charge (Theorem 4.2 and Corollary 4.3),how the monotonicity of the reversal potential in the permanent charge depends on the relativesizes of the diffusion constants together with the boundary conditions (Theorem 4.4), etc. Someare counterintuitive, including the specific dependence on the boundary concentrations of themonotonicity of the reversal potential in ρ (Proposition 4.6 and Remark 4.2). All these resultsare not known before and there are also several concrete open questions that we share our beliefbut could not verify. The well-known GHK equation for the reversal potential is briefly discussedand a short comparison with our result is provided.The rest of paper is divided as follows. In Section 1.1 we introduce the problem and providethe basic setup for our problem in Section 1.2. We apply the geometric singular perturbationtheory in Section 2 to derive the matching system of algebraic equations for the zero currentcondition. In Section 3, we discuss the reduced system for a simpler case and make preparationfor our main concern. The topics on reversal potential, its existence, uniqueness and dependenceon permanent charge and diffusion coefficients, are analyzed in Section 4. The topic on reversalpermanent charge is briefly discussed in Section 5. Section 6 is a short conclusion. The appendix(Section 7) details the reduction to the system of two equations for the zero current. The PNP system has been extensively studied by simulations and computations ([1, 2, 3, 4,6, 7, 25, 26]). It is clear from these simulations that macroscopic reservoirs – mathematicallyboundary conditions – must be included in the mathematical formulation to describe the actualbehavior of channels ([13, 38]). On the basis that ion channels have narrow cross-sectionsrelative to their lengths, 3-D PNP type models are further reduced to quasi-one-dimensionalmodels ([34, 37]): 1 A ( X ) ddX (cid:16) ε r ( X ) ε A ( X ) ddX Φ (cid:17) = − e (cid:16) n (cid:88) s =1 z s C s + Q ( X ) (cid:17) ,ddX J k = 0 , −J k = 1 k B T D k A ( X ) C k ddX µ k , k = 1 , , · · · , n (1.1)where X ∈ [0 , l ] is the coordinate along the longitudinal axis of the channel, A ( X ) is the area ofcross-section of the channel over the location X ; Q ( X ) is the permanent charge density, ε r ( X )is the relative dielectric coefficient, ε is the vacuum permittivity, e is the elementary charge, k B is the Boltzmann constant, T is the absolute temperature; Φ is the electric potential, and,for the k th ion species, C k is the concentration, z k is the valence (the number of charges perparticle), µ k is the electrochemical potential depending on Φ and C k , J k ( X ) is the flux densitythrough the cross-section over X , and D k is the diffusion coefficient.Equipped with system (1.1), we impose the following boundary conditions, for k = 1 , , · · · , n ,Φ(0) = V , C k (0) = L k >
0; Φ( l ) = 0 , C k ( l ) = R k > . (1.2)3or an analysis of the boundary value problem (1.1) and (1.2), we will work on a dimensionlessform. Let C be a characteristic concentration of the problems, for example, C = max ≤ k ≤ n (cid:8) L k , R k , sup X ∈ [0 ,l ] |Q ( X ) | (cid:9) . Set D = max ≤ k ≤ n { sup X ∈ [0 ,l ] D k ( X ) } and ¯ ε r = sup X ∈ [0 ,l ] ε r ( X ) . Let ε = ¯ ε r ε k B Te l C , ˆ ε r ( x ) = ε r ( X )¯ ε r , x = Xl , 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 . In terms of the new variables, BVP (1.1) and (1.2) becomes, 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.3)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.4)One often imposes the electroneutrality conditions on the concentrations to avoid sharpboundary layers n (cid:88) s =1 z s l s = n (cid:88) s =1 z s r s = 0 . (1.5)The electrochemical potential ˆ µ k ( x ) = ˆ µ idk ( x ) + ˆ µ exk ( x ) for the k th ion species consists of theideal component ˆ µ idk ( x ) and the excess component ˆ µ exk ( x ), where the ideal component isˆ µ idk ( x ) = z k φ ( x ) + ln c k ( x ) . (1.6)The classical PNP model only deals with the ideal component ˆ µ idk ( x ), which reflects the collisionbetween ion particles and water molecules and ignores the size of ions. The excess electrochemicalpotential ˆ µ exk ( x ) accounts for the finite size effect of ions. This component is essential for dealingwith properties of crowded ionic mixtures where concentrations exceed say 1M.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 (1.3) and (1.4), then the current I is I = I ( ε ) = n (cid:88) s =1 z s J s ( ε ) . (1.7)We will be interested in the zero order approximation of I = I (0) and J k = J k (0). Note that, J k depends on V , Q ( x ), l k ’s and r k ’s, so is I . As mentioned before, we will focus mainly on the4ependance of I = I ( V, Q ) on the electric potential V and permanent charge Q . Particularly,for fixed Q , the electric potential V so that I ( V, Q ) = 0 is the reversal potential . The rever-sal potential has been used to identify the type (i.e., selectivity) of ion channels in biologicalexperiments since 1949 ([18, 19]). Similarly, for fixed V , the permanent charge Q that makes I ( V, Q ) = 0 is called a reversal permanent charge as introduced in [11]. For the existence ofa reversal permanent charge Q of a general form, a necessary condition is that the quantities z k ( z k V + ln l k − ln r k ) , for k = 1 , , . . . , n , cannot have the same sign (Proposition 1.1 in [11]).In [11], the authors presented a necessary condition for the existence of reversal permanentcharges. More precisely, one has J k (cid:90) h ( x ) D k ( x ) c k ( x ) dx = ˆ µ (0) − ˆ µ k (1) . For the classical PNP model where ˆ µ k ( x ) = ˆ µ idk ( x ) = z k φ ( x ) + ln c k ( x ), J k (cid:90) h ( x ) D k ( x ) c k ( x ) dx = z k V + ln l k − ln r k . (1.8)Therefore, if the quantities z k ( z k V +ln l k − ln r k ) , for k = 1 , , . . . , n , have the same sign, then thecurrent I cannot be zero, independent of a permanent charge Q . A question is that under whatconditions on V , l k ’s and r k ’s, can current I be reversed for appropriate choices of permanentcharges Q ? Such a permanent charge Q is called a reversal permanent charge . We now specify the case treated in this paper. We will examine the question by workingon the simplest model, the classical PNP (cPNP) model (1.3) with the ideal electrochemicalpotential ˆ µ k = z k φ + ln c k , a simple profile of a permanent charge Q ( x ) (see (A2) below), andthe boundary condition (1.4). We will focus on the case of two ion species but allow differentdiffusion coefficients. More precisely, we will assume (A0) ˆ ε ( x ) = 1 and, for ≤ k ≤ n , D k ( x ) = D k is a constant;(A1) Electroneutrality boundary conditions (1.5);(A2) A piecewise constant permanent charge Q with one nonzero region; that is, for a partition x < x < x < x = 1 of [0 , , Q ( x ) = (cid:26) Q = Q = 0 , x ∈ ( x , x ) ∪ ( x , x ) ,Q , x ∈ ( x , x ) , (1.9) where Q is an arbitrary constant. For permanent charges Q of the form in (1.9) and for general n , under the condition of equaldiffusion coefficients D k ’s, the topics on the reversal potential and reversal permanent chargeswere examined completely in [11]. It turns out that the condition of equal diffusion coefficientsis highly degenerate (see Remark 2.1), even for n = 2. This is the main technical reason forus to limit to the case n = 2 in this work. As in [11], if the permanent charge Q in (1.9) is areversal permanent charge, we simply call Q a reversal permanent charge.5 GSP for the BVP (1.3) and the results on current reversal forthe case study ( n = 2 and z = − z > ). In [32], a geometric singular perturbation (GSP) framework, combining with special structures ofPNP systems, has been developed for studying the BVP (1.3) and (1.4). This general dynamicalsystem framework and the subsequent analysis have demonstrated the great power of analyzingPNP type problems with potential boundary and internal layers (see [10, 31, 32, 35] for studyon cPNP models, [30] for PNP with a local excess hard-sphere components, and [28, 33, 40] forPNP with nonlocal excess hard-sphere components).For convenience, we will give a brief account of the relevant results in [32] (with slightlydifferent notations) and refer the readers to the paper for details. We remind the readers thatwe will work on cPNP with ideal electrochemical potential µ k = z k φ + ln c k .Denote the derivative with respect to x by overdot and introduce u = ε ˙ φ and w = x .System (1.3) becomes, for k = 1 , , . . . , n , ε ˙ φ = u, ε ˙ u = − n (cid:88) s =1 z s c s − Q ( w ) − ε h w ( w ) h ( w ) u,ε ˙ c k = − z k c k u − ε J k D k h ( w ) , ˙ J k = 0 , ˙ w = 1 . (2.1)System (2.1) will be treated as a dynamical system with the phase space R n +3 and theindependent variable x is viewed as time for the dynamical system. The boundary condition (1.4)becomes, for k = 1 , , . . . , n , φ (0) = V, c k (0) = l k , w (0) = 0; φ (1) = 0 , c k (1) = r k , w (1) = 1 . Following the framework in [32], we convert the boundary value problem to a connectingproblem. To this end, we denote C = ( c , c , . . . , c n ) T and J = ( J , J , . . . , J n ) T , and for j = 1 ,
2, preassign values of φ and C at x j : φ ( x j ) = φ [ j ] and C ( x j ) = C [ j ] . Now for j = 0 , , ,
3, let B j be the subsets of the phase space R n +3 defined by B j = (cid:110) ( φ, u, C, J, w ) : φ = φ [ j ] , C = C [ j ] , w = x j (cid:111) , (2.2)Note that dim B j = n + 1. Then, the BVP (1.3) and (1.4) is equivalent to the following connecting orbit problem : finding an orbit of (2.1) from B to B . The construction would beaccomplished 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 connectingorbit for ε > B j − to B j for j = 1 , , J k and u at x and x (see again [32] for details). B j − to B j . By setting ε = 0 in system (2.1), we get the slow manifold Z j = (cid:110) u = 0 , n (cid:88) s =1 z s c s + Q j = 0 (cid:111) .
6n terms of the independent variable ξ = x/ε , we obtain the fast system of (2.1), for k =1 , , . . . , n , φ (cid:48) = u, u (cid:48) = − n (cid:88) s =1 z s c s − 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.3)where prime denotes the derivative with respect to ξ . The limiting fast system is, for k =1 , , . . . , n , φ (cid:48) = u, u (cid:48) = − n (cid:88) s =1 z s c s − Q j , c (cid:48) k = − z k c k u, J (cid:48) = 0 , w (cid:48) = 0 . (2.4)The slow manifold Z j is precisely the set of equilibria of (2.4) with dim Z j = 2 n + 1. Forthe linearization of (2.4) at each point on Z j , there are (2 n + 1) zero eigenvalues associated tothe tangent space of Z j and the other two eigenvalues are ± (cid:112)(cid:80) ns =1 z s c s . Thus, Z j is normallyhyperbolic (see [12, 17]). We will denote the stable and unstable manifolds of Z j by W s ( Z j )and W u ( Z j ), respectively.Let M [ j − , +] be the collection of all forward orbits from B j − under the flow of (2.4) andlet M [ j, − ] be the collection of all backward orbits from B j . Then the set of forward orbitsfrom B j − to Z j is N [ j − , +] = M [ j − , +] ∩ W s ( Z j ), and the set of backward orbits from B j to Z j is N [ j, − ] = M [ j, − ] ∩ W u ( Z j ). Therefore, the singular layer Γ [ j − , +] at x j − satisfiesΓ [ j − , +] ⊂ N [ j − , +] and the singular layer Γ [ j, − ] at x j satisfies Γ [ j, − ] ⊂ N [ j, − ] .All those important geometric objects are explicitly characterized in [32]. x and x . The limiting fast (layer) dynamics conserve electrochemical potentials, and hence, do not dependon diffusion constants (see, e.g. Proposition 3.3 in [32]). We thus can apply the result aboutthe fast dynamics from [11] directly and only point out the differences. The relevant results areLemmas 3.4 and 3.5 in [11]. The differences are that we have to keep φ [1 , − ] , φ [1 , +] , φ [2 , − ] and φ [2 , +] here in this paper, while in [11] it is known that φ [1 , − ] = V , φ [1 , +] = φ [2 , − ] (denoted by V ∗ there) and φ [2 , +] = 0. With this modification, these lemmas are cast below for n = 2. Lemma 2.1.
The fast layer dynamics over x = x provides, for k = 1 , ,(i) relative to (0 , x ) where Q = 0 , z c [1]1 e z ( φ [1] − φ [1 , − ] ) + z c [1]2 e z ( φ [1] − φ [1 , − ] ) = 0 , c [1 , − ] k = c [1] k e z k ( φ [1] − φ [1 , − ] ) ; (ii) relative to ( x , x ) where Q (cid:54) = 0 , z c [1]1 e z ( φ [1] − φ [1 , +] ) + z c [1]2 e z ( φ [1] − φ [1 , +] ) + Q = 0 , c [1 , +] k = c [1] k e z k ( φ [1] − φ [1 , +] ) ; (iii) the matching u [1] − = u [1]+ : c [1 , − ]1 + c [1 , − ]2 = c [1 , +]1 + c [1 , +]2 + Q ( φ [1] − φ [1 , +] ) . Lemma 2.2.
The fast layer dynamics over x = x provides, for k = 1 , , i) relative to ( x , x ) where Q (cid:54) = 0 , z c [2]1 e z ( φ [2] − φ [2 , − ] ) + z c [2]2 e z ( φ [2] − φ [2 , − ] ) + Q = 0 , c [2 , − ] k = c [2] k e z k ( φ [2] − φ [2 , − ] ) ; (ii) relative to ( x , where Q = 0 , z c [2]1 e z ( φ [2] − φ [2 , +] ) + z c [2]2 e z ( φ [2] − φ [2 , +] ) = 0 , c [2 , +] k = c [2] k e z k ( φ [2] − φ [2 , +] ) ; (iii) the matching u [2] − = u [2]+ : c [2 , − ]1 + c [2 , − ]2 + Q ( φ [2] − φ [2 , − ] ) = c [2 , +]1 + c [2 , +]2 . ( x j − , x j ) . The degeneracy of equal diffusion coefficients shows in the slow dynamics. We will point outthe exact place in the following construction of the slow orbits over the slow manifold Z j = (cid:110) u = 0 , n (cid:88) s =1 z s c s + Q j = 0 (cid:111) . Note that system (2.1) is degenerate at ε = 0 in the sense that all dynamical information on( φ, c , · · · , c n ) would be lost when setting ε = 0. In [32], the dependent variables are rescaled as u = εp, z n c n = − n − (cid:88) s =1 z s c s − Q j − εq. Replacing ( u, c n ) with ( p, q ), slow system (2.1) becomes, for k = 1 , , . . . , n − φ = p, ε ˙ p = q, ˙ c k = − z k pc k − J k D k h ( w ) , ˙ J = 0 , ˙ w = 1 ,ε ˙ q = (cid:16) n − (cid:88) s =1 ( z s − z n ) z s c s − z n Q j − εz n q (cid:17) p + 1 h ( w ) n (cid:88) s =1 z s J s D s . The limiting slow system is, for k = 1 , , . . . , n − φ = p, q = (cid:16) n − (cid:88) s =1 ( z s − z n ) z s c s − z n Q j (cid:17) p + 1 h ( w ) n (cid:88) s =1 z s J s D s = 0 , ˙ c k = − z k pc k − J k D k h ( w ) , ˙ J = 0 , ˙ w = 1 . (2.5)Therefore, on the new slow manifold S j = (cid:40) p = − h − ( w ) (cid:80) ns =1 z s J s D s (cid:80) n − s =1 ( z s − z n ) z s c s − z n Q j , q = 0 (cid:41) , system (2.5) reads, for k = 1 , , . . . , n − φ = − h − ( w ) (cid:80) ns =1 z s J s D s (cid:80) n − s =1 ( z s − z n ) z s c s − z n Q j , ˙ c k = h − ( w ) (cid:80) ns =1 z s J s D s (cid:80) n − s =1 ( z s − z n ) z s c s − z n Q j z k c k − h − ( w ) J k D k , ˙ J =0 , ˙ w = 1 . (2.6)8n S j where q = (cid:80) ns =1 z s c s + Q j = 0, it follows that n − (cid:88) s =1 ( z s − z n ) z s c s − z n Q j = n (cid:88) s =1 z s c s . Remark . Note that, with equal diffusion constant condition, the zero current I = (cid:80) s z s J s = 0reduces system (2.6) to ˙ φ =0 , ˙ c k = − J k h ( w ) , ˙ J = 0 , ˙ w = 1 . The system can be solved explicitly and the solution is simple enough which is the very reasonfor the authors in [11] to obtain their rather specific results for general n . This is NOT the caseif D k ’s are not the same. In order to get reasonably explicit solution that can lead to advancesof understanding of the physical problem, one has serious trouble to treat even the case with n = 2. In fact, we can only handle the case where n = 2 and z = − z at this moment.We now get back to system (2.6) and apply the assumption that n = 2, z = − z and I = z J + z J = 0 (so J = J ). In this case, z J D + z J D = D − D D D z J . (2.7)Applying zero current condition (2.7), the limiting slow system (2.6) becomes,˙ φ = − ( D − D ) J D D h ( τ )(2 z c + Q ) , ˙ c = − ( D + D ) z c + D QD D h ( τ )(2 z c + Q ) J , ˙ J = 0 , ˙ w = 1 . (2.8) Slow system (2.8) on ( x , x ) with Q = 0 : ˙ φ = − ( D − D ) J D D h ( τ ) z c , ˙ c = − D + D D D h ( τ ) J , ˙ J = 0 , ˙ w = 1 . (2.9)The solution of (2.9) with the initial condition ( V, l , J ,
0) is, w ( x ) = x, c ( x ) = l − D + D D D J H ( x ) , φ ( x ) = V − D − D z ( D + D ) ln c ( x ) l , where H ( x ) = (cid:82) x h ( s ) ds . Evaluating the solution at w = x = x we get the following lemma. Lemma 2.3.
Over (0 , x ) with z c ( x ) + z c ( x ) = − Q = 0 the slow dynamics system gives, c [1 , − ]1 = l − D + D D D J H ( x ) , φ [1 , − ] = V − D − D z ( D + D ) ln c [1 , − ]1 l . Slow system (2.8) on ( x , x ) with Q = Q (cid:54) = 0 : Note that h ( w ) >
0. Also, c k ’s are theconcentrations of ion species. Therefore, we will be interested in solutions with c k > k = 1 ,
2, and hence ( z − z ) z c − z Q = z c + z c > h ( w )(( z − z ) z c − z Q ) > ddy φ = D − D D D z J , ddy c = − ( D + D ) z J D D c + z QD J ,ddy J =0 , ddy w = h ( w ) (cid:0) z c + z Q (cid:1) . (2.10)9he solution of (2.10) with the initial condition ( φ [1 , +] , c [1 , +]1 , J , x ) is, φ ( y ) = φ [1 , +] + D − D D D z J y,c ( y ) = e − D D D D z J y c [1 , +]1 + D Q ( D + D ) z (cid:16) e − D D D D z J y − (cid:17) , (cid:90) wx h ( s ) ds = − D D z c [1 , +]1 ( D + D ) J (cid:16) e − D D D D z J y − (cid:17) − D z QD + D (cid:16) D D e − D D D D z J y − D + D ) z J + y (cid:17) + z yQ. (2.11)Assume w ( y ∗ ) = x for some y ∗ >
0, then φ ( y ∗ ) = φ [2 , − ] and c ( y ∗ ) = c [2 , − ]1 . Then, from (2.11)one has the following result. Lemma 2.4.
Over ( x , x ) with z c ( x ) + z c ( x ) + Q = 0 the slow dynamics system gives, φ [2 , − ] = φ [1 , +] + D − D D D z J y ∗ ,c [2 , − ]1 = e − D D D D z J y ∗ c [1 , +]1 + D Q ( D + D ) z (cid:0) e − D D D D z J y ∗ − (cid:1) ,J = − D D (cid:0) c [2 , − ]1 − c [1 , +]1 (cid:1) − ( φ [2 , − ] − φ [1 , +] ) Q ( D + D ) (cid:0) H ( x ) − H ( x ) (cid:1) . Slow system (2.8) on ( x , x ) with Q = 0 : The slow dynamics system is (2.9) and the so-lution with the initial condition ( φ [2 , +] , c [2 , +]1 , J , x ) is, c ( x ) = c [2 , +]1 − D + D D D J (cid:0) H ( x ) − H ( x ) (cid:1) , φ ( x ) = φ [2 , +] − D − D z ( D + D ) ln c ( x ) c [2 , +]1 . Evaluating the solution at w = x = 1 we get the following lemma. Lemma 2.5.
Over ( x , with z c ( x ) + z c ( x ) = 0 the slow dynamics system gives, c [2 , +]1 = r + D + D D D J (cid:0) H (1) − H ( x ) (cid:1) , φ [2 , +] = D − D z ( D + D ) ln r c [2 , +]1 . .2 Matching for zero current and singular orbits on [0 , . The matching conditions are u [1] − = u [1]+ , u [2] − = u [2]+ , and J has to be the same on all subintervals.It follows from Lemmas 2.1 and 2.2, and equations in Lemma 2.3 to Lemma 2.5, z c [1]1 e z ( φ [1] − φ [1 , − ] ) + z c [1]2 e z ( φ [1] − φ [1 , − ] ) = 0 ,z c [2]1 e z ( φ [2] − φ [2 , +] ) + z c [2]2 e z ( φ [2] − φ [2 , +] ) = 0 ,z c [1]1 e z ( φ [1] − φ [1 , +] ) + z c [1]2 e z ( φ [1] − φ [1 , +] ) + Q = 0 ,z c [2]1 e z ( φ [2] − φ [2 , − ] ) + z c [2]2 e z ( φ [2] − φ [2 , − ] ) + Q = 0 , − z − z z c [1 , − ]1 = c [1]1 e z ( φ [1] − φ [1 , +] ) + c [1]2 e z ( φ [1] − φ [1 , +] ) + Q ( φ [1] − φ [1 , +] ) , − z − z z c [2 , +]1 = c [2]1 e z ( φ [2] − φ [2 , − ] ) + c [2]2 e z ( φ [2] − φ [2 , − ] ) + Q ( φ [2] − φ [2 , − ] ) ,J D D = − c [1 , − ]1 − l )( D + D ) H ( x ) = − r − c [2 , +]1 )( D + D )( H (1) − H ( x )) , = − z ( c [2 , − ]1 − c [1 , +]1 ) − ( φ [2 , − ] − φ [1 , +] ) z Q ( D + D )( H ( x ) − H ( x )) ,φ [2 , − ] = φ [1 , +] + D − D D D z J y ∗ ,c [2 , − ]1 = e − D D D D z J y ∗ c [1 , +]1 + D Q ( D + D ) z (cid:16) e − D D D D z J y ∗ − (cid:17) , (2.12)where, c [1 , − ]1 = c [1]1 e z ( φ [1] − φ [1 , − ] ) = (cid:113) c [1]1 c [1]2 , c [1 , − ]2 = c [1]2 e z ( φ [1] − φ [1 , − ] ) = (cid:113) c [1]1 c [1]2 ,c [2 , +]1 = c [2]1 e z ( φ [2] − φ [2 , +] ) = (cid:113) c [2]1 c [2]2 , c [2 , +]2 = c [2]2 e z ( φ [2] − φ [2 , +] ) = (cid:113) c [2]1 c [2]2 ,c [1 , +]1 = c [1]1 e z ( φ [1] − φ [1 , +] ) , c [1 , +]2 = c [1]2 e z ( φ [1] − φ [1 , +] ) ,c [2 , − ]1 = c [2]1 e z ( φ [2] − φ [2 , − ] ) , c [2 , − ]2 = c [2]2 e z ( φ [2] − φ [2 , − ] ) ,φ [1 , − ] = V − D − D z ( D + D ) ln c [1 , − ]1 l , φ [2 , +] = D − D z ( D + D ) ln r c [2 , +]1 . (2.13) Remark . In (2.12), the unknowns are: φ [1] , φ [2] , c [1]1 , c [1]2 , c [2]1 , c [2]2 , J , φ [1 , +] , φ [2 , − ] , y ∗ and Q that is, there are eleven unknowns that matches the total number of equations on (2.12).It follows from last two equations of (2.13), φ [1] − φ [1 , − ] = D + D D (cid:0) φ [1] − V (cid:1) + D − D z D ln c [1]1 l ,φ [2] − φ [2 , +] = D + D D φ [2] − D − D z D ln r c [2]1 . (2.14) Remark . The matching conditions in (2.12) are (kind of) equivalent to the matching con-dition in [10], in a sense that, here we added zero current conditions and also have diffusion11oefficients in our equations in (2.12). We went through the same procedure here as the authorsof [10] did to find F ( A ) = 0 in their paper though. Recall that to find F ( A ) = 0, the authors of[10] first obtained (See the last equation of (44) in [10]), (cid:113) Q + B − Q = e − ( J + J ) y ∗ (cid:0)(cid:113) Q + A − Q (cid:1) − Q J J + J (cid:0) − e − ( J + J ) y ∗ (cid:1) . Now, applying zero current z J + z J = 0 with z = − z = 1 one has J = J and the aboveequation becomes, e − J y ∗ (cid:113) Q + A − (cid:113) Q + B = 0 . One can see that the same equation will be obtained here from the last equation of (7.4), withthe constraints, D = D = 1 and z = 1. z = − z > . The matching system (2.12) is nonlinear and challenging to analyze in general. In [11], for equaldiffusion constants D k ’s, the study of reversal potential and reversal permanent charges hasbeen successfully carried out for a general n . It is a little bit surprising that with general D k ’sthe problem becomes overwhelmingly harder, at least, technically, even for the case that we willtreat here where n = 2 with z = − z .In [10], the authors introduced two intermediate variables that allow a significant reductionof the governing system of matching (2.12) without zero current assumption. We will use thesame intermediate variables for our reduction. Thus, we set A = (cid:113) c [1]1 c [1]2 and B = (cid:113) c [2]1 c [2]2 . (3.1)It will be shown in (7.11) that B = B ( A ) = 1 − βα ( l − A ) + r . We will thus treat B as a functionof A instead of an independent variable from now on. We denote l = l = l, r = r = r, Q = 2 Q , α = H ( x ) H (1) , β = H ( x ) H (1) , θ = D − D D + D . (3.2)Note that 0 < α < β < − < θ <
1, which will always be assumed. The vector( Q , V, θ, α, β, l, r ) contains major parameters of the system, which affect the behavior of thesystem through their nonlinear interactions. In the sequel, we will always fix the parameters α , β , l and r , and focus on the roles of ( V, Q , θ ). One can see that the roles of ( α, β, l, r ) can bestudied within our analysis framework. For ease of notation, we also introduce S a := (cid:113) Q + z A , S b := (cid:113) Q + z B , N := β − αα z ( A − l ) + S a − S b . (3.3)The most critical ingredient for our analysis is the following result on a reduced system ofthe matching system (2.12). Proposition 3.1.
The matching system (2.12) for zero current I = 0 can be reduced to G ( A, Q , θ ) = z V and G ( A, Q , θ ) = 0 , (3.4)12 here G ( A, Q , θ ) = θ (cid:16) ln S a + θQ S b + θQ + ln lr (cid:17) − (1 + θ ) ln AB + ln S a − Q S b − Q ,G ( A, Q , θ ) = θQ ln S a + θQ S b + θQ − N. (3.5) Proof.
We defer the proof to the appendix Section 7.At this moment, we would like to make some comments on the above reduction.
Remark . The reduction of (2.12) to system (3.4) is critical for the remaining analysis. Wecomment that there is no practical principle to lead the reduction and no criterion for a ‘good’final form of a reduction. In general, there could be infinitely many different forms of thereduction. It turns out the above reduction works well.For the special case where h = 1, x = 1 / x = 2 / z = 1 = − z , and D = D , a reducedsystem consists of F ( A ) = 0 in (48) in [10] and I = 0. One can get different equivalent formsand, as expected, one equivalent reduced system can be put into exactly the same as the onestated in Proposition 3.1. We also note that, for a given Q , one cannot solve for A from either F ( A ) = 0 or I = 0 uniquely. But, we will show that one can solve for A from G = 0 uniquely– a critically important indication that the specific form of system (3.4) is special.We now prepare several properties of the functions G and G to be used later on. Lemma 3.2.
One has(i) ∂ A G ( A, Q , θ ) has the same sign as that of Q ,(ii) ∂ Q G ( A, Q , θ ) has the same sign as that of l − r ,(iii) ∂ θ G ( A, Q , θ ) has the same sign as that of l − r ,(iv) ∂ A G ( A, Q , θ ) < ,(v) if θQ > , then ∂ Q G ( A, Q , θ ) has the same sign as that of ( l − r ) Q ,(vi) ∂ θ G ( A, Q , θ ) has the same sign as that of ( l − r ) Q .Proof. Partial 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 , θ ) = g ( S a ) − g ( S b ) + ln lr − ln AB ,∂ A G ( A, Q , θ ) = − − βα z BS b + θQ − z AS a + θQ − β − αα z ,∂ 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 ( g ( S a ) − g ( S b )) , (3.6)13here g ( X ) := ln( X + θQ ) + θQ X + θQ . (3.7)All statements except those for signs of ∂ θ G k ’s follow directly from (3.6). For signs of ∂ θ G k ’s,note that g (cid:48) ( X ) = X ( X + θQ ) > X > . So g ( S a ) − g ( S b ) has the same sign as that of S a − S b . It is obvious that S a − S b has the samesign as that of A − B and it will be shown in Theorem 3.4 that l − r and A − B have the samesign too. The statements on the signs of ∂ θ G k ’s then follow. A = A ( Q , θ ) of G ( A, Q , θ ) = 0 . Recall that B = 1 − βα ( l − A ) + r . One has B = A if and only if A ∗ = (1 − β ) l + αr − β + α . It is clearthat l < A ∗ < r if l < r and l > A ∗ > r if l > r . Theorem 3.3.
For any given ( Q , θ ) , G ( A, Q , θ ) = 0 has a unique solution A = A ( Q , θ ) .Proof. For any ( Q , θ ), it follows from Lemma 3.2 that ∂ A G ( A, Q , θ ) <
0, and hence, G ( A, Q , θ )is strictly decreasing in A . Let A M = l + αr/ (1 − β ) be the largest value for A (when B = 0)and let B M = (1 − β ) l/α + r be the largest value for B (when A = 0).Set x = (cid:113) Q + z B M > | Q | and y = (cid:113) Q + z A M > | Q | . Then, G (0 + , Q , θ ) = f ( x ) := θQ ln | Q | + θQ x + θQ + β − αα z l − | Q | + x,G ( A − M , Q , θ ) = f ( y ) := θQ ln y + θQ | Q | + θQ − β − αα z ( A M − l ) − y + | Q | . It is easy to check that f (cid:48) ( t ) > > f (cid:48) ( t ) for t >
0, and hence, f ( x ) > f ( | Q | ) = β − αα z l > f ( y ) < f ( | Q | ) = − β − αα z ( A M − l ) < . Thus, for any ( Q , θ ) there is a unique A = A ( Q , θ ) such that G ( A ( Q , θ ) , Q , θ ) = 0.In the following, we also denote B ( A ( Q , θ )) by B ( Q , θ ). Theorem 3.4.
The solution A = A ( Q , θ ) of G ( A, Q , θ ) = 0 satisfies(a) A (0 , θ ) = (1 − α ) l + αr and lim Q →±∞ A ( Q , θ ) = l ,(b) if l > r , then l > A ( Q , θ ) > A ∗ > B ( Q , θ ) > r ,(c) if l < r , then l < A ( Q , θ ) < A ∗ < B ( Q , θ ) < r ,(d) if θQ ≥ , then ∂ Q A ( Q , θ ) has the same sign as that of ( l − r ) Q . roof. (a). The value A (0 , θ ) can be deduced from G ( A, , θ ) = − β − αα z ( A − l ) − z ( A − B ) = 0 and B = 1 − βα ( l − A ) + r. For the claim about the limits, one has, from G ( A ( Q , θ ) , Q , θ ) = 0,lim Q →±∞ θQ ln S a + θQ S b + θQ = lim Q →±∞ (cid:16) β − αα z ( A − l ) + S a − S b (cid:17) = β − αα z lim Q →±∞ ( A − l ) . On the other hand, apply L’Hospital rule to getlim Q →±∞ Q ln S a + θQ S b + θQ = − lim Q →±∞ ( Q S a + θ )( S b + θQ ) − ( Q S b + θ )( S a + θQ )( S a + θQ )( S b + θQ ) Q = 0 . Thus, lim Q →±∞ A ( Q , θ ) = l .(b). Recall that, for A = A ∗ = (1 − β ) l + αr − β + α , B = A ∗ . Thus, G ( A ∗ , Q , θ ) = β − αα z ( l − A ∗ ) = β − α − β + α z ( l − r ) . Note that, for some S ∗ between S a and S b ,ln S a + θQ S b + θQ = ln( S a + θQ ) − ln( S b + θQ ) = S a − S b S ∗ + θQ . Thus, for some S ∗ between S a and S b , G ( l, Q , θ ) = θQ ln S a + θQ S b + θQ − ( S a − S b ) = − ( S a − S b ) S ∗ S ∗ + θQ . If l > r , then G ( A ∗ , Q , θ ) >
0, which yields A ∗ < A ( Q , θ ) since G is decreasing in A . The latter implies B ( Q , θ ) < A ∗ < A ( Q , θ ), and hence, S a > S b , which then implies G ( l, Q , θ ) <
0. Due to again that G is deceasing in A , r < A ∗ < A ( Q , θ ) < l if l > r (independent of Q ).(c). Similarly, if l < r , then G ( A ∗ , Q , θ ) < < G ( l, Q , θ ), and hence, A ∗ > A ( Q , θ ) > l .(d). It follows from (3.6) that, if θQ > θ = 0, then ∂ Q G has the same sign as thatof ( S a − S b ) Q . The latter has the same sign as that of ( l − r ) Q . The statement then followsfrom ∂ A G < ∂ Q A = − ∂ Q G /∂ A G . Remark . Note that, with zero current condition I = 0, we have that A ( Q , θ ) always liesbetween l and r for any Q . This is not true without zero current condition (see [41]).We believe that, if l (cid:54) = r , then A ( Q , θ ), or equivalently, G ( A ( Q , θ ) , Q , θ ) has a uniquecritical point in Q . It is true if D = D (so θ = 0) but we could not establish it in general.Figure 1 shows numerical simulations of A ( Q , θ ) for α = , β = , l = 2 and r = 3 withdifferent choices of ( D , D ). Both graphs support our belief that A ( Q , θ ) has a unique criticalpoint in Q . Theorem 3.5.
The quantity ∂ θ A ( Q , θ ) has the same sign as that of ( l − r ) Q .Proof. It follows from G = 0 in (3.5) that, ∂ θ A = − ∂ θ G /∂ A G . The statement then followsfrom (iv) and (vi) in Lemma 3.2. 15igure 1: Graph of A ( Q , θ ) with different choices of ( D , D ) : left panel for D > D and rightpanel for D < D . For the case of zero current with z = − z , one has J = J . Denote the equal fluxes by J thatwe call it zero current flux . Once a solution ( A, V ) of G = z V and G = 0 is obtained, itfollows from (7.10) that J is given by J = − D D ( A − l )( D + D ) αH (1) = − D D ( r − B )( D + D )(1 − β ) H (1) . (3.8)Note that the function A = A ( Q , θ ) depends on θ = ( D − D ) / ( D + D ), or equivalently,on ρ = D /D , but J = J ( Q , D , D ) is not homogeneous of degree zero in ( D , D ).The following result is a direct consequence of Theorems 3.4 and 3.5. Corollary 3.6.
The zero current flux J = J ( Q , D , D ) satisfies(a) if θQ ≥ , then ∂ Q J and ( l − r ) Q have opposite signs,(b) if Q > , then ∂ D J and l − r have the same sign,(c) if Q < , then ∂ D J and l − r have the same sign.Proof. Direct calculations from (3.8) give ∂ Q J = − D D ( D + D ) αH (1) ∂ Q A, ∂ D J = (1 + θ ) αH (1) (cid:16) l − A ( Q , θ ) + (1 − θ ) ∂ θ A (cid:17) ,∂ D J = (1 − θ ) αH (1) (cid:16) l − A ( Q , θ ) − (1 + θ ) ∂ θ A (cid:17) . The statement follows from the above formulas and Theorems 3.4 and 3.5.
Remark . We do not know the signs of ∂ Q J , ∂ D J , and ∂ D J under conditions other thanthose in the statement of Corollary 3.6 in general.16 Reversal potential V rev = V rev ( Q , θ ) . We are searching for the value V = V rev of the transmembrane potential V = φ (0) − φ (1) thatproduces zero current I . For the case we considered, we will show the existence and uniquenessof reversal potentials. Theorem 4.1.
Consider ionic flow of two ( n = 2 ) ion species with z = − z . For any given ( Q , θ ) , there exists a unique reversal potential V rev = V rev ( Q , θ ) .Proof. It follows from Theorem 3.3 that for any given ( Q , θ ), there exists a unique A = A ( Q , θ )such that G ( A ( Q , θ ) , Q , θ ) = 0. A reversal potential V = V rev is then determined from z V rev = G ( A ( Q , θ ) , Q , θ ) and is given by V rev = θz (cid:0) ln S a + θQ S b + θQ + ln lr (cid:1) − θz ln A ( Q , θ ) B ( Q , θ ) + 1 z ln S a − Q S b − Q , (4.1)where S a and S b are given in terms of Q and A ( Q , θ ) as in (3.3). V rev on Q . We will consider how the reversal potential V rev = V rev ( Q , θ ) depends on Q . Recall that wedenote J = J by J . Theorem 4.2.
For the reversal potential V rev = V rev ( Q , θ ) , one has(i) if l > r , then J > , and hence, − z ln lr < V rev ( Q , θ ) < z ln lr ;(ii) if l < r , then J < , and hence, z ln lr < V rev ( Q , θ ) < − z ln lr ;(iii) V rev (0 , θ ) = θz ln lr and lim Q →±∞ V rev ( Q , θ ) = ± z ln lr .Proof. (i) It follows from part (b) in Theorem 3.4 and the formula for J in (3.8) that, if l > r ,then J >
0. The range for V rev is a consequence of that fact that J k has the same sign as thatof z k V + ln l/r . (ii) can be established similarly.(iii) The value of V rev (0 , θ ) follows from (4.1) directly. To show the limits, we recall fromTheorem 3.3 that lim Q →±∞ A ( Q , θ ) = l (and hence, lim Q →±∞ B ( Q , θ ) = r ). Note also that − < θ <
1. Therefore,lim Q → + ∞ ln A ( Q , θ ) B ( Q , θ ) = ln lr , lim Q →±∞ ln S a + θQ S b + θQ = 0 , and lim Q → + ∞ ln S a − Q S b − Q = 2 ln lr , lim Q →−∞ ln S a − Q S b − Q = 0 . Using (4.1), one then haslim Q → + ∞ z V rev ( Q , θ ) = θ ln lr − (1 + θ ) ln lr + 2 ln lr = ln lr , lim Q →−∞ z V rev ( Q , θ ) = θ ln lr − (1 + θ ) ln lr = − ln lr . The proof is completed. 17he next result is a direct consequence of Theorem 4.2, whose proof will be omitted.
Corollary 4.3.
One has,(i) If ( D − D )( l − r ) > , then, for some Q < , V rev ( Q , θ ) = 0 ;(ii) if ( D − D )( l − r ) < , then, for some Q > , V rev ( Q , θ ) = 0 . We now provide remarks on the physical basis for results in Theorem 4.2 and Corollary 4.3.
Remark . The statements (i) and (ii) in Theorem 4.2 can be obtained in a direct way asfollows. Note that, in general, J k has the same sign as that of z k V + ln l/r . Thus, if l > r ,then for V ≤ V = − z ln lr , one has J ≤ z V + ln l/r ≤ z V + ln l/r = 0, and J > z V + ln l/r = − z V + ln l/r ≥ − z V + ln l/r = 2 ln l/r >
0. Therefore, if l > r , then V rev > − z ln lr and, similarly, V rev ( Q , θ ) < z ln lr .In [41], it shows that, as Q → + ∞ , J ( Q ) →
0. Thus, J ( Q ) → I ( Q ) = 0. Since J is proportional to − z V rev ( Q , θ ) + ln l/r with a positive proportional constant in general (see[41]), it follows that − z V rev ( Q , θ ) + ln l/r → Q → + ∞ , which is exactly what claimed inTheorem 4.2 for this limit. The other claim follows from the same argument.Statement (i) in Theorem 4.2 says that, if D < D and l > r , then V rev (0 , θ ) >
0. Thismakes sense since, for V = 0 and l > r , J > J >
0, and, with D < D , J < J .To help J more than J to get J = J , one needs to increase V and this is why, in this case, V rev (0 , θ ) >
0. The latter often implies that, if V = 0, then I ( V = 0) <
0, or equivalently, J > J . Thus, intuitively, in order for the zero potential to be a reversal potential, a permanentcharge helping J more than J is needed; that is, the permanent charge should be negative,which agrees with statement (i) in Corollary 4.3. Other statements in Corollary 4.3 can beexplained similarly.Concerning the monotonicity of V rev = V rev ( Q , θ ), we have Theorem 4.4.
For any given θ ∈ ( − , , one has(i) if θ = 0 , then V rev ( Q , θ ) is increasing in Q for l > r and decreasing in Q for l < r ;(ii) if θ > , then, for Q ≥ , V rev ( Q , θ ) is increasing in Q for l > r and decreasing in Q for l < r ;(iii) if θ < , then, for Q ≤ , V rev ( Q , θ ) is increasing in Q for l > r and decreasing in Q for l < r .Proof. It follows from z V rev = G ( A ( Q , θ ) , Q , θ ) and ∂ Q A = − ∂ Q G /∂ A G that ∂ Q V rev = 1 z ∂ A G (cid:16) ∂ Q G ∂ A G − ∂ A G ∂ Q G (cid:17) . (4.2)The statements then follow from Lemma 3.2.We conjecture that V rev ( Q , θ ) is always monotonic in Q but could not prove it. Numericalsimulations in Figure 2 support our conjecture. Figure 2 has two graphs of V = V rev ( Q , θ ) for − ≤ Q ≤
10 with α = 1 / β = 2 / l = 2 < r = 3, and θ = − / θ = 1 / | Q | small, we have 18igure 2: The function V = V rev ( Q , θ ) decreases when l < r as conjectured: left panel for θ < and right panel for θ > . Theorem 4.5.
Near Q = 0 , the reversal potential V rev ( Q , θ ) is approximated by V rev ( Q , θ ) = θz ln lr + 1 − θ z ( β − α )( l − r ) (cid:0) (1 − α ) l + αr (cid:1)(cid:0) (1 − β ) l + βr (cid:1) Q + O ( Q ) . (4.3) Proof.
It follows from (3.6) that ∂ A G ( A, , θ ) = 0 and ∂ Q G ( A (0 , θ ) , , θ ) = 1 − θ z A (0 , θ ) − B (0 , θ ) A (0 , θ ) B (0 , θ ) . Recall that A (0 , θ ) = (1 − α ) l + αr and B (0 , θ ) = (1 − β ) l + βr . One then has ∂ Q V rev (0 , θ ) = 1 z ∂ A G ( A (0 , θ ) , , θ ) ∂ Q A (0 , θ ) + 1 z ∂ Q G ( A (0 , θ ) , , θ )= 1 − θ z ( β − α )( l − r ) (cid:0) (1 − α ) l + αr (cid:1)(cid:0) (1 − β ) l + βr (cid:1) . The expansion (4.3) then follows from that V rev (0 , θ ) = θz ln lr . V rev on ρ . Recall that ρ := D /D and θ = ( ρ − / ( ρ + 1). It is more convenient to view V rev as a functionof ρ instead of θ when we consider its dependence on D and D . Thus we will abuse notationto denote V rev ( Q , θ ) by V rev ( Q , ρ ). Proposition 4.6.
One has ∂ ρ V rev ( Q , ρ ) has the same sign as that of l − r .Proof. Direct calculations from (4.1) give ∂ ρ V rev ( Q , ρ ) = 2 z (1 + ρ ) (cid:16) g ( S a ) − g ( S b ) + ln lB ( Q , θ ) rA ( Q , θ ) (cid:17) , where g ( X ) is defined in (3.7) and is increasing in X for X >
0. In particular, if l > r then g ( S a ) − g ( S b ) >
0. Moreover, it follows from Theorem 3.4 that r < B ( Q , θ ) < A ( Q , θ ) < l .The proof is thus complete for l > r . The case for l < r is similar.19 emark . Proposition 4.6 shows how diffusion coefficients affect reversal potential and revealsa fascinating attribute that may not be completely intuitive at first glance. Indeed, recall theobservation in [11] that, for k = 1 , J k D k (cid:90) h ( x ) c k ( x ) dx = z k V + ln lr . (4.4)The relation of course holds true for the zero current condition: J = J with V = V rev . Now,if we fix D and increase D (so ρ is increasing), then J increases since all but J D in (4.4)are independent of D ([32]), and hence, to satisfy zero current condition, we should increase J . Intuitively increasing V rev seems to accomplish the latter. But this intuition agrees withProposition 4.6 only for l > r and is the exactly opposite for l < r . That is, for l < r , Proposition4.6 says, as ρ increases, V rev ( Q , ρ ) decreases. This counterintuitive behavior could be explainedby the fact that c ( x ) actually depends on V rev and reducing V rev could increase J .From the above discussion, we feel that it is almost impossible to deduce the whole truthwithout analysis. V rev . We will first recall Goldman-Hodgkin-Katz (GHK) equation for the reversal potential V rev andthen make a comparison with our result.Based on essentially the assumption that the electric potential φ ( x ) is linear in x (or theelectric field is constant), Goldman ([14]), and Hodgkin and Katz ([19]) derived an equation (theGHK equation) for the reversal potential, which extends that of Nernst equation for a single ionspecies. Under the assumption, the I-V (current-voltage) relation is given by I = V n (cid:88) k =1 z k D k r k − l k e z k V − e z k V . For the case where n = 2 and z = − z , the GHK equation for the reversal potential is V GHKrev ( ρ ) = 1 z ln r + ρll + ρr . (4.5)The assumption that the electric potential φ ( x ) is linear in x is thought to probably makesense without channel structure; in particular, Q = 0. This is not correct either. In fact, when Q = 0, from Theorem 4.2 or the expansion of V rev ( Q , θ ) near Q = 0 in Theorem 4.5, thereversal potential is V rev (0 , ρ ) = 1 z ρ − ρ + 1 ln lr , which is different from that in (4.5). In our opinion, what is more important is that our resulton the reversal potential is the first for general Q (cid:54) = 0 with different diffusion coefficient. Thus,for n = 2 with z = − z , the GHK equation for reversal potential should be replaced by V rev ( Q , θ ) = 1 z G ( A ( Q , θ ) , Q , θ )with A ( Q , θ ) being the solution of G ( A, Q , θ ) = 0.It is very important to generalize this result to mixtures with more than two ion species.20 Reversal permanent charge Q rev ( V, θ ) . In view of the duality of reversal potential V and the reversal permanent charge Q ∗ , we nowpresent a general result for reversal permanent charge with a given electric potential V . Wecomment that there are differences between these two problems. On one hand, as probablyexpected, reversal potentials should always exist. On the other hand, there is a simple necessarycondition for the existence of the reversal permanent charge Q rev as discussed above. This isindeed established below for the special case of permanent charges Q in (A3). Theorem 5.1.
For n = 2 with z = − z , there exists a reversal permanent charge Q rev if andonly if (cid:0) z V + ln lr (cid:1)(cid:0) z V + ln lr (cid:1) > . (5.1) Proof.
Since J k , for k = 1 ,
2, has the same sign as that of z k V + ln lr and z = − z , the conditionin (5.1) is necessary for a zero current I , and hence, for the existence of a reversal permanentcharge. To show the condition is also sufficient, we set g ( Q ) := G ( A ( Q ) , Q , θ ) − z V . FromTheorem 3.3 one obtains,lim Q → + ∞ ln S a − Q S b − Q = 2 lim Q → + ∞ ln AB = 2 ln lr , lim Q →−∞ ln S a − Q S b − Q = lim Q →±∞ ln S a + θQ S b + θQ = 0 . Then from above, the equation for G in (3.5), Lemma 3.3 and above one has,lim Q → + ∞ g ( Q ) = − z V + ln lr , lim Q →−∞ g ( Q ) = − z V − ln lr . (5.2)The condition (5.1) implies that the above values have opposite signs. By the IntermediateValue Theorem, there is at least one Q = Q rev ( V, θ ) such that g ( Q ) = 0.This existence result can be viewed as a duality of Theorem 4.1 together with (i) and (ii) inTheorem 4.2. The next result is a duality to (iii) of Theorem 4.2, whose proof will be omitted. Theorem 5.2.
For any ( V, l, r ) that satisfies the condition (5.1) one has, lim z V → θ ln l/r Q rev ( V, θ ) = 0 and lim z V →± ln l/r Q rev ( V, θ ) = ±∞ . Recall we could not show but conjecture that V rev ( Q , θ ) is monotone in Q in Section 4.Should the conjecture be shown, Q rev ( V, θ ) would be monotone in V . In this paper, we work on the cPNP model allowing unequal diffusion constants and for a singleprofile of permanent charges, to answer the specific questions about reversal potentials andreversal permanent charges that are among the central issues of biological functions. Our studyrelies on a modern general geometric singular perturbation theory and on some special structuresof the cPNP models. These allow us to obtain a nonlinear matching system of algebraic equations(2.12) for zero current condition that include both the reversal potential and reversal permanentcharge topics. Using an intermediate variable introduced in [10], the matching system is furtherreduced to an effective system of two algebraic equations with two unknowns. A number ofinteresting properties of biological significance are resulted from analysis of these governingequations, some are not totally intuitive. 21
Appendix: Proof of Proposition 3.1.
We consider a special case where z = − z . Set c [1]1 c [1]2 = A and c [2]1 c [2]2 = B . We will use thenotion l , r , Q , α , β , S a , S b and N introduced in (3.2) and (3.3).With above terms, from (2.13) we get c [1 , − ]1 = c [1 , − ]2 = A, c [2 , +]1 = c [2 , +]2 = B, c [1 , +]1 = S a − Q z , c [2 , − ]1 = S b − Q z . (7.1)From the third and fourth equations in (2.12) one has, φ [1] − φ [1 , +] = 1 z ln S a − Q z c [1]1 , φ [2] − φ [2 , − ] = 1 z ln S b − Q z c [2]1 . (7.2)Then, from the first two equations of (2.12), (2.13) and (2.14) give φ [1] = V + 2 D (cid:0) D + D (cid:1) z ln( z A ) + ( D − D ) (cid:0) D + D (cid:1) z ln( z l ) − z ln( z c [1]1 ) ,φ [2] = 2 D (cid:0) D + D (cid:1) z ln( z B ) + ( D − D ) (cid:0) D + D (cid:1) z ln( z r ) − z ln( z c [2]1 ) . (7.3)The rest of system (2.12) becomes, c [1]2 = c [1]1 e z ( φ [1] − φ [1 , − ] ) , c [2]2 = c [2]1 e z ( φ [2] − φ [2 , +] ) ,z A = S a + Q ln S a − Q z c [1]1 , z B = S b + Q ln S b − Q z c [2]1 ,J D D = − A − l )( D + D ) αH (1) = − r − B )( D + D )(1 − β ) H (1)= − c [2 , − ]1 − c [1 , +]1 ) − ( φ [2 , − ] − φ [1 , +] ) Q ( D + D )( β − α ) H (1) ,D − D D D z J y ∗ = φ [1] − φ [2] + 1 z ln c [1]1 (cid:0) S b − Q (cid:1) c [2]1 (cid:0) S a − Q (cid:1) S b − Q = e − D D D D z J y ∗ (cid:0) S a − Q (cid:1) − D Q D + D (cid:0) − e − D D D D z J y ∗ (cid:1) . (7.4)From third and fourth equations in (7.4), c [1]1 = S a − Q z exp (cid:110) S a − z AQ (cid:111) , c [2]1 = S b − Q z exp (cid:110) S b − z BQ (cid:111) . (7.5)The equations (7.3) and (7.5) give φ [1] = V + 2 D (cid:0) D + D (cid:1) z ln( z A ) + D − D ( D + D ) z ln( z l ) − z ln( S a − Q ) − S a − z Az Q ,φ [2] = 2 D (cid:0) D + D (cid:1) z ln( z B ) + D − D ( D + D ) z ln( z r ) − z ln( S b − Q ) − S b − z Bz Q . (7.6)22hus φ [2] − φ [1] = − V + 2 D (cid:0) D + D (cid:1) z ln BA − D − D ( D + D ) z ln lr − z ln S b − Q S a − Q + S a − S b + z ( B − A ) z Q . (7.7)Now, the equation (7.2) and y ∗ equation in sixth line of (7.4) give, φ [2 , − ] − φ [1 , +] = φ [2] − φ [1] − z ln c [1]1 (cid:0) S b − Q (cid:1) c [2]1 (cid:0) S a − Q (cid:1) . But, from third and fourth equations of (7.4),ln c [1]1 (cid:0) S b − Q (cid:1) c [2]1 (cid:0) S a − Q (cid:1) = 1 Q (cid:0) S a − S b + z ( B − A ) (cid:1) , and hence, φ [2 , − ] − φ [1 , +] = − V + 2 D (cid:0) D + D (cid:1) z ln BA − D − D ( D + D ) z ln lr − z ln S b − Q S a − Q . (7.8)Furthermore, it follows from above that, φ [2 , − ] − φ [1 , +] = φ [2] − φ [1] − z Q (cid:0) S a − S b + z ( B − A ) (cid:1) = D − D D D z J y ∗ . (7.9)Thus, J equations in (7.4), with equations in (7.9) and (7.1) give, J D D = − A − l )( D + D ) αH (1) = − r − B )( D + D )(1 − β ) H (1)= − B − A − Q ( φ [2] − φ [1] )( D + D )( β − α ) H (1) . (7.10)Now, from the equations in (7.10), B = 1 − βα ( l − A ) + r, φ [2] − φ [1] = − A − l + α ( l − r ) αQ . (7.11)Thus, the equations in (7.9) and (7.11) give, J y ∗ D D = 1 z (cid:0) D − D (cid:1) Q N ( A, Q ) , (7.12)where N = N ( A, Q ) = β − αα z ( A − l ) + S a − S b . is defined in (3.3). On the other hand, from(7.7) and (7.11) we obtain an equation in terms of A and Q ,( β − α )( A − l ) αQ = V − D (cid:0) D + D (cid:1) z ln BA + D − D ( D + D ) z ln lr + 1 z ln S b − Q S a − Q − S a − S b z Q . Now, it follows from above equation and the expression for N ( A, Q ) that, NQ − z V − D D + D ln AB − D − D D + D ln lr + ln S a − Q S b − Q = 0 . (7.13)23ubstituting (7.12) into the last equation of (7.4) we get the other equation for A and Q , (cid:113) Q + z B − Q = e ( D D D − D Q N (cid:0)(cid:113) Q + z A − Q (cid:1) − D Q D + D (cid:0) − e ( D D D − D Q N (cid:1) , that is equivalent to ( D − D ) Q ( D + D ) ln S a + ( D − D ) Q ( D + D ) S b + ( D − D ) Q ( D + D ) − N = 0 . This equation is G = 0 in Proposition 3.1. Also, adding G Q to equation (7.13) one obtains G = z V in Proposition 3.1. This completes the proof of Proposition 3.1. Acknowledgement.
Hamid Mofidi thanks Mingji Zhang and both authors thank Bob Eisen-berg for helpful discussions. Weishi Liu’s research is partially supported by Simons FoundationMathematics and Physical Sciences-Collaboration Grants for Mathematicians
References [1] N. Abaid, R. S. Eisenberg, and W. Liu, Asymptotic expansions of I-V relations via a Poisson-Nernst-Planck system.
SIAM J. Appl. Dyn. Syst. (2008), 1507-1526.[2] V. Barcilon, Ion flow through narrow membrane channels: Part I. SIAM J. Appl. Math. (1992),1391-1404.[3] V. Barcilon, D.-P. Chen, and R. S. Eisenberg, Ion flow through narrow membrane channels: PartII. SIAM J. Appl. Math. (1992), 1405-1425.[4] V. Barcilon, D.-P. Chen, R. S. Eisenberg, and J. W. Jerome, Qualitative properties of steady-state Poisson-Nernst-Planck systems: Perturbation and simulation study. SIAM J. Appl. Math. (1997), 631-648.[5] 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.[6] D. P. Chen and R.S. Eisenberg, Charges, currents and potentials in ionic channels of one confor-mation. Biophys. J. (1993), 1405-1421.[7] S. Chung and S. Kuyucak. Predicting channel function from channel structure using Browniandynamics simulations. Clin. Exp. Pharmacol Physiol. (2001), 89-94.[8] B. Eisenberg, Ion Channels as Devices. J. Comp. Electro. (2003), 245-249.[9] 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.[10] B. Eisenberg and W. Liu, Poisson-Nernst-Planck systems for ion channels with permanent charges.
SIAM J. Math. Anal. (2007), 1932-1966.[11] B. Eisenberg, W. Liu, and H. Xu, Reversal permanent charge and reversal potential: case studiesvia classical Poisson-Nernst-Planck models. Nonlinearity (2015), 103-128.[12] N. Fenichel, Geometric singular perturbation theory for ordinary differential equations. J. Differ-ential Equations (1979), 53-98.[13] D. Gillespie, W. Nonner, and R. S. Eisenberg, Coupling Poisson-Nernst-Planck and density func-tional theory to calculate ion flux. J. Phys.: Condens. Matter (2002), 12129-12145.
14] D. E. Goldman, Potential, impedance, and rectification in membranes,
J. Gen. Physiol. (1943),37-60.[15] B. Hille, Ion Channels of Excitable Membranes (Third Edition).
Sinauer Associates, Inc., Sunder-land, Massachusetts, USA 2001.[16] B. Hille, Transport Across Cell Membranes: Carrier Mechanisms, Chapter 2. Textbook of Phys-iology (ed. H. D. Patton, A. F. Fuchs, B. Hille, A. M. Scher and R. D. Steiner). Philadelphia,Saunders (1989), 24-47.[17] M. Hirsch, C. Pugh, and M. Shub, Invariant Manifolds.
Lecture Notes in Math. , Springer-Verlag, New York, 1976.[18] A. L. Hodgkin, A. Huxley, and B. Katz, Ionic currents underlying activity in the giant axon of thesquid.
Arch. Sci. Physiol. (1949), 129-150.[19] A. L. Hodgkin and B. Katz, The effect of sodium ions on the electrical activity of the giant axonof the squid. J. Physiol. (1949), 37-77.[20] A. L. Hodgkin, The ionic basis of electrical activity in nerve and muscle.
Biological Rev. (1951),339-409.[21] A. L. Hodgkin and A. F. Huxley, Currents carried by sodium and potassium ions through themembrane of the giant axon of Loligo . J. Physol. (1952), 449-472.[22] A. L. Hodgkin and A. F. Huxley, The components of membrane conductance in the giant axon of
Loligo . J Physiol. (1952), 473-496.[23] A. L. Hodgkin and A. F. Huxley, A quantitative description of membrane current and its applicationto conduction and excitation in nerve.
J. Physiol. (1952), 500-544.[24] A. L. Hodgkin and R. D. Keynes, The potassium permeability of a giant nerve fibre.
J. Physiol. (1955), 61-88.[25] U. Hollerbach, D.-P. Chen, and R.S. Eisenberg, Two- and Three-Dimensional Poisson-Nernst-Planck Simulations of Current Flow through Gramicidin-A.
J. Comp. Science (2001), 373-409.[26] 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.[27] S. Ji, B. Eisenberg, and W. Liu, Flux Ratios and Channel Structures.
J. Dynam. DifferentialEquations (2017). https://doi.org/10.1007/s10884-017-9607-1.[28] 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.[29] 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.[30] G. Lin, W. Liu, Y. Yi, and M. Zhang, Poisson-Nernst-Planck systems for ion flow with a localhard-sphere potential for ion size effects. SIAM J. Appl. Dyn. Syst. (2013), 1613-1648.[31] W. Liu, Geometric singular perturbation approach to steady-state Poisson-Nernst-Planck systems. SIAM J. Appl. Math. (2005), 754-766.[32] W. Liu, One-dimensional steady-state Poisson-Nernst-Planck systems for ion channels with multipleion species. J. Differential Equations (2009), 428-451.[33] W. Liu, X. Tu, and M. Zhang, Poisson-Nernst-Planck systems for ion flow with density functionaltheory for hard-sphere potential: I-V relations and critical potentials. Part II: Numerics.
J. Dynam.Differential Equations (2012), 985-1004.
34] W. Liu and B. Wang, Poisson-Nernst-Planck systems for narrow tubular-like membrane channels.
J. Dynam. Differential Equations (2010), 413-437.[35] W. Liu and H. Xu, A complete analysis of a classical Poisson-Nernst-Planck model for ionic flow. J. Differential Equations (2015), 1192-1228.[36] N. F. Mott, The theory of crystal rectifiers.
Proc. Roy. Soc. A (1939), 27-38.[37] 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.[38] W. Nonner, L. Catacuzzeno, and B. Eisenberg, Binding and selectivity in L-type Calcium channels:A mean spherical approximation. Biophys. J. (2000), 1976-1992.[39] 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.[40] 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.[41] 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.